Sergiy Bogomolov

Christian Schilling

Ezio Bartocci

Gregory Batt

Hui Kong

Radu Grosu

Abstraction-based Parameter Synthesis for Multiaffine Systems

Multiaffine hybrid automata (MHA) represent a powerful formalism to model complex dynamical systems. This formalism is particularly suited for the representation of biological systems which often exhibit highly non-linear behavior. In this paper, we consider the problem of parameter identification for MHA. We present an abstraction of MHA based on linear hybrid automata, which can be analyzed by the SpaceEx model checker. This abstraction enables a precise handling of time-dependent properties. We demonstrate the potential of our approach on a model of a genetic regulatory network and a myocyte model.

Introduction

Hybrid automata can model systems from a wide range of real-world domains. Due to its behavioral complexity, the biological domain can particularly benefit from the expressiveness of hybrid automata [START_REF] Bartocci | Modeling and simulation of cardiac tissue using hybrid I/O automata[END_REF]. However, biological models mostly have highly non-linear dynamics.

Parameter identification is the problem where we want to find a parameter set for which a given property is satisfied by the system. In the biological domain, this problem is of large importance considering the current limitations on experimental measurement techniques [START_REF] Myers | Engineering genetic circuits[END_REF].

In this paper, we present a novel approach to solve the parameter identification problem for the class of multiaffine hybrid automata (MHA). We reduce the parameter identification problem to solving multiple verification problems. In short, the algorithm consists of the following steps: We partition the parameter space into a number of equivalence classes. Given an equivalence class, we show how the system behavior can be approximated with a linear hybrid automaton (LHA), which can be analyzed by the hybrid model checker SpaceEx [START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF]. In addition, we utilize a hierarchical search to start the analysis with coarser regions and iteratively refine the partition based on the model structure. We are also able to prune the search when we detect that our analysis will not find any parameters in a subregion. We have implemented our approach and show its potential on a genetic regulatory network and a myocyte model.

Outline. The rest of the paper is organized as follows. In Section 2, we introduce some preliminary notions. Then we present our new approach, first the construction of the relevant parts in Section 3, followed by the hierarchical search procedure in Section 4. In Section 5, we evaluate the approach on two biological models. We discuss related work in Section 6 and conclude in Section 7.

Preliminaries

In this section, we introduce the notions used in the rest of the paper.

Multiaffine function. A multiaffine function f : R n → R q (n, q ∈ N) is a polynomial in the variables x 1 , . . . , x n with the property that the degree of f in any of the variables is less than or equal to 1 [START_REF] Kloetzer | Reachability analysis of multi-affine systems[END_REF]. Formally, f has the following form:

f (x 1 , . . . , x n ) = i1,...,in∈{0,1} c i1,...,in x i1 1 • • • x in n ,
with c i1,...,in ∈ R q for all i 1 , . . . , i n ∈ {0, 1} and the convention that x 0 k = 1. Hybrid automaton. A hybrid automaton (HA) [START_REF] Alur | The algorithmic analysis of hybrid systems[END_REF] is a mathematical model with both continuous and discrete behavior. It is represented by the tuple H = (Loc, Var , Inv , Flow , Trans, Init). Loc is a set of discrete locations. Var is a set of real-valued variables x 1 , . . . , x n . Each ∈ Loc is associated with a set of differential equations (or inclusions) Flow ( ) that defines the time-driven evolution of the continuous variables. A state s ∈ Loc × R n consists of a location and values of the continuous variables x 1 , . . . , x n . The set of discrete transitions Trans defines how the state can jump between locations when inside the transition's guard set. The system can remain in a location while the state is inside the invariant set Inv ( ). All behavior originates from the set of initial states Init. A trajectory is a function which defines the state of the HA for every time moment. In the verification setting, we are interested in whether there exists a trajectory from the set Init to a set Bad which defines the bad states to be avoided.

Let x(t) ∈ R n denote the values of the continuous variables at time t. We consider continuous dynamics Flow of the following two forms. If ẋ(t) = f (x, t) where f (x, t) is a multiaffine function, then the HA is called a multiaffine hybrid automaton (MHA). If ẋ(t) ∈ P where P is a polytope, then the HA is called a linear hybrid automaton (LHA). We always consider convex polytopes and omit the dependence of f on t in what follows.

Genetic regulatory network.

A genetic regulatory network [START_REF] Batt | Temporal logic analysis of gene networks under parameter uncertainty[END_REF] is defined by the dynamics of the following form:

ẋi = f i (x, p) = j∈Pi κ ij r P ij (x) - j∈Di γ ij r D ij (x) x i , i = 1, . . . , n (1)
r the piecewise-affine (PA) models proposed in [START_REF] Grosu | From cardiac cells to genetic regulatory networks[END_REF] (see also [START_REF] Dang | Sensitive state-space exploration[END_REF]). However, contrary to the step functions used in PA models, ramp functions capture the graded response of gene expression to continuous changes in effector concentrations.

-(x i , θ i , θ i ) r + (x i , θ i , θ i ) θ i x i 0 θ i 1 θ i x i 0 θ i 1 (a) A B b a (b) ẋa = κar -(xb, θ 1 b , θ 2 b ) -γa xa ẋb = κbr -(xa, θ 1 a , θ 2 a ) -γb xb (c) φ1 = (xa < θ 1 a ∧ xb > θ 2 b → G (xa < θ 1 a ∧ xb > θ 2 b )) ∧(xb < θ 1 b ∧ xa > θ 2 a → G (xb < θ 1 b ∧ xa > θ 2 a )) (e) κi γi θ 1 i θ 2 i xa [0, 40] 2 8 12 xb [0, 20] 1 8 12 (d)

Embedding transition systems

The specific form of the PMA functions f suggests a division of the state space X into hyperrectangular regions (Figure 2(a) for our example). For every i ∈ {1, . . . , n}, let Λ i = {λ j i } j∈{1,...,li} be the ordered set of all threshold constants in f , and of all atomic proposition constants in Π, associated with gene i, together with 0 and max xi . The cardinality of Λ i is l i . Then, we define R as the following set of n-dimensional hyperrectangular polytopes R ⊆ X , simply called rectangles:

R = {R c | c = (c 1 , . . . , c n ) and ∀i ∈ {1, . . . , n} : c i ∈ {1, . . . , l i -1}}, where R c = {x ∈ X | ∀i ∈ {1, . . . , n} : λ ci i < x i < λ ci+1 i }.
c is the coordinate of the rectangle R c . The union of all rectangles in X is denoted by X R : X R = ∪ R∈R R. Note that X R = X . Notably, threshold hyperplanes are not included in X R . Two rectangles R and R , are said adjacent, denoted R R , if they share a facet. coord : R → n i=1 {1, . . . , l i -1} maps every rectangle R ∈ R to its coordinate, and rect : X R → R maps every point x in X R to the rectangle R such that x ∈ R. For the cross-inhibition network, the set R = {R 11 , . . . , R 33 } of all rectangles is represented in Figure 2(b). R 11 and R 21 are adjacent, whereas R 11 and R 22 are not.

Formally, the semantics of a PMA system Σ is defined by means of an embedding transition system. Here x i is the i-th component of the state vector x ∈ X ⊂ R n . P i and D i are sets of indices. κ ij and γ ij are production and degradation rate parameters, respectively. We assume that some parameters are uncertain, i.e., are defined on proper intervals. We denote the number of uncertain parameters by m. Therefore, the system is parametrized by the vector p = (p 1 , ..., p m ) ∈ D, where D is the hyper-rectangular domain of uncertain parameters.

The terms r ij are continuous piecewise-multiaffine functions arising from products of ramp functions r + and r -of the form shown in Fig. 1(a). Each r ij captures the combined impact of several regulatory proteins in the sets P i and D i , respectively, on the control of the production or degradation of protein i. Assuming protein i does not regulate its own degradation, i.e., x i does not occur in r D ij (x) for j ∈ D i , function f = (f 1 , . . . , f n ) is multiaffine in x and affine in p.

We note that the ramp functions induce a partition of the state space X into a grid of hyper-rectangular regions H. The values of the separating hyperplanes are called thresholds θ. Let θi be the number of thresholds in dimension i. Definition 1. Let H := H c | c = (c 1 , . . . , c n ), c i ∈ {1, . . . , θi -1}, i = 1, . . . , n be the set of hyper-rectangles H c with coordinates c in the grid. Furthermore, let coord : H → n i=1 {1, . . . , θi -1} map hyper-rectangle H c to its coordinate c. Problem statement. A property ϕ specifies the desired system behavior. In this paper, we consider properties of the form R Init → ¬♦R Bad where the region R Init denotes the initial system states and the region R Bad denotes the states to be avoided. Note that ϕ belongs to the class of safety properties. Our goal is to identify a subset of the parameter domain D which ensures that the property ϕ holds for a given MHA, i.e., R Bad is avoided when starting in R Init .

Note that an MHA provides a semantically equivalent representation of the dynamics [START_REF] Alur | The algorithmic analysis of hybrid systems[END_REF]. We consider an abstraction of a parametric MHA M (p) to an infinite transition system. Definition 2. Let H x be the hyper-rectangle strictly containing x. Given a parameter p and an MHA M (p), the embedding transition system T M (p) is defined as follows. The states X are the same as the continuous states of M (p). There is a transition x → x iff: 1) H x and H x are either equal or adjacent. 2) There is a solution ξ of (1) and time points t 0 < t 1 such that ξ(t 0 ) = x, ξ(t 1 ) = x and for all t in [t 0 , t 1 ], ξ(t) stays within H x or H x . Almost all trajectories in M (p) are represented in T M (p). The exception are trajectories not passing through a common facet of two hyper-rectangles. In what follows, we consider the representation of the system T M (p).

Example (two-genes network). In the following, we illustrate our approach on the two-genes network [START_REF] Batt | Robustness analysis and tuning of synthetic gene networks[END_REF] (also called toggle switch or cross-inhibition network ).

ẋa = κ a • r -(x a , θ 4 a , θ 5 a ) • r -(x b , θ 2 b , θ 3 b ) -γ a x a κ a ∈ [0, 30], γ a = 1 ẋb = κ b • r -(x a , θ 2 a , θ 3 a ) -γ b x b κ b ∈ [0, 40], γ b = 2 (θ 1 a , θ 2 a , θ 3 a , θ 4 a , θ 5 a , θ 6
a ) = (0, 8, 12, 18, 22, 30) (θ 1 b , θ 2 b , θ 3 b , θ 4 b ) = (0, 8, 12, 20) Here x a and x b define the concentrations of the proteins A and B, respectively. The uncertain parameters κ a and κ b define the range of their production rates in the given intervals. As Fig. 1(b) shows, protein A inhibits the production of both proteins A and B, while protein B only inhibits the production of protein A. We are interested in checking whether the protein concentrations cannot reach some specific threshold values when starting in a given initial region.

In addition, we consider an extended version of the dynamics (1) which features a stimulus. The stimulus is a time-dependent function which models an external influence on the system.

Example (two-genes network with stimulus). We extend the previous example with the first equation now featuring some stimulus u:

ẋa = κ a • r -(x a , θ 4 a , θ 5 a ) • (1 -r + (x b , θ 2 b , θ 3 b ) • (1 -u)) -x a u = r -(t, t 2 , t 3 ) (t 1 , t 2 , t 3 , t 4 ) = (0, 0.29, 0.3, 1)
We use a stimulus which is 1 at the beginning up until 0.29 ms and then drops linearly to 0 within 0.01 ms, expressed by the ramp function of time r -(t, 0.29, 0.3). This stimulus regulates the production of the protein A together with the protein B. The term (1

-r + (x b , θ 1 b , θ 2 b ) • (1 -u))
encodes the logical formula ¬(x b ∧¬u), which is equivalent to ¬x b ∨u. Thus, this term contributes to the production of the protein A whenever the protein B is absent or the stimulus is present. Since the stimulus is time-dependent and decreasing, this means that the inhibitory effect of the protein B is only relevant for the production of the protein A after 0.29 ms.

Abstraction of MHA

In this section, we first introduce an LHA L M (p) which overapproximates the behavior of the transition system T M (p) for a particular parameter value p ∈ D. In order to provide efficient exploration of the parameter state space, we then lift this definition to parameter sets. We use the two-genes model as a running example throughout this section.

Pointwise LHA abstraction

State space and invariants. We use the same continuous variables for L M (p) as in the MHA M (p). Thresholds θ partition the continuous state space into a grid of hyper-rectangular regions, which naturally induces a discrete structure of the LHA and location invariants. In particular, we map every hyper-rectangular region to a location in the LHA and use the bounds on the regions as invariants. Example. In the two-genes example we get a 2D-partition by the planes at x a = θ i a and x b = θ j b . Parts of the state space and the associated locations are shown in Fig. 2

x b xa H1 H2 H3 H4 H1 H2 H3 H4 θ 1 a θ 2 a θ 3 a θ 1 b θ 2 b θ 3 b (a) x b xa κ b κa H P κ 1 a κ 2 a κ 1 b κ 2 b θ 1 a θ 2 a θ 1 b θ 2 b d1 d2 d3 d4 (b) x b xa Q(d 1 ) Q(d 2 ) Q(d 4 ) Q(d 3 ) Q ∀ (P ) Q ∃ (P ) κ 1 a - θ 2 a κ 2 a - θ 2 a κ 1 a - θ 1 a κ 2 a - θ 1 a κ 1 b -2θ 2 b κ 2 b -2θ 2 b κ 1 b -2θ 1 b κ 2 b -2θ 1 b (c)

(a). The invariant of location

H 3 is θ 1 a ≤ x 1 ≤ θ 2 a ∧ θ 2 b ≤ x 2 ≤ θ 3 b . Discrete transitions.
We use the quotient of T M (p) with respect to the state space partition to define the discrete transitions. Definition 3. Let H and H be the regions associated with the location and , respectively. L M (p) has a transition from location to if and are adjacent, there is a solution ξ of (1) and time points t 0 < t 1 < t 2 such that ξ(t) ∈ Inv ( ) for all t in [t 0 , t 1 ], and ξ(t) ∈ Inv ( ) for all t in [t 1 , t 2 ].

We do not add any guards on the transitions as the chosen invariants already ensure that a transition between two locations can only be taken on the common facet of two adjacent regions H and H .

In order to effectively construct the transitions, we use the facts that the dynamics f (x, p) are multiaffine in x and we consider only hyper-rectangular regions for the locations in L M (p). In the following, let hull denote the convex hull operator. Theorem 1. [START_REF] Belta | Controlling a class of nonlinear systems on rectangles[END_REF] Let f : R n → R n be a multiaffine function and H ⊂ R n be a hyper-rectangle with corner set C H . Then

f (H) ⊆ hull ({f (v) | v ∈ C H }).
Intuitively, this theorem says that the behavior of f inside a hyper-rectangle H is completely determined by the behavior of f in the corners of H. As a consequence, the following proposition along the lines of a similar proposition for Kripke structures [START_REF] Batt | Temporal logic analysis of gene networks under parameter uncertainty[END_REF] can be proven. Proposition 1. L M (p) has a transition from location to associated with hyper-rectangles H and H only if the projection of f (x, p) on the H → H direction is positive in at least one corner of the facet separating H from H .

Direction and strength of the derivative ẋi = f i (x, p) in a corner v of a hyperrectangle depends linearly on parameter vector p. As a consequence, f i (v, p) = 0 is the hyperplane separating parameter values p where ẋi is positive from the ones where ẋi is negative. Thus, Proposition 1 allows us to construct the transitions of L M (p) based on the sign of the function f at the vertices of the hyper-rectangles.

Example. The transitions for the excerpt shown in Fig. 2(a) are determined by the direction of the derivatives in the corners (shown in blue). For instance, we add a transition from H 1 to H 2 because there is a corner, e.g., (θ 2 a , θ 1 b ), which point to this direction, but there is no transition from H 2 to H 1 .

Continuous flows. For computing the flows of the LHA we again use the multiaffine dependence of f (x, p) on x and the affine dependence on p.

For a fixed hyper-rectangle H with corner set C H and a fixed parameter vector p, by Theorem 1 we know that f (x, p) is included in the convex hull of the hyper-rectangle corners v ∈ C H . Therefore, we can bound the flow of L M (p) by a polytope, i.e., the dynamics can be represented in the form of differential inclusion.

Definition 4. The flow of L M (p) is defined as Q(p) := hull ({f (v, p) | v ∈ C H }).

Set-based LHA abstraction

In order to handle infinite sets of parameters, we lift the pointwise definition of L M (p) to sets of parameters. In particular, given a parameter polytope P , we introduce an LHA L ∃ M (P ) which overapproximates the behavior of L M (p) for all p ∈ P . Therefore, if L ∃ M (P ) satisfies a property ϕ, we can conclude that P is a valid parameter set. Otherwise, we partition the parameter set P into two subsets P 1 and P 2 and proceed with their analysis. In order to prune the parts of the parameter space where our analysis will not provide any valid parameters, we introduce a further LHA called L ∀ M (P ) which underapproximates the LHA L M (p) for the parameter class P . In the following, we assume a parameter p ∈ P .

State space and invariants. We use the same state space and invariants for both L ∃ M (P ) and L ∀ M (P ) as defined for L M (p). Discrete transitions. The theorem below provides an effective way to compute the image of f (v, p) for a particular state space corner v and a parameter vector p ∈ P . Theorem 2. [START_REF] Habets | Reachability and control synthesis for piecewise-affine hybrid systems on simplices[END_REF] Let f : R m → R m be an affine function and P ⊂ R m be a convex polytope with corner set C P . Then

f (P ) = hull ({f (d) | d ∈ C P }).

Now, L ∃

M (P ) has a transition from location to if there is a transition from to in L M (p) for some p ∈ P . Analogously, L ∀ M (P ) has a transition from location to if there is a transition from to in L M (p) for all p ∈ P . Definition 5. Let c := coord (H) and c := coord (H ) be the coordinates of two adjacent hyper-rectangles. Furthermore, let V := C H ∩ C H be the corners on the separating facet and let and be the locations associated to H and H , respectively.

We define g( ,

) := v∈V {p ∈ P | f i (v, p) • (c i -c i ) > 0}
, where i ∈ {1, . . . , n} such that c i -c i = 0. A transition → belongs to the LHA -L ∃ M (P ) if g( , ) = ∅, and to -L ∀ M (P ) if g( , ) = P . Note that Theorem 2 allows us to construct the parameter set satisfying the constraint f i (v, p) • (c i -c i ) > 0 by only considering the vertices of P . The term c i -c i = ±1 is used to express the direction of the transition. The construction uses the union operation and test for equality and emptiness for polytopes.

Continuous flows. For L ∃ M (P ) to be an overapproximation of all L M (p) and L ∀ M (P ) to be an underapproximation, the tightest definition we can find is the union and intersection of all Q(p), respectively. Let Q * ∃ (P ) := p∈P Q(p) and Q * ∀ (P ) := p∈P Q(p). Computing Q * ∃ (P ) and Q * ∀ (P ) is, however, infeasible. Therefore, similarly to the transition construction, we propose an approximation which relies on the values of the derivatives in the corners of state space partitions H and parameter sets P .

Theorem 3. Let f (x, p) : R n × R m → R n with x ∈ R n
and p ∈ R m be a multiaffine function which is affine in p, H ⊂ R n be a hyper-rectangle with a corner set C H and P ⊂ R m be a convex polytope with a corner set C P . Then the following holds:

-p∈P Q(p) ⊆ hull d∈C P Q(d) , -p∈P Q(p) ⊆ d∈C P Q(d).
Note that the left-hand sides are Q * ∃ (P ) and Q * ∀ (P ), respectively. Based on this theorem, we define the flows in the following way. Definition 6. The flow of L ∃ M (P ) is defined as Q ∃ (P ) := hull d∈C P Q(d) . The flow of L ∀ M (P ) is defined as Q ∀ (P ) := d∈C P Q(d). We obtain an algorithm for computing Q ∃ (P ) and Q ∀ (P ) by first traversing all vertices v ∈ C H and d ∈ C P and collecting f (v, d) in order to compute the polytope Q(d). In the end, we take the finite union and intersection of those polytopes. Note that similar to L M (p) we end up with LHA whose dynamics are defined by differential inclusions.

The following proposition relates the flow representations we have introduced.

Proposition 2. The sets Q(p), Q * ∃ (P ), Q * ∀ (P ), Q ∃ (P ) and Q ∀ (P ) are related as follows.

-

Q * ∀ (P ) ⊆ Q ∀ (P ). -Q * ∀ (P ) ⊆ Q(p) ⊆ Q * ∃ (P ) ⊆ Q ∃ (P ) for all p ∈ P .
Thus, Q ∃ (P ) is indeed an overapproximation of L M (p) as required. However, while Q * ∀ (P ) is an underapproximation of L M (p), Q ∀ (P ) is not necessarily an underapproximation of L M (p) for p ∈ P . We discuss this issue in the next section.

We define the automaton L ∃ * M (P ) by replacing Q ∃ (P ) with Q * ∃ (P ) in the continuous flow of the automaton L ∃ M (P ). We derive L ∀ * M (P ) from L ∀ M (P ) in the analogous way, i.e., by replacing Q ∀ (P ) with Q * ∀ (P ). Example. Consider the state space rectangle H and the parameter space rectangle P in Fig. 2(b). Recall that the state equation ẋ = f (x, p) is given as

ẋa = f a (x, p) = κ a -x a ẋb = f b (x, p) = κ b -2x b .
Hence the dynamics f (v, d), for v and d in the corner sets C H and C P of H and P , respectively, are of the form (κ i a -θ j a , κ k b -2θ b ). Let the corners of the parameter space rectangle P be denoted in anticlockwise order as d 1 , d 2 , d 3 and d 4 . Now construct the state space rectangles Q(d 1 ), Q(d 2 ), Q(d 3 ) and Q(d 4 ). The intersection of all these rectangles results in the rectangle Q ∀ (P ), while the union is the rectangle Q ∃ (P ). Since it is already convex, hull is the identity operation in this case. The results are visualized in Fig. 2(c).

We observe that, by construction, L M (p) is a conservative abstraction of T M (p), and L ∃ M (P ) is a conservative abstraction of L M (p), i.e., if L ∃ M (P ) satisfies a safety property, then so does T M (p). This fact ensures the soundness of our approach. Proposition 3. L ∃ M (P ) is an overapproximation of T M (p) for any p ∈ P .

Hierarchical parameter search

In this section, we first show that by using sampling techniques we can compute the automaton L ∀ * M (P ) with arbitrary precision. Afterwards, in order to leverage different levels of abstractions during the parameter space exploration, we introduce a discrete abstraction of MHA. Finally, we describe an abstractionbased parameter search procedure which explores the parameter domain in a hierarchical fashion.

Computation of underapproximative abstractions

As outlined in the previous section, the role of the underapproximation L ∀ * M (P ) is to detect parameter regions which our approach cannot classify as valid ones. The LHA L ∀ M (P ), an overapproximated version of L ∀ * M (P ), does not generally underapproximate L M (p) for all p ∈ P . The reason is that the flows Q ∀ (P ) do not necessarily underapproximate the flows Q(p) of all L M (p). In particular, as we take an intersection only over corner points of a considered parameter set, there might exist some p ∈ P such that Q ∀ Note that as we only use L ∀ M (P ) for pruning purposes, soundness of our approach is not affected by imprecision. We might at most ignore some parameter region which could have been classified as valid by our approach. Still, in order to improve the precision of L ∀ M (P ), we can randomly sample parameter vectors p ∈ P and consider the intersection of Q(p) with Q ∀ (P ).

(P ) ⊆ Q(p). L M (p) L ∀ * M (P ) L ∀ M (P ) L ∃ * M (P ) L ∃ M (P ) K ∀ M (P ) K M (p) K ∃ M (P ) T M (p)
The following theorem describes the possible improvements by sampling. For notational convenience, let • denote the Euclidean norm, ∂S denote the border of a closed set S, d (t, S) := min s∈S t -s denote the distance of t to the border of S, and Ker (S) := {s ∈ S | d (s, ∂S) ≥ }.

Theorem 4. Let f (x, p) : R n × R m → R n with x ∈ R n and p ∈ R m be a multiaffine function which is affine in p, H ⊂ R n be a hyper-rectangle and P ⊂ R m be a convex polytope. Then the following formulae hold:

∀p, p ∈ P. lim p-p →0 Q(p) = Q(p ) (2)
∀p ∈ P , > 0. ∃δ > 0. ∀p ∈ P. p -p < δ

=⇒ Ker 2 (Q(p)) ⊆ Q(p) ∩ Q(p ) ∧ Ker 2 (Q(p )) ⊆ Q(p) ∩ Q(p ) (3)
This theorem asserts that if we sample sufficiently many points in a uniform way in parameter space, then we can approximate L ∀ * M (P ) infinitely closely for non-degenerate cases, i.e., when Ker 2 (Q(p)) = ∅ for all p ∈ P . However, for practical purposes, sampling can clearly shrink the overapproximation even in degenerate cases.

Fig. 3 shows the relationship between the systems considered in this paper.

Discrete abstraction of MHA

For later discussion, it is useful to define the induced Kripke structures (KS) of the LHA L M (p), L ∃ M (P ) and L ∀ M (P ). Basically, we drop the continuous behavior and map initial and bad states to the locations with non-empty intersection. We denote by K M (p), K ∃ M (P ), K ∀ M (P ) the Kripke structures induced by L M (p), L ∃ M (P ), L ∀ M (P ), respectively. Clearly, the induced Kripke structure is a conservative abstraction of the LHA as it allows for additional trajectories to the bad states for two reasons. The behavior in the states is unconstrained due to the absence of flows, and the initial and bad states of the Kripke structure overapproximate their LHA counterparts.

Proposition 4. K M (p) (K ∃ M (P ), K ∀ M (P ), respectively) is an overapproxima- tion of L M (p) (L ∃ M (P ), L ∀ M (P )
, respectively) for any p ∈ D (P ⊆ D, respectively).

We incorporate the Kripke structures into our approach in the following way. Whenever we analyze an LHA, we first analyze the respective KS. If the KS satisfies the given property, we skip the LHA analysis. This is justified because by Proposition 4 we know that the respective LHA will also satisfy the property. In this way, we improve analysis performance as the reachability problem for LHA is computationally harder to solve. We note that the construction of the KS does not impose any further computational efforts since we need to construct the locations and transitions for the LHA anyway.

Parameter identification

Given an MHA M with ẋ = f (x, p), a safety property ϕ and the domain of uncertain parameters D, we explore the parameter space in a hierarchical way. Recall from the construction of the transitions that the constraints f (x, p) = 0 are the separating hyperplanes responsible for adding transitions. Based on those constraints, the instantiation of f (x, p) in every corner v of the state space leads to a parameter space partition into polytopes.

We now explain the algorithm with the help of the pseudocode given in Fig. 4.

In a preprocessing step, the algorithm examines the corners of the state space partition and collects the constraints Ψ over the parameters in the function CollectConstraintsList (line 3) such that f (x, p) = 0.

Next, the algorithm moves on to the function Explore (line 5) which actually implements the search in the parameter space. This function successively builds a number of abstractions of the MHA for a considered parameter set in order to find valid subsets. It takes a list CL of constraints which encodes hyperplanes used to define the current parameter set. We initially call the function with CL = ε as we first consider the whole parameter space. Now we look at the function Explore in more detail.

We start by calling the function PolytopeFromConstraintsList (line 7) which builds a parameter polytope P based on the provided list CL of constraints over parameters and the parameter space domain D. In line 8 we compute the KS K ∃ M (P ), K ∀ M (P ) and the LHA L ∃ M (P ), L ∀ M (P ). Note that on the implementation level we compute them only on demand. The computed approximations are analyzed in the following way: e l s e i f ( ¬Reach ( L ∃ M , ϕ ) )

V := V ∪ P ; return ; % v a l i d s e t found by t h e LHA a n a l y s i s

13 e l s e i f ( ¬b ∨ Reach ( K ∀ M , ϕ ) ) i f ( Reach ( L ∀ M , ϕ ) ) 16
return ; % no f u t u r e v a l i d s e t s f o r t h e LHA a n a l y s i s 1. If the property ϕ holds for the KS K ∃ M (P ) already (line 9), we conclude that the current parameter set P is valid. Therefore, we add P to the set of valid parameters V and stop considering the current branch in the search tree (line 10). 2. If the discrete abstraction K ∃ M (P ) was too coarse to prove the validity of P , we continue with the finer analysis using L ∃ M (P ) (line 11). Similar to step 1, in the case of property satisfaction we add P to the valid parameters (line 12). 3. If the parameter set validity has not been shown up to now, we proceed to the pruning phase by considering K ∀ M (P ) and L ∀ M (P ). If both of them violate the property ϕ, we prune the search tree as we expect that no valid parameter sets can be found for any subset of P . Note that due to efficiency reasons we first analyze K ∀ M (P ) (line 13) and move on to L ∀ M (P ) (line 15) only if the KS is not safe with respect to the property ϕ. If L ∀ M (P ) is not safe either, we prune the current subtree. However, if L ∀ M (P ) is safe, we continue with the search. In this case, we can omit the KS analysis for all nodes in the current subtree as it will always give the same result. We assume the conditions in lines 9 and 13 are evaluated in a lazy fashion and therefore we only use the KS analysis based on the value of the Boolean switch b. . . . . . .

4.

If K ∀ M (P ) or L ∀ M (P ) are safe, we partition the parameter set P into two subsets by considering a further constraint from the list Ψ (line 18). Those two subsets correspond to the positive and negative values of the chosen constraint, respectively. We proceed by recursively analyzing both subsets (lines 19-20).

Search tree implications. In Fig. 5, we illustrate the potential impact of the LHA L ∃ M (P ) and L ∀ M (P ) on the structure of the search tree compared to an approach only using the KS analysis. We observe that L ∃ M (P ) satisfies the property ϕ in the node P 1 , whereas K ∃ M (P ) violates the property. The LHA analysis benefits in two ways from this result. Firstly, it finds a large parameter set P 1 , whereas the KS analysis can at most find valid sets in some of the child nodes, which are subsets of P 1 . Secondly, the LHA analysis does not explore the children of P 1 in the search tree, which improves the algorithm performance. Moreover, L ∀ M (P ) satisfies the property ϕ in the node P 0 , whereas K ∀ M (P ) violates the property. Therefore, the KS analysis prunes the subtree P 0 , but a valid parameter set P 00 can be found by the LHA analysis.

Evaluation

We have implemented the algorithm in MATLAB in the tool Hydentify. We use the library PPL [START_REF] Bagnara | The parma polyhedra library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF] for the operations on polytopes and the SpaceEx model checker [START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF] for the analysis of the LHA. Note that the default version of SpaceEx does not stop immediately after having found a property violation. Therefore, we have modified the version of SpaceEx so that it stops as soon as a property violation has been detected. This adjustment lets us improve the analysis performance. We apply the PHAVer scenario of SpaceEx which uses constraint polytopes to represent reachable regions to precisely analyze LHA.

Batt et al. [START_REF] Batt | Temporal logic analysis of gene networks under parameter uncertainty[END_REF] have presented a parameter identification approach for multiaffine systems implemented in a tool called RoVerGeNe. They approximate the system on the level of the induced Kripke structures. In the following evaluation, we compare the parameter identification results of our approach and RoVerGeNe. The integration of the induced KS into our algorithm allows for both a qualitative and a quantitative comparison of the two approaches. The implementation and models we used for the evaluation are available online Two-genes network model. We evaluate our tool on a number of models from the class of genetic regulatory networks. The experiments have been performed on a notebook with an Intel Core 2 Duo @ 2.26 GHz processor and 4 GB RAM.

For the evaluation purpose, we consider two classes of the two-genes network model introduced in Section 2. The first model class is the original system, while the second class is augmented by a stimulus. Note that our LHA framework enables an easy modeling and analysis of models with time dependent stimuli.

In particular, we model a stimulus as a ramp function of an auxiliary variable t defined by the differential equation ṫ = 1. For every model class, we present two model instances. We look for parameters which lead to the repression of a given protein. For every instance and parameter identification algorithm, we report the following data: the coverage of the parameter domain, the number of the valid parameter sets found, the number of nodes in the search tree considered, the number of KS and LHA analyzed, and the runtime in seconds. By the term parameter coverage we denote the relation of the volume of the found valid parameters to the volume of the whole parameter domain D. The results are provided in Table 1. The instances 1-2 correspond to the model class without a stimulus, whereas the other two instances belong to the class with a stimulus.

We first observe that the valid parameter regions found by our algorithm are usually much larger than the ones found by RoVerGeNe for both the models with and without the stimulus. Instance 2 provides a particularly illustrative example for the difference. Here, RoVerGeNe does not find any valid parameters, whereas our approach discovers valid parameter regions covering 38 % of the whole parameter domain. This behavior can be justified as follows. On the one hand, RoVerGeNe reports that both K ∃ M (P ) and K ∀ M (P ) at the root level reach the bad states. This results in analysis termination of RoVerGeNe. On the other hand, our algorithm proceeds in-depth with the analysis of the parameter space and detects 5 valid parameter regions. In instance 3, we see similar impact of taking LHA into account. In particular, both approaches consider 5 K ∀ M (P ). However, our approach additionally considers 3 L ∀ M (P ) which allow to extra unfold the parameter space. In this way, our approach analyzes 15 nodes and finds 5 valid regions compared to 9 nodes and 3 valid regions for RoVerGeNe, respectively. At the same time, the refined precision of LHA can shrink the search space. For example, in instance 4, the new algorithm achieves the parameter coverage of 84 % vs. 60 % by RoVerGeNe having considered only 9 nodes vs. 13 nodes in case of RoVerGeNe. We note that the valid parameter sets which are near the search tree root lead to larger parameter coverage with only a few parameter sets. This fact is confirmed by instance 4 where our approach finds 3 valid sets which cover a bigger region than the 4 valid sets found by RoVerGeNe.

Myocyte model. A fundamental question in the treatment of cardiac disorders, such as tachycardia and fibrillation [START_REF] Cherry | Visualization of spiral and scroll waves in simulated and experimental cardiac tissue[END_REF], is the identification of circumstances under which such a disorder arises. Cardiac contraction is electrically regulated by particular cells, known as myocytes. For each electric stimulus originating in the sino-atrial node of the heart (its natural pacemaking unit), the myocytes propagate this stimulus and enforce the contraction of the cardiac muscle, known as a heart beat. Grosu et al. [START_REF] Grosu | From cardiac cells to genetic regulatory networks[END_REF] have identified an MHA model for human ventricular myocytes and recast the biological investigation of lack of excitability to a computational investigation of the parameter ranges for which the MHA accurately reproduces lack of excitability. We apply our algorithm to this model and compare its performance with RoVerGeNe. The model has 4 continuous variables and 4 parameters. In our setting, a valid parameter set ensures that the myocyte is not excited.

We remark that our parameter identification approach has a large potential with respect to parallelization as the LHA and KS can be analyzed independently. We made use of this property and utilized a parallel version of our implementation for the analysis of the myocyte model. The experiments have been run on a Linux cluster with 32 AMD @ 2.3 GHz cores and 256 GB RAM. The model behavior is analyzed within a biologically reasonable time span of 1 ms. We note that the stimulus and particularly its duration require a special treatment as it strongly impacts the myocyte behavior. The stimulus in our model starts with the value 1 and linearly drops to 0. We explore the impact of the stimulus length on the myocytes excitement.

Our approach empirically shows that the whole parameter domain is valid for all stimuli of length up to approximately 0.12 ms. In other words, we can provide a lower bound on the stimulus length which makes the myocyte model excitable. For this purpose, we have discretized the stimulus length with a step of 0.1 ms, i.e., we have considered stimuli of the length 0, 0.1, . . . , 1 ms. Having identified an interval of interest [0.1; 0.2], we have discretized it in a finer way with a step of 0.02 ms. The new analysis takes 187 seconds and detects that the whole parameter domain is valid for the stimulus of length 0.12 ms, whereas RoVerGeNe reports the coverage of 29 % after 48 seconds. The parameter coverage computed by our algorithm drops to 30 % for the stimulus length of 0.14 ms and the analysis takes 1785 seconds. We note that the coverage computed by RoVerGeNe stays the same for all stimulus lengths as it cannot reason about time. This is a conceptual improvement over RoVerGeNe.

A number of approaches have been developed to solve the parameter identification problem for hybrid automata. First, as already outlined in the previous section, Batt et al. [START_REF] Batt | Temporal logic analysis of gene networks under parameter uncertainty[END_REF] presented a parameter identification approach based on the abstraction of MHA by Kripke structures. By using our LHA abstraction, we improve the abstraction precision and in this way find more valid parameters. Dang et al. [START_REF] Dang | Sensitive state-space exploration[END_REF] introduced a "sensitive barbarian" approach. Bartocci et al. [START_REF] Bartocci | A temporal logic approach to modular design of synthetic biological circuits[END_REF] consider a modular version of this approach. The main idea is to combine numerical simulation with sensitivity analysis to reduce the considered parameter space. A crucial difference to our approach lies in the fact that we utilize a symbolic analysis of the reachable states. In a further approach, Dreossi et al. [START_REF] Dreossi | Parameter synthesis for polynomial biological models[END_REF] provide a parameter synthesis algorithm for polynomial dynamical systems. Their synthesis technique uses the Bernstein polynomial representation and recasts the synthesis problem as a linear programming problem. Note that they consider only discrete time dynamical systems, whereas we treat time as a continuous entity. The work by Liu et al. [START_REF] Liu | Parameter synthesis for cardiac cell hybrid models using delta-decisions[END_REF] tackles the parameter synthesis problem using δ-complete decision procedures [START_REF] Gao | δ-complete decision procedures for satisfiability over the reals[END_REF] for first-order logic (FOL) formulae to overcome undecidability issues. In this setting, a FOL formula describes the states reachable with a finite number of steps. Therefore, the parameter identification problem is reduced to finding a satisfying valuation of the parameters for this formula. This approach requires enumerating all the discrete paths of a particular length, which leads to performance degradation for large models. In our approach, we employ the symbolic model checker SpaceEx, which prunes the state space exploration by checking whether the currently considered states have already been visited.

Conclusion

We have presented a novel parameter identification algorithm for multiaffine hybrid automata. In our algorithm, we compute equivalence classes in the parameter space and explore them in a hierarchical way. The approximation of the system dynamics with linear hybrid automata lets us keep the timing information in our abstraction. This allows us to precisely treat time-dependent properties such as a stimulus.

Given a parameter polytope P , we compute an LHA which overapproximates the system behavior for P . Furthermore, we compute another LHA which enables us to prune the search tree. We have evaluated our approach on a model of a genetic regulatory network and a myocyte model and demonstrated its improvement over RoVerGeNe, a tool for parameter identification based on a purely discrete abstraction.

In the future, we plan to investigate the application of hybrid model checkers which support more expressive continuous dynamics. This enables approximating the parametrized system dynamics with a hybrid automaton class featuring dynamics beyond the ones of LHA.

Fig. 1 .

 1 Fig. 1. (a) Ramp functions r + and r -. θi and θ i are threshold parameters. (b) Gene network comprising two genes, a and b, coding for two repressor proteins, A and B. Each protein represses the expression of the other gene, forming a cross-inhibition network. (c) PMA model of the network in (b). Because of its simplicity, this model is actually piecewise-affine. (d) Known and uncertain parameter values. (e) Bistability property expressed in LTL.

Fig. 1 .

 1 Fig. 1. (a) Ramp functions r + and r -. (b) The two-genes network model.

Fig. 2 .

 2 Fig. 2. (a) From state space partition to transition system. Example for κa = 10 and κ b = 15 (arrows normalized). (b-c) Flow computation for P = [8, 12] × [15, 20].

Fig. 3 .

 3 Fig. 3. Relations of the systems presented in this paper. Let p ∈ P . An arrow S1 -→ S2 indicates that S2 is an overapproximation of S1.

Definition 7 .

 7 Given an LHA H with set of locations Loc, let S be a set of discrete states with |S| = |Loc| and let disc : Loc → S be a bijection which maps every location to a discrete state. In addition, let Bad be the bad states of H. The pair (H, Bad ) induces a Kripke structure K = (S, S 0 , S B , T ), where S = {disc( ) | ∈ Loc} is the set of states, S 0 = {disc( ) | inv ( ) ∩ Init = ∅} is the set of initial states, S B = {disc( ) | inv ( ) ∩ Bad = ∅} are the bad states, and T = {disc( ) → K disc( ) | → H ∈ Trans} is the set of transitions.

17 %Fig. 4 .

 174 Fig. 4. The algorithm in pseudocode.

Fig. 5 .

 5 Fig.5. Search tree in the parameter space. yellow: L ∃ M (P ) satisfies the property ϕ while K ∃ M (P ) violates it. blue: L ∀ M (P ) satisfies the property ϕ while K ∀ M (P ) violates it. gray: The node is only explored by the LHA analysis. red: The node is only explored by the KS analysis.

  valid # sets # nodes # ∃-KS/∃-LHA # ∀-KS/∀-LHA runtime [s]

	model ID	RoVerGeNe	Hydentify	RoVerGeNe	Hydentify	RoVerGeNe	Hydentify	RoVerGeNe	Hydentify KS LHA	RoVerGeNe	Hydentify KS LHA	RoVerGeNe	Hydentify
	1 60 85	7 5	23 23	23	9	21	16	5	10	11	32
	2	0 38	0 5	1	79	1	1	79	1	1	62	2	95
	3 65 73	3 5	9	15	9	9	12	5	5	3	6	22
	4 60 84	4 3	13	9	13	7	7	9	4	1	8	18
	Table 1. model ID: 1-2: without stimulus; 3-4: with stimulus; % valid: percentage of
	parameter space verified; # sets: number of parameter sets found; # nodes: number
	of nodes in the search tree; # ∃-KS/∃-LHA: number of K ∃ M /L ∃ M analyzed; # ∀-
	KS/∀-LHA: number of K ∀ M /L ∀ M analyzed; runtime: runtime in seconds

1

. %

b := ⊥ ; % no f u t u r e v a l i d s e t s f o r t h e KS a n a l y s i s

http://swt.informatik.uni-freiburg.de/tool/spaceex/hydentify

Acknowledgments. This work was partly supported by the European Research Council (ERC) under grant 267989 (QUAREM), by the Austrian Science Fund (FWF) under grants S11402-N23, S11405-N23 and S11412-N23 (RiSE/SHiNE) and Z211-N23 (Wittgenstein Award), and by the German Research Foundation (DFG) as part of the Transregional Collaborative Research Center "Automatic Verification and Analysis of Complex Systems" (SFB/TR 14 AVACS, http://www.avacs.org/).

V := V ∪ P ; return ; % v a l i d s e t found by t h e KS a n a l y s i s


Federico Cacciafesta

Eric Séré

Éric S Ér

LOCAL SMOOTHING ESTIMATES FOR THE MASSLESS DIRAC-COULOMB EQUATION IN 2 AND 3 DIMENSIONS

whether they are published or not. The documents may come

Introduction and generalities

The massless Dirac equation with an electric Coulomb potential reads (1.1)

   iu t + D n u - ν |x| u = 0, u(t, x) : R t × R n x → C N u(0, x) = f (x)
where the massless Dirac operator D n is defined in terms of the Pauli matrices (1.2)

σ 1 = 0 1 1 0 , σ 2 = 0 -i i 0 , σ 3 = 1 0 0 -1 as D 2 = -i(σ 1 ∂ x + σ 2 ∂ y ) = 0
-i∂ z -i∂ z 0 (we denote ∂ z = ∂ x -i∂ y and ∂ z = ∂ x + i∂ y ) in dimension n = 2 with N = 2, and

D 3 = -i 3 k=1 α k ∂ k = -i(α • ∇)
where the 4 × 4 Dirac matrices are given by (1.3)

α k = 0 σ k σ k 0 , k = 1, 2, 3 in dimension n = 3, with N = 4.
The Dirac equation is widely used in physics to describe relativistic particles of spin 1/2 (see e.g. [START_REF] Landau | Quantum mechanics -Relativistic quantum theory[END_REF], [START_REF] Thaller | The Dirac Equation[END_REF]). In particular, the 2D massless equation is used as a model for charge carriers in graphene, a layer of carbon atoms arranged in a honeycomb lattice (see e.g. the physics survey [START_REF] Neto | The electronic properties of graphene[END_REF] and the mathematical papers [START_REF] Fefferman | Wave packets in honeycomb structures and twodimensional Dirac equations[END_REF], [START_REF] Hainzl | Ground state properties of graphene in Hartree-Fock theory[END_REF], [START_REF] Maier | Stability of impurities with Coulomb potential in graphene with homogeneous magnetic field[END_REF]). Note that there have been many rigorous results in the recent years on the stationary (massive) Dirac equation, essentially in dimension 3 (see e.g. [START_REF] Esteban | Variational methods in relativistic quantum mechanics[END_REF] for references). Comparatively, its dynamics has been less investigated up to now.

The σ j and α j matrices were introduced in view of making the Dirac operator a square root of the Laplace operator: therefore, they satisfy by construction the following anticommutating relations σ j σ k + σ k σ j = 2δ ik I 2 , j, k = 1, 2; α j α k + α k α j = 2δ ik I 4 , j, k = 1, 2, 3, (we recall that the spectrum of D n is the whole line R). These conditions ensure that (i∂ t -D n )(i∂ t + D n ) = (∆ -∂ 2 tt )I N , which strictly links the dynamics of the free massless Dirac equation to a system of N decoupled wave equations. Therefore, dynamical properties for the free Dirac equation can directly be derived from their wave counterparts: dispersive properties of the flow can be mainly encoded in the celebrated family of Strichartz estimates, which are given by (1.4) e itDn f

L p t Ḣ 1 q -1 p -1 2 x f L 2
where the exponents (p, q) are wave admissible, i.e. satisfy

2 p + d -1 q = d -1 2 , 2 ≤ p ≤ ∞, 2(d -1) d -3 ≥ q ≥ 2
and e itDn represents the propagator for the solutions to (1.1) with V = 0 for n = 2, 3. We stress the fact that the estimate corresponding to the couple (p, q) = (2, ∞) in dimension 3, the so called endpoint, that is

(1.5) e itD 3 f L 2 t L ∞ x f Ḣ1
is known to fail as the corresponding one for the 3D wave flow. In [START_REF] Machihara | Endpoint Strichartz estimates and global solutions for the nonlinear Dirac equation[END_REF] a refined version of this estimate involving angular spaces was given. As a consequence, the well-posedness of the celebrated cubic nonlinear Dirac equation was established, assuming initial data in Ḣ1 with slight additional regularity on the angular variable. The problem was afterwards completely solved in [START_REF] Bejenaru | The cubic Dirac equation: Small initial data in H 1 (R 3 )[END_REF]. See also [START_REF] Bejenaru | The cubic Dirac equation: Small initial data in H 1/2 (R 2 )[END_REF] for the 2D case.

Another family of a priori estimates encoding informations about dispersion is given by the so called local-smoothing type estimates, which turn to be particularly useful to handle potential-type perturbations. Different versions of these estimates are available: to the best of our knowledge the most general one in this setting is (1.4) and (1.6) can be found in [START_REF] Fanelli | Strichartz and smoothing estimates of dispersive equations with magnetic potentials[END_REF] (see also [START_REF] Escobedo | A semilinear Dirac equation in H s (R 3 ) for s > 1[END_REF], [START_REF] Boussaid | Stable directions for small nonlinear Dirac standing waves[END_REF]). We are not aware of any corresponding result for the 2D case; nevertheless, the strategy developed in [START_REF] Fanelli | Strichartz and smoothing estimates of dispersive equations with magnetic potentials[END_REF] can be suitably adapted to this context to obtain similar results.

(1.6) w -1/2 σ e itD 3 f L 2 t L 2 x ≤ C f L 2 where w σ (x) = |x|(1 + | log |x||) σ , σ > 1. Both
In recent years a lot of effort has been spent in order to investigate what happens to dispersive flows when potential perturbations come into play. In [START_REF] Fanelli | Strichartz and smoothing estimates of dispersive equations with magnetic potentials[END_REF] both (1.4) and (1.6) were proved for the flow e it(D 3 +V ) under the assumption

|V (x)| ≤ δ w σ (x) , σ > 1, n = 3,
for some δ sufficiently small, which by the way is a sufficient condition to guarantee that the operator D 3 + V is selfadjoint (see [START_REF] Thaller | The Dirac Equation[END_REF]). Similar results are available for magnetic potentials as well (see e.g. [START_REF] Boussaid | Virial identity and weak dispersion for the magnetic Dirac equation[END_REF], [START_REF] Cacciafesta | Virial identity and dispersive estimates for the n-dimensional Dirac equation[END_REF]). The angular endpoint estimates proved in [START_REF] Machihara | Endpoint Strichartz estimates and global solutions for the nonlinear Dirac equation[END_REF] have been generalized in [START_REF] Cacciafesta | Global small solutions to the critical Dirac equation with potential[END_REF], [START_REF] Cacciafesta | Endpoint estimates and global existence for the nonlinear Dirac equation with a potential[END_REF] to small potential perturbations, by introducing some new mixed Strichartz-smoothing estimates. At the same time, some effort has been spent in order to find examples of potentials such that the correspoding flows do not disperse (in the sense discussed above). This problem has been tackled for the magnetic Dirac equation in [START_REF] Arrizabalaga | On the lack of dispersion for a class of magnetic Dirac flows[END_REF] in 3D, essentially showing that for vector potentials of the form A(x) ∼ = |x| -δ Mx with δ ∈ (1, 2) most of the mass of the solution is localized around a non-dispersive function. These arguments suggest that, heuristically, the degree of homogeneity of the operator works as a threshold for the validity of dispersive estimates: therefore, for the Dirac equation, the Coulomb potential is to be thought of as a critical case (we stress the fact that the restriction to the massless case seems unavoidable now). In the very last years the problem of understanding dispersive estimates for scaling invariant potentials has been approached for other dynamics: in [START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with the inverse-square potential[END_REF] and [START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with potentials of critical decay[END_REF] the authors have proved indeed that Morawetz and Strichartz estimates can be recovered for both the Schrödinger and wave equations with inverse square potentials, and in [START_REF] Fanelli | Time decay of scaling critical electromagnetic Schrödinger flows[END_REF] the L 1 → L ∞ time decay estimate is proved in the electromagnetic case. The aim of this paper is to adapt the techniques of [START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with the inverse-square potential[END_REF] to the massless Dirac Coulomb equation, and thus to prove local smoothing estimates for solutions to (1.1). The strategy of proof of [START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with the inverse-square potential[END_REF] is fairly straightforward: the basic idea is to use a spherical harmonics decomposition to reduce the equation to a radial one and then rely on the Hankel transform to diagonalize it. Several problems arise when trying to adapt these ideas to the Dirac setting. The first difficulty is given by the fact that the Dirac operator does not preserve radiality, and therefore one needs a more complex, but still classical, spherical harmonics decomposition related to the SU 2 group action (see [START_REF] Thaller | The Dirac Equation[END_REF], section 4.6). Another issue is that we cannot use the Hankel transform, and we have to build up a relativistic analogue. To do this, we need to study the continuous spectrum of the Dirac-Coulomb model, both in 2 and 3 dimensions. The corresponding generalized eigenfunctions are well known (see for instance [START_REF] Landau | Quantum mechanics -Relativistic quantum theory[END_REF] and [START_REF] Novikov | Elastic scattering theory and transport in graphene[END_REF]), but their rich structure will raise some additional technical difficulties. We stress the fact that, to the best of our knowledge, this is the first result concerning dispersive properties of Dirac-Coulomb, an operator which plays a fundamental role in relativistic quantum chemistry. We hope that our result will be useful in the study of the dynamics of nonlinear models involving, for instance, several relativistic electrons in a molecule.

Before stating our main result, we introduce some notations that will be used throughout the paper.

Notations.

We shall use the standard notation Ḣs for the homogeneous Sobolev space with the norm f Ḣs = |D| s f L 2 where |D| = (-∆) 1/2 , and L p t L q x = L p (R t ; L q (R n x )) for the mixed space-time Strichartz spaces. We denote with Ω s the operator

(Ω s φ)(x) = |x| s φ(x)
and with a little abuse of notation we will use the same symbol to indicate the operators which are pointwise equal for all times, (Ω s ψ)(t, x) = |x| s ψ(t, x).

We introduce the following sets of indices

(1.7) B 2 (k) := {k : k = 1/2 + m, m ∈ Z} (1.8) B 3 (k) := {k : k ∈ Z\{0}}.
As most of the forthcoming objects will be significantly different in 2D or 3D, we will often use the apex n to distinguish the dimensions; this will not have to be confused with standard powers, but interpretation will be clear from time to time. For the definition of the angular spaces H n k , n = 2, 3 we refer to forthcoming Proposition 2.1 and subsequent Remark 2.1.

We are now ready to state our main Theorem.

Theorem 1.1. Let n = 2, 3, ν n ∈ -n-1 2 , n-1 2
, and let k n ∈ B n (k). Let u be a solution of (1.1). Then, for any

1/2 < α < k 2 n -ν 2 n + 1/2 and any f ∈ L 2 ((0, ∞)) ⊗ H n
≥kn there exists a constant C = C(ν, α, k) such that the following estimate holds

(1.9) Ω -α D n - ν n |x| 1/2-α u L 2 t L 2 x ≤ C f L 2 .
Remark 1.1. The dependence on the dimension of ν, which will be neglected in the rest of the paper, has been added in the statement of the Theorem to allow us to remark the different ranges of admissible charges in dimensions 2 and 3. These ranges are in fact imposed by the condition on the Dirac Coulomb operator to be self-adjoint, which is indeed guaranteed in those intervals (see Remark 1.4). We stress also the fact that for ν approaching the endpoints, respectively ±1 in 3D and ±1/2 in 2D, the range of admissible exponents for the lowest value of k n in estimate (1.9) shrinks to zero.

Remark 1.2. It seems to be possible to generalize our proof of Theorem 1.1 to any space dimension n ≥ 2. As it will be clear, the key ingredients are given indeed by the construction of a suitable integral operator which transforms the problem into an ODE, and then by the careful analysis of some interaction integrals between generalized eigenstates of the Dirac-Coulomb equation that naturally appear. Therefore, the crucial step consists in finding the generalized eigenstates of the Dirac-Coulomb equation: this has been performed in [START_REF] Gu | Exact solutions to the Dirac equation for a Coulomb potential in D + 1 dimensions[END_REF] in higher space dimension n ≥ 4 by generalizing the spherical wave decomposition (see Proposition 2.1), i.e. by carefully analyzing the symmetries of the SO(n) group. Anyway, we prefer not to deal with this problem here both because technicalities would lead us too far and because of the scarce physical interest of it.

Remark 1.3. The endpoint case α = 1/2 in (1.9) would permit to recover with a standard argument the full set of Strichartz estimates for the perturbed flow (see e.g. [START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with the inverse-square potential[END_REF]); even though we don't have a concrete counterexample, this estimate seems to fail. Forthcoming remark 3.1 gives more details on the issue.

Remark 1.4. In order to have a unitary flow for equation (1.1) one needs the operator D n -ν |x| to be selfadjoint; many papers have been devoted to the study of this property (see e.g. [START_REF] Schmincke | Distinguished self-adjoint extensions of Dirac operators[END_REF], [START_REF] Wüst | A convergence theorem for self-adjoint operators applicable to Dirac operators with cut-off potentials[END_REF], [START_REF] Nenciu | Self-adjointness and invariance of the essential spectrum for Dirac operators defined as quadratic forms[END_REF], [START_REF] Klaus | Characterization and uniqueness of distinguished self-adjoint extensions of Dirac operators[END_REF], [START_REF] Esteban | Self-adjointness for Dirac operators via Hardy-Dirac inequalities[END_REF], [START_REF] Esteban | Self-adjointness via partial Hardy inequalities Mathematical results in quantum mechanics[END_REF], [START_REF] Arrizabalaga | Self-adjoint extensions of Dirac operators with Coulomb type singularity[END_REF]) in the 3D setting (note that the methods introduced in these papers are quite general and can be adapted to other dimensions). In 3D the essential self-adjointness is in fact guaranteed if |ν| ≤ √ 3 2 , while in the interval √ 3

2 < |ν| ≤ 1 it is still possible to build a distinguished selfadjoint extension. In the 2D case, as soon as ν is nonzero, the Dirac-Coulomb operator D 2 -ν |x| defined on C ∞ 0 (R 2 ) is not essentially selfadjoint (this is an immediate consequence of Theorem 4.16 in [START_REF] Thaller | The Dirac Equation[END_REF]). But C. Warmt ([41], Satz 2.2.6) recently proved the existence of a distinguished self-adjoint extension if |ν| < 1/2. Therefore, in the range of ν n for which we prove our Theorem 1.1, the operator D n -ν |x| is self-adjoint (in the sense just discussed). The plan of the paper is the following: in section 2 we build the setup, reviewing the theory of partial wave subspaces and the spectrum of the Dirac-Coulomb operator in 2D and 3D; then we define our analogue of the Hankel transform used in [START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with the inverse-square potential[END_REF]. Section 3 is devoted to the proof of Theorem 1.1. A short appendix containing some generalities on special functions is included at the end for the reader's convenience.

(Italy). The second named author acknowledges support from the ANR "NoNAP" of the french ministry of research.

2. The setup.

We devote this section to recall some classical facts about the radial Dirac operator and his spectrum, including a short discussion on how to explicitly calculate the generalized eigenstates of the continuous spectrum. For the sake of brevity we will omit several details, referring to [START_REF] Landau | Quantum mechanics -Relativistic quantum theory[END_REF] and [START_REF] Thaller | The Dirac Equation[END_REF] for a complete discussion and for more informations on the topic. Also, we construct an integral transform that will be crucial in our proof.

2.1. Partial wave decomposition and continuous spectrum. The main ingredient we need to introduce is the so called partial wave decomposition. We collect in Proposition 2.1 the tools we will need in what follows that is, essentially, the fact that the Dirac-Coulomb operator can be seen as a radial operator with respect to some suitable decomposition, both in dimensions 2 and 3. We will denote for brevity with A 2 and A 3 the following families of indices (2.1)

A 2 := k ∈ Z + 1 2 ;
(2.2)

A 3 := j = 1 2 , 3 2 , . . . , m j = -j, -j + 1, . . . , +j, k j = ±(j + 1/2)
Moreover, we define the 4-vectors

(2.3) Ξ + m j ,∓(j+1/2) = iΩ m j j∓1/2 0 , Ξ - m j ,∓(j+1/2) = 0 Ω m j j±1/2 , where Ω m j j-1/2 = 1 √ 2j j + m j Y m j -1/2 j-1/2 j -m j Y m j +1/2 j-1/2 Ω m j j+1/2 = 1 √ 2j + 2 j + 1 -m j Y m j -1/2 j+1/2 -j + 1 -m j Y m j +1/2 j+1/2
and Y m l are the usual spherical harmonics. Proposition 2.1. Let n = 2, 3. Then, we can define unitary isomorphisms between the Hilbert spaces

L 2 (R 2 ) 2 ∼ = L 2 ((0, ∞), dr) k∈A 2 h 2 k , L 2 (R 3 ) 4 ∼ = L 2 ((0, ∞), dr) j,m j ,k j ∈A 3 h 3 j,m j ,k j
by the following decompositions

(2.4) Φ 2 (x) = k∈A 2 1 2 √ π F k (r)e i(k-1/2)θ G k (r)e i(k+1/2)θ , (2.5) Φ 3 (x) = j,m j ,k j ∈A 3 1 2 √ π F j,m j ,k j (r)Ξ + j,m j ,k j (θ, φ) G j,m j ,k j (r)Ξ - j,m j ,k j (θ, φ)
which hold for any

Φ n ∈ L 2 (R n , C N ).
Moreover, the Dirac-Coulomb operator leaves invariant the partial wave subspaces C ∞ 0 ((0, ∞)) ⊗ h n k and, with respect to the basis e i(k-1/2)θ , e i(k+1/2)θ and Ξ + j,m j ,k j , Ξ - j,m j ,k j is respectively represented by the radial matrices

(2.6) d 2 k = -ν r -d dr + k r d dr + k r -ν r , d 3 j,m j ,k j = -ν r -d dr + 1+k j r d dr - 1-k j r -ν r . The Dirac-Coulomb opeator D n -ν r on C ∞ 0 (R n ) N is unitary equivalent to the direct sum of "partial wave Dirac operators" d j,m j ,k j , D 2 - ν r ∼ = k d 2 k , D 3 - ν r ∼ = j,m k ,k k d 3 j,m k ,k k For fixed k 2 ∈ B(k 2 ) and k 3 ∈ B(k 3 ) we denote with (2.7) H 2 ≥k 2 = k:|k|≥|k 2 | h k , H 3 ≥k 3 = j:j≤|k 3 |-1/2,m j ,k j : h j,m j ,k j .
Proof. See [START_REF] Thaller | The Dirac Equation[END_REF].

Remark 2.1. Notations in Proposition 2.1 are very heavy; this is partly due to the richness of the Dirac operator itself and partly to the fact that the spherical coordinates decomposition changes in dimensions 2 and 3. We should stress the fact that formulas (2.6) only depend on k: therefore, in what follows we will abuse notations neglecting the dependence on the other parameters in the case n = 3 in order to have a unified presentation.

We now show how the functions of the continuous spectrum for the Dirac-Coulomb operator are obtained. We want to deal with the eigenvalue equation

(2.8) D n - ν |x| Ψ n = εΨ n , ε > 0, n = 2, 3.
The application of Proposition 2.1 suggests to seek for solutions in the form, for a fixed value of ε,

(2.9) Ψ 2 k,ε (r, x) = F 2 k,ε (r)e i(k-1/2)θ G 2 k,ε (r)e i(k+1/2)θ , n = 2 (2.10) Ψ 3 k,ε (r, x) = F 3 k,ε (r)Ξ + k (x) G 3 k,ε (r)Ξ - k (x)
, n = 3 so that equation (2.8) will in fact reduce to a sytem of ODE for the radial components F n k (r), G n k (r) which is the same in both dimensions, apart for some constants, and can be expliticly solved (for the details on these calculations we refer to [START_REF] Landau | Quantum mechanics -Relativistic quantum theory[END_REF] section IV for the 3D and to [START_REF] Novikov | Elastic scattering theory and transport in graphene[END_REF] for the 2D case). We thus have the following formulas, for fixed values of k and ε > 0, and n = 2, 3, (2.11)

F n k,ε = √ 2|Γ(γ + 1 + iν)| Γ(2γ + 1) e πν/2 (2εr) γ-(n-1)
2

Re e i(εr+ξ)

1 F 1 (γ -iν, 2γ + 1, -2iεr) (2.12) G n k,ε = i √ 2|Γ(γ + 1 + iν)| Γ(2γ + 1) e πν/2 (2εr) γ-(n-1)
2

Im e i(εr+ξ) 1 F 1 (γ -iν, 2γ + 1, -2iεr)
where 1 F 1 (a, b, z) are confluent hypergeometric functions (see Appendix for the definition and some properties), γ = √ k 2 -ν 2 and e -2iξ = γ-iν k is a phase shift. We will denote with (2.13)

ψ n k,ε (r) = F n k,ε (r) G n k,ε (r)
, n = 2, 3 the vector of radial coordinates of the generalized eigenfunctions.

Notice that with this choice the following normalization condition holds for (2.9)-(2.10)

(2.14) +∞ 0 ψ n ǫ (r)ψ n ǫ ′ (r)r n-1 dr = δ(ǫ -ǫ ′ ) ǫ n-1 , n = 2, 3
meaning with that that the function

K(r, r ′ ) = r n-1 e -lε ψ n ǫ (r)ψ n ǫ ′ (r)ε n-1 dε is a δ-approximation for l → 0.
Remark 2.2. We point out that, due to our choice of normalization (2.14), the homogeneity of the functions ψ n (r) with respect to ε and r is the same; this fact, that is not true anymore in the massive case, will be crucial in what follows, and therefore prevents the application of our strategy to the massive case. Remark 2.4. Formulas (2.11)-(2.12) give the "eigenstates" of the continuous spectrum corresponding to positive energies ε > 0. Using a charge conjugation argument (or directly working on the explicit radial equations), it is possible to obtain the corresponding formulas for the radial components (2.9) for negative energies: it turns indeed out that the function ψ -k,ε (r) solves

(2.15) dk ψ -k,ε (r) = -εψ -k,ε , n = 2, 3,
where the operator dk is the one defined in (2.6) with opposite sign of the charge ν. This shows that the "eigenstates" corresponding to negative energies for equation (2.8) can be understood by mean of (2.15) to be "eigenstates" corresponding to positive energies, with opposite value of j and opposite charge. Therefore, denoting with F , G the functions obtained by (2.11)-(2.12) by changing the sign of ν, we obtain for the radial coordinates of the "eigenstates" corresponding to negative energies as in (2.13) the following relation

ψ n k,-ε (r) = F n k,-ǫ (r) G n k,-ǫ (r) = F n -k,ǫ (r) Gn -k,ǫ (r)
.

2.2. The integral transform. We now introduce the crucial integral transform that will be used in the proof of the main result, consisting in a projection on the continuous spectrum of the Dirac-Coulomb operator. Throughout this subsection we will omit the dependance on the dimensions of the functions, as the calculations carried out will be the same.

Definition 2.2. Let Φ ∈ L 2 ((0, ∞), dr) ⊗ h k for some fixed k and let ϕ(r) = (ϕ 1 (r), ϕ 2 (r)) be the vector of its radial coordinates in decomposition (2.4)-(2.5). We define the following integral transform (2.16)

P k ϕ(ǫ) = P + k ϕ(ǫ) P - k ϕ k (ǫ) = +∞ 0 ψ k,ε (r)ϕ(r)r n-1 dr C +∞ 0 ψ k,-ε (r)ϕ(r)r n-1 dr = +∞ 0 H k (εr) • ϕ(r)r n-1 dr
where we have introduced the matrix (keep in mind Remark 2.4)

(2.17)

H k = F k,ε (r) G k,ε (r) F k,-ε (r) G k,-ε (r)
We collect in the following proposition some important properties of the operator P k . Proposition 2.3. For n = 2, 3 and any Φ ∈ L 2 ((0, ∞), dr) ⊗ h k the following properties hold:

(1) P k is an L 2 -isometry.

(

) P k d k = σ 3 ΩP k . 2
(3) The inverse transform of P k is given by

(2.18) P -1 k ϕ(r) = +∞ 0 H * k (ǫr) • ϕ(ǫ)ǫ n-1 dǫ where H * k = F k,ε (r) G k,-ε (r) G k,ε (r) F k,-ε (r) .
(4) For every σ ∈ R we can define the fractional operators

(2.19) A σ k ϕ k (r) = P k σ 3 Ω σ P -1 k ϕ k (r) = +∞ 0 S σ k (r, s) • ϕ k (s)s n-1 ds.
where the integral kernel S k (r, s) is the 2 × 2 matrix given by Proof. Property ( 1) and ( 2) come from the definition of P k , once noticed that

(2.20) S k (r, s) = +∞ 0 H k (ǫr) • H * k (ǫs)ǫ n-
P ± k (d k ϕ) = ψ k,±ǫ , d k ϕ = d k ψ k,±ǫ , ϕ = ±ǫ ψ k,±ǫ , ϕ and from normalization relation (2.14). Property (3) is a consequence of (1) (H * k is the adjoint of H k ).
To prove property (4) we use the definition of A k to write

A σ k ϕ k (r) = P k Ω σ P -1 k ϕ k (r) = +∞ 0 H k (ǫr)ǫ n-1+σ +∞ 0 H * k (ǫs)ϕ k (s)s n-1 ds dǫ.
Exchanging the order of the integrals yields (4).

Proof of Theorem 1.1

We here prove our main result. Following [START_REF] Burq | Strichartz estimates for the wave and Schrödinger equations with the inverse-square potential[END_REF], we rely on partial wave decomposition and on Proposition 2.3 to explicitly write the solutions on the single subspaces; then, a detailed and careful analysis of the kernel S k (r, s) of the operator A k defined in (2.19) will be performed to sum back.

3.1. The strategy. We start with decomposition (2.4)-(2.5) and work on a single spherical space which will be identified by k ∈ A 2 as in (2.1) and by the triple {j, m k , k j } ∈ A 3 as in (2.2) in the case n = 3. With a slight abuse of notations, we will again neglect the dependence on everything but k in the 3D case (as noticed in Remark 2.1, the radial Dirac-Coulomb operator in dimension 3 only depends on k = k j , see (2.6)), as well as the dependence on the dimension for the various objects, in order to perform much readable and unified calculations. We thus rely on Proposition 2.1 aiming to prove estimate (1.9) on a fixed partial wave subspace with a constant C bounded with respect to k, and eventually sum back. We thus take an initial condition f with angular part in h k for some fixed k, and denote with L k f the solution to the initial value problem

(3.1) iu t + d k u = 0, u(0, x) = f (x)
with d k as in (2.6). We apply operator P k , which is an L 2 -isometry, to the LHS side of estimate (1.9) and use Proposition 2.3 to obtain (the application of the matrix σ 3 does not alter the L 2 norm)

P k Ω -α |d k | 1/2-α L k f L 2 t L 2 x = A -α k Ω 1/2-α P k L k f L 2 t L 2 x
where we have used (2.19). Now,

P k L k f solves (see (2) of Proposition 2.1) (3.2) i∂ t P k L k f + σ 3 ΩP k L k f = 0, P k L k f (0, ξ) = P k f (ξ),
so that the solution to this problem is explicitly given by

P k L k f (t, ξ) = e itσ 3 ξ P k f (ξ).
Fourier transforming in time (since F t→τ is an L 2 -isometry which commutes with A σ k ) gives (with a minor abuse, we will neglect the σ 3 matrix from now on as we will only need to pick the L 2 norm and thus the sign will be unessential)

(F t P k L k f )(τ, ξ) = (P k f )(ρ)δ(τ + ξ).
We can then write

(A -α k Ω 1/2-α F t P k L k f )(τ, ξ) = +∞ 0 S -α k (ξ, s)δ(τ + s)P k f (s)s 2n-2α-1 2 ds = -S -α k (ξ, τ )P k f (τ )τ 2n-2α-1 2
.

In view of proving (1.9) we now need to calculate the L 2 norm in time and space of the quantity above which then gives (notice that, as the angular part in decomposition (2.4)-(2.5) is L 2 -unitary, we only need to consider the radial integrals)

(3.3) +∞ 0 +∞ 0 ((P k f ) * (τ )S -α k (ρ, τ ) T ) • (S -α k (ρ, τ )(P k f )(τ ))τ 2n-2α-1 ρ 2 dρ.
Since S p k (ρ, τ ) T = S p k (τ, ρ), the integral in dρ yields S -2α k (τ, τ ) and we are therefore left with

(3.4) +∞ 0 (P k f ) * (τ )S -2α k (τ, τ )(P k f )(τ )τ 2n-2α-1 dτ.
To conclude the proof of the estimate on the k-th space we need to bound (3.4) with the L 2 -norm of P k f (which we recall to be the L 2 norm of f ). Using forthcoming Proposition 3.1 we thus have

(3.5) +∞ 0 (P k f ) * (τ )S -2α k (τ, τ )(P k f )(τ )τ 2n-2α-1 dτ ≤ +∞ 0 Tr(S -2α k (τ, τ ))|(P k f )(τ )| 2 τ 2n-2α-1 dτ ≤ C k +∞ 0 |(P k f )(τ )| 2 (τ )τ n-1 dτ = C k f L 2 .
with C k a bounded constant of k. Expanding f in partial waves and using the triangle inequality together with the orthogonality of spherical harmonics concludes the proof.

The main integral. The crucial step consists thus in evaluating

Tr(S -2α (τ, τ )) = Re +∞ 0 F k,ε (τ ) 2 + G k,ε (τ ) 2 ǫ n-1-2α dǫ + Re +∞ 0 F k,-ε (τ ) 2 + G k,-ε (τ ) 2 ǫ n-1-2α dǫ (3.6) = Re I 1 (τ ) + I 2 (τ ) .
More than giving an explicit solution of such integrals, which seems to be possible but significantly complicated, we are more interested in proving the following Proposition 3.1. Let n = 2, 3 k n ∈ B(k) as defined in (1.7)-(1.8) and ν n as in Theorem 1.1. Then for every

1/2 < α < k 2 n -ν 2 n + 1/2 there exists a constant C kn such that Tr(S -2α k (τ, τ )) = C kn τ -n+2α . Moreover, sup kn |C kn | < ∞ Proof.
By recalling the structure of the generalized eigenstates (2.11)-(2.12), we have that integral (3.6) takes the form

(3.7) I 1 (τ ) = C k,ν τ 2γ-(n-1) × +∞ 0 ǫ 2γ-2α e 2iǫτ 1 F 1 (γ -iν, 2γ + 1, -2iǫτ ) 1 F 1 (γ -iν, 2γ + 1, -2iǫτ )dǫ with C k,ν = 2 2γ+2 π e πν |Γ(γ + 1 + iν| 2 Γ(2γ + 1) 2 e 2iξ where we recall that γ = √ k 2 -ν 2 .
As the only difference between I 1 and I 2 is encoded in the sign of ν (see Remark 2.4) which is virtually unessential, we will limit our discussion to I 1 , the other case being completely analogous.

Integrals like (3.7) have already been object of study in literature (see e.g. [START_REF] Landau | Quantum mechanics -Nonrelativistic theory[END_REF], [START_REF] Gargaro | Real and Virtual Radiation in Electron-Nucleus Scattering[END_REF], [START_REF] Biedenharn | Some Exact Radial Integrals for Dirac-Coulomb Functions[END_REF]) since they are connected with some different physical problems, as for instance matrix multipole elements for the Dirac Coulomb model, and they present some difficulties being not convergent at infinity as well as in the origin (see (4.3)). Singularity at the origin can be dealt with by a standard integration by parts argument (as the one used to define Γ function for negative real parts) and is discussed in [START_REF] Sun | Radial integrals with finite energy loss for Dirac Coulomb functions[END_REF] provided 2γ -2α + 1 is not zero or a negative integer, while singularity at infinity is slightly more difficult to be handled. To make things work, we need first of all to define I k (r, s) off-diagonal and then use a limiting procedure. Also, the introduction of a term of the form e -δǫ with δ ∈ R + in the integral is needed in order to make it convergent; we will then denote with I δ k the integral

I δ k (r, s) = C k,ν τ 2γ-(n-1) × +∞ 0 ǫ 2γ-2α e ǫ(i(r+s)-δ) 1 F 1 (γ -iν, 2γ + 1, -2iǫr) 1 F 1 (γ -iν, 2γ + 1, -2iǫs)dǫ.
This modified integral can be explicitly computed, and its solution is written in terms of a double Appel series, which unfortunately turns out to be not convergent for the values we are interested in. To deal with this problem, it is necessary to introduce an analytic continuation of such a series: we don't include the proofs of these results here as they can be found in corresponding references. We begin with the following where F 2 is defined in (4.5).

Proof. See e.g. Lemma 1 in [START_REF] Saad | Integrals containing confluent hypergeometric functions with applications to perturbed singular potentials[END_REF].

We should now point out that the ray of convergence of F (α; a 1 , a 2 ; b 1 ; b 2 ; x, y) is |x| + |y| < 1, which provides a first obstacle to the application of this result to our case. To overcome this difficulty, we rely on the following analytic continuation of

F 2 . Lemma 3.3. If Re(b 1 -a 1 ) > 0, Re(a 1 ) > 0 (3.9) F 2 (d; a 1 , a 2 , b 1 , b 2 ; 1, 1) = e -iπ(b-a-d) Γ(b 1 )Γ(d -b 1 + 1) Γ(a 1 )Γ(d -a 1 + 1) 3 F 2 (d, b 2 -a 2 , d -b 1 + 1; b 2 , d -a 1 + 1; 1) +e iπa Γ(b 1 )Γ(d -b 1 + 1) Γ(b 1 -a 1 )Γ(d + a 1 + 1 -b 1 ) 3 F 2 (d, a 2 , d + 1 -b 1 ; b 2 , d + a + 1 -b 1 ; 1)
Proof. See [START_REF] Biedenharn | Some Exact Radial Integrals for Dirac-Coulomb Functions[END_REF].

For the sake of generality we have stated this Lemma with general values of the parameters; we list here for convenience the necessary substitutions to treat our case:

d = 2γ + 1 -2α; (3.10) a 1 = a 2 = γ -iν; b 1 = b 2 = 2γ + 1.
Since necessary conditions of Lemma 3.3 are fullfilled with these values, this result allows us to pass to the limit I δ k (r, s) for δ → 0, s → r =: τ , and thus leaves us with

I 1 (τ ) = C k,ν τ 2γ-(n-1) (2iτ ) -(2γ+1-2α) Γ(2γ + 1 -2α) ×F 2 (2γ + 1 -2α; γ -iν, γ -iν; 2γ + 1, 2γ + 1; 1, 1) with C k,ν = 2 2γ+2 π e πν |Γ(γ + 1 + iν| 2 Γ(2γ + 1
) 2 e 2iξ . Using now representation (3.9) for F 2 and plugging in values (3.10), we are led to (3.11)

I 1 (τ ) = τ 2α-n C k,ν (2i) -(2γ+1-2α) Γ(2γ + 1 -2α)(A + B)
where

A = e iπ(γ-2α-iν) Γ(2γ + 1)Γ(1 -2α) Γ(γ -iν)Γ(γ + 2 -2α + iν) × m (2γ + 1 -2α) m (γ + 1 + iν) m (1 -2α) m (2γ + 1) m (γ + 2 -2α + iν) m m! , B = e iπ(γ-iν) Γ(2γ + 1)Γ(1 -2α) Γ(γ + 1 + iν)Γ(γ + 1 -2α -iν) × m (2γ + 1 -2α) m (γ -iν) m (1 -2α) m (2γ + 1) m (γ + 1 -2α -iν) m m! .
Both the 3 F 2 functions above turn to be convergent: we can indeed explicitly check convergence condition (4.2) to give in both cases α > 0. Notice also that assumption α < γ + 1/2 guarantees boundedness on the term Γ(2γ + 1 -2α). To conclude with, we now need to take the real part of (3.11) and show that it is a bounded function of γ. After some simplifications, we have

Re(I 1 (τ )) = 2 1+2α e 2πν π Γ(2γ + 1 -2α)Γ(1 -2α) Γ(2γ + 1) × Re e -i π 2 (2α+1) Γ(γ + 1 + iν)(γ -iν) Γ(γ + 2 -2α -iν) m (2γ + 1 -2α) m (γ + 1 + iν) m (1 -2α) m (2γ + 1) m (γ + 2 -2α + iν) m m! + e i π 2 (2α-1) Γ(γ + 1 -iν) Γ(γ + 1 -2α -iν) m (2γ + 1 -2α) m (γ -iν) m (1 -2α) m (2γ + 1) m (γ + 1 -2α -iν) m m!
where we have written i -(2γ+1-2α) = e -i π 2 (2γ+1-2α) . Using now the basic property of the Pochhammer symbol (x) m = Γ(x + m)/Γ(x) many other terms simplify, leaving us with

Re(I 1 (τ )) = 2 1+2α e 2πν π ×Re e -i π 2 (2α+1) (γ -iν) m Γ(2γ + 1 -2α + m)Γ(γ + 1 + m + iν)Γ(1 -2α + m) Γ(2γ + 1 + m)Γ(γ + 2 -2α + m -iν)m! +e i π 2 (2α-1) (γ -iν) m Γ(2γ + 1 -2α + m)Γ(γ + m -iν)Γ(1 -2α + m) Γ(2γ + 1 + m)Γ(γ + 1 -2α + m + iν)m! .
Stirling formula for Gamma function

Γ(z) = 2π z z e z 1 + O 1 z
shows that both the series above are asymptotic (in γ) to γ -1 , and therefore the proof is concluded.

Remark 3.1. It should be noticed that when α → 1 the term Γ(1 -2α) produces a singularity of order 1, as it is seen by applying well known property Γ(z)Γ(1 -z) = π/ sin(πz). Anyway, it can be checked that this singularity vanishes when taking the real part. Calculations are very easy in the unperturbed case (i.e. when ν = 0): in this case indeed everything is real except for the exponential terms, which become respectively e -i 3 2 π and e i 1 2 π which are purely imaginary. In the perturbed case things are a bit more complicated, but the problem has already been dealt with in literature by some limiting argument (see [START_REF] Biedenharn | Some Exact Radial Integrals for Dirac-Coulomb Functions[END_REF] and references therein). On the other hand, as α approaches 1/2 the singularity produced by the term Γ(1 -2α) is not balanced anymore (exponential terms have a nonzero real part): this rules out the case α = 1/2 and therefore the chance of recovering Strichartz estimates with standard arguments from our proof. The analysis of this problem will be the object of forthcoming works. 4. Appendix: special functions, properties and integrals.

For reader's convenience we include a short appendix in which we recall definitions and basic properties of the special functions coming into play. We refer to [START_REF] Slater | Geralized hypergeometric functions[END_REF] and [START_REF] Erdelyi | Higher transcendental functions[END_REF] for further details and for a deeper insight of the argument. The first thing to investigate is, of course, under which conditions these series converge. As a very general framework, we can summarize the situation as follows

• if p ≤ q the series is convergent for all values of z;

• if p = q + 1 the series is convergent for |z| < 1, and for the special values z = 1 if • if p > q + 1 the series never converges except for z = 0. In this work we are especially interested in the confluent hypergeometric function 1 F 1 (a, c, z), which is thus the case p = q = 1. This series is therefore known to be convergent for every finite value of z (if c is not zero nor a negative integer). Moreover, if a is a negative integer or zero, 1 F 1 reduces to a polynomial of degree |a|. Confluent hypergeometric functions can be thought of as a limiting case of the celebrated Gauss hypergeometric serie, that namely is case p = 2, q = 1 in (4.1), via the relation Among the possible generalizations of hypergeometric function a very useful one is given by the so called Appell functions (see [START_REF] Appel | Sur les fonctions hypergomtriques de plusieurs variables[END_REF]). Historically, the motivation for introducing series of this type is connected to investigating different kind of products of two different Gauss functions. In particular, we are here interested in the following, which is sometimes called Appell function of the second type, It can be seen that convergence condition for F 2 is given by |x| + |y| < 1.

Remark 2 . 3 .

 23 The sign of the Coulomb potential plays in fact a virtually unessential role in what is above, and the wave functions for repulsive fields can be obtained simply by changing the sign of ν in (2.11)-(2.12).

Lemma 3 . 2 .

 32 Let α, β such that |α| + |β| < |h| and d not a negative integer nor zero1 e -ht F (a 1 , b 1 , αt)F (a 2 , b 2 , βt)dt = h -d Γ(d)F 2 (d; a 1 , a 2 ; b 1 , b 2 ;

Definition 4 . 1 .(a 1 )

 411 Given a 1 , . . . , a p and b 1 , . . . , b q complex numbers, we define the generalized hypergeometric function as(4.1) p F q (a 1 , . . . , a p ; b 1 , . . . , b q ; z) := ∞ n=0 n . . . (a p ) n (b 1 ) n . . . (b q ) n z n n!provided the sum of the series is finite, where we are using the Pochhammer symbols (a) 0 = 1, (a) n = a(a + 1) . . . (a + n -1), n ≥ 1.

1 F 1 2 F 1

 1121 (a, c; z) = lim b→∞ (a, b; c; z/b). These functions satisfy the differential equationzu ′′ + (c -z)u ′ -au = 0.The asymptotic behaviour of 1 F 1 can be derived by writing(4.3) 1 F 1 (a, b, z) ∼ = Γ(b) Γ(b -a) (-z) -a G(a, a-b+ 1, z) + Γ(b) Γ(a) e z z a-b G(b-a, 1 -a, z)where G has the asymptotic series representation G(a, b, z)Among the properties that confluent hypergeometric functions satisfy we recall the following one (4.4) F (a, c, z) = e z F (c -a, c, -z).

(4. 5 )

 5 F 2 (d; a 1 , a 2 ; b 1 ; b 2 ; x, y) = m+n (a 1 ) m (a 2 ) n (b 1 ) m (b 2 ) n m!n!x m y n .

Federico Cacciafesta :

 Cacciafesta SAPIENZA -Università di Roma, Dipartimento di Matematica, Piazzale A. Moro 2, I-00185 Roma, Italy E-mail address: cacciafe@mat.uniroma1.it URL: www.mat.uniroma1.it/people/~cacciafesta Éric Séré: CEREMADE (UMR CNRS 7534), Université Paris-Dauphine, Place de Lattre de Tassigny, F-75775 Paris Cedex 16 E-mail address: sere@ceremade.dauphine.fr URL: www.ceremade.dauphine.fr/~sere

  1+α dǫ Remark 2.5. Property (1) allows, by standard arguments, to extend the definition of the operator P k to functions in L 2 .

	Remark 2.6. When summing on k, property (2.19) defines in a standard way
	fractional powers of the operator |D n -ν |x| |, which are used in the statement of Theorem (1.1).

Acknowledgments. The first named author is supported by the FIRB 2012 "Dispersive dynamics, Fourier analysis and variational methods" funded by MIUR