Xiaoying Dai

Yvon Maday

STABLE PARAREAL IN TIME METHOD FOR FIRST AND SECOND ORDER HYPERBOLIC SYSTEM *

Keywords: parareal in time algorithm, parallelisation, time discretization, evolution equations, hyperbolic system, wave equation AMS subject classifications. 65L05, 65P10, 65Y05

The parareal in time algorithm allows to perform parallel simulations of time dependent problems. This algorithm has been implemented on many types of time dependent problems with some success. Recent contributions have allowed to extend the domain of application of the parareal in time algorithm so as to handle long time simulations of Hamiltonian systems. This improvement has managed to avoid the fatal large lack of accuracy of the plain parareal in time algorithm consequence of the fact that the plain parareal in time algorithm does not conserve invariants. A somehow similar difficulty occurs for problems where the solution lacks regularity, either initially or in the evolution, like for the solution to hyperbolic system of conservation laws. In this paper we identify the problem of lack of stability of the parareal in time algorithm and propose a simple way to cure it. The new method is used to solve a linear wave equation and a non linear Burger's equation, the results illustrate the stability of this variant of the parareal in time algorithm.

1. Introduction. Parallel in time algorithms represent a new competitive way of using the ever increasing number of cores in today's supercomputers platforms in cases where classical domain decomposition methods either are of no relevance (for systems of differential equations) or in cases where the scalability of the sole domain decomposition (for evolution PDE's) gets to saturation. This kind of algorithm in the 4 th direction has received too little attention in the literature in the past due to the inherent sequential nature of the time direction: it seems indeed difficult to simulate what will arrive in the far future, e.g. next week, without knowing in details what will arrive in the close future, e.g. tomorrow and each day of the current week. Due to the paramount importance of finding new ways for parallelisation, going together with the ever increasing number of processors present in the new generations of architectures, algorithms such as parareal in time, parallel deferred corrections or parallel exponential integrators have become an active research topic showing promising applicability for large-scale computations and perhaps one way to go to exascale computations. We refer e.g. to the book of K. Burrage [START_REF] Burrage | Parallel and sequential methods for ordinary differential equations[END_REF] for a synthetic approach on the subject (see also [START_REF] Burrage | Parallel methods for ODEs[END_REF]). In this book, the various techniques of parallel in time algorithms are classified into three categories: (i) parallelism across the system, (ii) parallelism across the method and (iii) parallelism across the time.

The parareal in time algorithm pioneered in [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF] extended later in [START_REF] Baffico | Parallel in time molecular dynamics simulations[END_REF][START_REF] Bal | A "parareal" time discretization for non-linear PDE's with application to the pricing of an American put[END_REF][START_REF] Maday | A parareal in time procedure for the control of partial differential equations C[END_REF] under a form much better tuned to the treatment of nonlinear problems enters in the third category above (see [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF] for a first algorithm in this category). It combines two propagators, one fine and one coarse, in a predictor-corrector manner that allows to use the fine propagator in parallel on the various processors, while the coarse propagator is solved in a sequential manner. The algorithm provides good, provable convergence, for diffusive PDEs such as parabolic type problem. However, several studies have shown some instabilities for hyperbolic equations [START_REF] Farhat | Time-decomposed parallel time-integratorspart I: theory and feasibility studies for fluid, structure, and fluidstructure applications[END_REF][START_REF] Fischer | A parareal in time semi-implicit approximation of the navier-stokes equations[END_REF][START_REF] Mercerat | Application of the parareal Algorithm for Acoustic Wave Propagation[END_REF][START_REF] Chouly | An Enhanced parareal Algorithm for Partitioned Parabolic-Hyperbolic Coupling[END_REF] and the numerical analysis performed on simple examples or partial differential equations in [START_REF] Bal | A "parareal" time discretization for non-linear PDE's with application to the pricing of an American put[END_REF][START_REF]On the convergence and the stability of the parareal algorithm to solve partial differential equations[END_REF][START_REF] Staff | Stability of the parareal Algorithm In Domain decomposition methods in science and engineering[END_REF] allows us to understand this behavior.

A challenging field of application for parallel in time propagators is for long time simulation of Hamiltonian systems. Here also the plain method lacks the geometric properties that are essential to guarantee the quality of the approximation. Let us remind that these geometric properties lead to the convergence of the algorithm towards the solution of the original systems on long time simulations. These properties -either the symplecticity or the symmetry of the scheme -are not preserved for the parareal scheme even if the basics coarse and fine propagators are symplectic or symmetric. A cure to this lack of geometric properties has been proposed in [START_REF] Dai | parareal algorithms for Hamiltonian systems[END_REF], based on two ingredients: the first one is the symmetrization of the parareal in time algorithm that leads to a new multistep schemes. This ingredient alone is not sufficient though to provide a correct algorithm for Hamiltonian systems. Indeed, some resonances are artificially introduced by the parareal strategy itself as is explained in [START_REF] Dai | parareal algorithms for Hamiltonian systems[END_REF]; they prevent the symmetric variant of behaving well. The second ingredient that has been introduced in [START_REF] Dai | parareal algorithms for Hamiltonian systems[END_REF] consists simply, after each correction of the predictor corrector scheme, to project the solution coming out the parareal in time algorithm, on the manifold expressing the preservation of some basic invariant quantities of the Hamiltonian system. This is quite simple to implement, at least in its basic formulation (there exists a symmetric version based on a more involved implicit resolution of projection on the manifold).

Following this second idea advocated to stabilize the parareal in time algorithm for Hamiltonian system, we propose in this paper an extension of it allowing to cure the lack of stability that occurs when solving hyperbolic type problem. Another strategy, based on providing a much more involved correction procedure re-utilizing previously computed information has been proposed in [START_REF] Farhat | Time-parallel implicit integrators for the near-real-time prediction of linear structural dynamic responses[END_REF][START_REF] Cortial | A time-parallel implicit method for accelerating the solution of nonlinear structural dynamics problemsInternational[END_REF] and already improved the performance of the parareal in time algorithm for second-order hyperbolic systems. The method proposed here seems much more simple to implement and more general in scope and applications.

The remainder of this paper is organized as follows. In Section 2, we remind the basics on the plain parareal in time algorithm including the current state of the art on the stability and convergence analysis of the algorithm; we use this to introduce the proposed stable variant of the parareal scheme. In Section 3 we address the case of the wave equation, a simple second-order hyperbolic equation. The correction we propose is discussed in full details and the numerical simulation illustrates the good properties of the variant. In Section 4 we address the case of the Burgers' equation, with small viscosity, here again the variant performs well. In Section 5 we present some elements of numerical analysis proving the stability of the proposed variant of the parareal scheme and finally some conclusions are drawn in Section 6.

2. Some elements on the parareal in time algorithm. where y(t) takes its values in a Banach space B and A is a partial differential operator defined over B. We assume here that this problem is well posed, at least for some time: the solution y(t) exists for t ∈ [0, T ). More precisely, we assume that there exists a propagator E such that y(t) = E t 0 (y 0 ), named exact propagator. We note that, with obvious notations, we have

y(t) = E t s (y(s)), ∀s, t, 0 ≤ s ≤ t ≤ T. (2.2)
Let us assume now that we want to simulate this PDE on the interval [0, T ), we thus introduce a discrete propagator F that is assumed to be a fine and precise approximation of the exact solution based on a discretization in time only. The propagator F can be based on an approximation using a time discretization and, e.g. the use of an Euler-type or a more involved scheme based on a small enough time step δt.

For instance, let us denote as T 0 = 0, T N = T and T n = n∆T (with ∆T = T N ) special times at which we are interested to consider snapshots of the solution y(T n ), then, with notations similar as above (2.2), we have the approximation

Y n = F Tn T0 (y 0 ), n = 1, ..., N
we also have

Y n+1 = F Tn+1 Tn (Y n ).
In what follows we present the parareal in time algorithm that allows to define a sequence

{Y k } k , Y k ≡ (Y k 1 , Y k 2 , ..., Y k N )
, such that, for any n, Y k n converges to Y n in the space B as k goes to infinity. The definition of this sequence requires the availability of another approximated propagator, denoted as G. This second propagator is assumed to be cheaper than F at the price of some inaccuracy. One can think immediately about having G based on a similar Euler scheme as F with the larger time step dT, δt < < dT ≤ ∆T . But other possibilities are also offered: in addition to the solution scheme in time, a discretization in space has to be performed, for instance, the spacial discretization can be based on a finite element method, the coarse solver G can then be based both on a coarser time step and a coarser spacial discretization, we refer to [START_REF] Fischer | A parareal in time semi-implicit approximation of the navier-stokes equations[END_REF] and [START_REF] Maday | The "parareal in time" algorithm[END_REF] where such a strategy is followed. Another possibility consists in providing F with all the physics of the phenomenon while G will be based on a simplified physics (see eg [START_REF] Maday | parareal in time algorithm for kinetic systems based on model reduction[END_REF] for an example).

Assume F and G are given, we can proceed to the definition of the sequence. Let us initiate Y 0 by

Y 0 n = G Tn T0 (y 0 ), n = 1, ..., N (2.3)
Assume now Y k is known, the iterative process proceeds by setting Y k 0 = y(0), and,

Y k+1 n+1 = G Tn+1 Tn (Y k+1 n ) + F Tn+1 Tn (Y k n ) -G Tn+1 Tn (Y k n ), n = 1, ..., N. (2.4)
We first comment on the parallelization of the above algorithm.

(Y k n ) -G Tn+1
Tn (Y k n ) can then be built, for any n and the step can proceed. We finally solve the predictor based on the coarse solver G Tn+1 Tn (Y k+1 n ) -sequentially -and immediately correct it by adding the previously computed corrections until convergence. Note that in this process, the fine solutions are only used on small propagation intervals and can be performed in parallel.

2.2. Numerical analysis of the parareal in time algorithm. We review in this subsection some elements on the stability and convergence of the parareal in time algorithm that have been first presented in [START_REF]On the convergence and the stability of the parareal algorithm to solve partial differential equations[END_REF]. Let us assume that there exists a sequence of Banach spaces B j such that B j+1 ⊂ B j . It is also assumed that the problem 2.1 is stable in these spaces, in the sense that there exists a constant C > 0, such that, for all j, 1 ≤ j ≤ J

∥E t 0 (y 0 )∥ Bj ≤ C∥y 0 ∥ Bj (2.5)
The algorithm converges, under the following natural hypothesis: denoting by δF (resp. δG) the difference δF t s = E t s -F t s (resp. δG t s = E t s -G t s ), we assume that there exists a constant C > 0, such that, for all j, 1 ≤ j ≤ J ∀t ≥ 0, ∀τ ≥ 0, ∥G t+τ t (x) -G t+τ t (y)∥ Bj ≤ (1 + Cτ )(∥x -y∥ Bj ).

(2.6)

In addition let us assume that there exists a constant C > 0, such that, for all j,

1 ≤ j ≤ J ∀t ≥ 0, ∀τ ≥ 0, ∥δF t+τ t (x) -δF t+τ t (y)∥ Bj ≤ Cτ η∥x -y∥ Bj+1 , (2.7
)

∥δG t+τ t (x) -δG t+τ t (y)∥ Bj ≤ Cτ ε∥x -y∥ Bj+1 , (2.8)
where η and ε are small quantities, typically, for the Euler scheme these are η = δt and ε = ∆T and C is a generic constant, independent of the quantities on the right hand side of the previous inequalities and on j ≤ J, but may depend on J. Remark 1. For Euler schemes based on a time step δt for F and ∆T G, hypothesis (2.7) and (2.8) hold with η = δt and ε = ∆T respectively.

The convergence is stated in the following (see theorem 1 in [START_REF]On the convergence and the stability of the parareal algorithm to solve partial differential equations[END_REF]) Theorem 2.1. Assume that the discrete propagators F and G satisfy (2.6), (2.7) and (2.8), then the error between the exact solution and the solution provided by the parareal scheme (2.4) satisfies for any k ≥ 0:

∀n ≤ N, ∥Y k n -y(T n )∥ B0 ≤ C(ε k+1 + η)(1 + ∥y 0 ∥ B k ).
(2.9)

Proof: The proof is done by induction over j ≥ 0 for the statement

∀i ≥ 0, ∥Y j n -y(T n )∥ Bi ≤ C(ε j+1 + η)(1 + ∥y 0 ∥ Bi+j+1 ), (2.10)
which is obvious for j = 0 since it corresponds to the analysis of the plain coarse operator. Let us assume that it is true for j and let us prove it for j + 1. From the definition of Y j+1 n+1 , we write

Y j+1 n+1 -y(T n+1 ) = G Tn+1 Tn (Y j+1 n ) + F Tn+1 Tn (Y j n ) -G Tn+1 Tn (Y j n ) -E Tn+1 Tn [y(T n )] = G Tn+1 Tn (Y j+1 n ) -G Tn+1 Tn (y(T n )) +[E Tn+1 Tn -G Tn+1 Tn ](Y j n ) -[E Tn+1 Tn -G Tn+1 Tn ](y(T n )) -[E Tn+1 Tn -F Tn+1 Tn ](Y j n ) + [E Tn+1 Tn -F Tn+1 Tn ](y(T n )) -[E Tn+1 Tn -F Tn+1 Tn ](y(T n ))
or again with the notations introduced above

Y j+1 n+1 -y(T n+1 ) = G Tn+1 Tn (Y j+1 n ) -G Tn+1 Tn (y(T n )) +δG Tn+1 Tn (Y j n ) -δG Tn+1 Tn (y(T n )) -δF Tn+1 Tn (Y j n ) + δF Tn+1 Tn (y(T n )) -[E Tn+1 Tn -F Tn+1 Tn ](y(T n ))
hence by triangular inequality

∥Y j+1 n+1 -y(T n+1 )∥ Bi ≤ ∥G Tn+1 Tn (Y j+1 n ) -G Tn+1 Tn (y(T n ))∥ Bi +∥δG Tn+1 Tn (Y j n ) -δG Tn+1 Tn (y(T n ))∥ Bi +∥δF Tn+1 Tn (Y j n ) + δF Tn+1 Tn (y(T n ))∥ Bi +∥[E Tn+1 Tn -F Tn+1 Tn
](y(T n ))∥ Bi that from hypothesis (2.6), (2.7) and (2.8) can be upper bounded as follows

∥Y j+1 n+1 -y(T n+1 )∥ Bi ≤ (1 + C∆T )∥Y j+1 n -y(T n )∥ Bi +C∆T ε∥Y j n -y(T n )∥ Bi+1 +C∆T η∥Y j n -y(T n )∥ Bi+1 + C∆T η∥y(T n )∥ Bi+1
the stability hypothesis (2.5) together with the induction hypothesis (2.10) allow to get

∥Y j+1 n+1 -y(T n+1 )∥ Bi ≤ (1 + C∆T )∥Y j+1 n -y(T n )∥ Bi +C∆T (ε + η)(ε j+1 + η)∥y 0 ∥ Bi+j+2 + C∆T η∥y 0 ∥ Bi+1 ≤ (1 + C∆T )∥Y j+1 n -y(T n )∥ Bi +C∆T (ε j+2 + η)∥y 0 ∥ Bi+j+2
from which we get (2.10) for j + 1 from the discrete Gronwall lemma, and the proof is complete.

It is interesting to notice two things at this level. First the constant C in (2.9) depends on k and T and this convergence result can be polluted on a long time range by a constant C that becomes too large. Second, the convergence results implies, of course, stability. Hence, provided that the initial condition is regular enough and remains regular the parareal in time algorithm is stable.

The previous remark explains actually the problem of the use of the parareal scheme when non dissipative PDE or PDE with small dissipation are simulated. This is the case for small viscosity parabolic problem as is explained in the next subsection. This is also the case for hyperbolic systems. We can illustrate this statement on two equations: the periodic wave equation and the periodic Burger's equation for which, starting with regular enough initial condition solves or corrects the instability of the parareal in time algorithm in this situation. This allows to rectify somehow the statements that were done before in e.g. [START_REF] Cortial | A time-parallel implicit method for accelerating the solution of nonlinear structural dynamics problemsInternational[END_REF] where the stability problem occurs for second-order hyperbolic problems and not for parabolic and first-order hyperbolic problems. The problem is not really on the type of the equations (and certainly not on the order of the equations) but on the regularity of the solution, including the regularity of the initial condition.

2.3.

Behavior on a parabolic case. We have performed a series of simple test on the linear parabolic equation, provided with periodic boundary conditions

   ∂u ∂t -µ ∂ 2 u ∂x 2 + ∂u ∂x = 0, u(x, 0) = u 0 (x), u(x, t) = u(x + 2π, t), (2.11)
in the case where the viscosity µ is positive yet very small µ = 10 -10 . We have chosen a standard spectral Fourier approximation in space (truncated series with N = 256 and chosen a second order in time implicit Euler scheme for a propagation over [0, T = 2]. The parareal in time algorithm is implemented with ∆T = 2.10 -2 , dT = 10 -2 and δt = 10 -4 .

The simulation are done in order to illustrate the behavior of the scheme as regards the regularity of the solution. We have thus chosen different types of initial conditions with different regularities by choosing u 0 (x) =

∑ ûℓ e iℓx with ûℓ =

   0.5i |ℓ| p , if ℓ > 0 0, if ℓ = 0 -0.5i |ℓ| p , if ℓ < 0,
In the case of low regularity: p = 0.5 or p = 1, the relative error after 15 parareal iterations is about 100! There is no blow up, the iterations can continue, the error eventually going to zero. For p = 4 or even better p = 10, the convergence is achieved after very few iterations, as is expected from Theorem 2.1.

These simple computations provide a preliminary illustration that the behavior of the parareal in time algorithm is deeply linked to the regularity of the solution.

3. The periodic wave equation.

The basic Fourier discretization.

We consider the following one dimensional wave equation with periodic boundary conditions: Find u such that

   ∂ 2 u ∂t 2 -c 2 ∂ 2 u ∂x 2 = 0, u(x, 0) = f (x), ∂u ∂t (x, 0) = g(x), u(x, t) = u(x + T, t), ∂u ∂x (x, t) = ∂u ∂x (x + T, t) (3.1)
where f and g are periodic given initial condition, c stands for the speed of the wave, and T ∈ R + represents the period. The weak form of equation (3.1) is as follows: Find u such that ∀t > 0, u(., t

) ∈ S = {v ∈ H 1 (0, T), v(0) = v(T)} ( ∂ 2 u ∂t 2 , v) + c 2 ( ∂u ∂x , ∂v ∂x ) = 0, ∀v ∈ S, (3.2)
here, (u, v) = ∫ T 0 uvdx. We use now the Fourier spectral method to semi discretize in space problem (3.2). That is, we choose

S N = span{e iℓ 2π T x , ℓ = -N, • • • , N }, (3.3)
and approximate u(x, t) by

u N (x, t) = N ∑ ℓ=-N ûℓ (t)e iℓ 2π T x . (3.4)
Note that ûℓ (t) = 1 T (u N (x, t), e iℓ 2π T x ), as results from the equality (e iℓ 2π T x , e ik 2π T x ) = Tδ ℓ,k where δ ℓ,k denotes the standard Kroneker symbol.

The solution u N is defined by a Galerkin approximation so as to satisfy

( ∂ 2 u N ∂t 2 , v N ) + c 2 ( ∂u N ∂x , ∂v N ∂x ) = 0, ∀v N ∈ S N . (3.5)
As is classical, this Galerkin approximation reduces to a system of differential equations in the Fourier coefficients ûk (t). More precisely, let us choose iteratively v N = e ik 2π T x , k ∈ [-N, N ] in (3.5). We then obtain the following system of independent equations

T ∂ 2 ûℓ ∂t 2 (t) + (2π) 2 T ℓ 2 c 2 ûℓ (t) = 0, ∀ℓ ∈ [-N, N ], (3.6)
that can be exactely solved. After resolution, we know all the coefficients ûℓ (t) in a closed form, then from equation (3.4) we get the Galerkin approximation of u(x, t) in S N . Next, using respectively v = ∂u ∂t in (3.2) and v N = ∂u N ∂t in (3.5) allows to derive, for any time t

∥ ∂u ∂t (., t)∥ 2 L 2 (0,T) + c 2 ∥ ∂u ∂x (., t)∥ 2 L 2 (0,T) = ∥ ∂u ∂t (., 0)∥ 2 L 2 (0,T) + c 2 ∥ ∂u ∂x (., 0)∥ 2 L 2 (0,T) (3.7)
and

∥ ∂u N ∂t (., t)∥ 2 L 2 (0,T) + c 2 ∥ ∂u N ∂x (., t)∥ 2 L 2 (0,T) = ∥ ∂u N ∂t (., 0)∥ 2 L 2 (0,T) + c 2 ∥ ∂u N ∂x (., 0)∥ 2 L 2 (0,T)
(3.8) respectively. If we now introduce the following Hamiltonian

H(u N )(t) = ∥ ∂u N ∂t (., t)∥ 2 L 2 (0,T) + c 2 ∥ ∂u N ∂x (., t)∥ 2 L 2 (0,T) (3.9)
The conservation (3.8) reads for any time t

H(u N )(t) = H(u N )(0) (3.10)
Of course the numerical discretization in time should preserve this invariant . Remark 2. As is often the case for Hamiltonian system, this basic invariant (which is here the total energy for the wave equation) is not the only one available, for instance it is immediately seen on (3.6) that by defining

h ℓ (v) = | ∂v ∂t | 2 + ( 2π T ) 2 ℓ 2 c 2 |v| 2 , (3.11)
for any ℓ, ℓ = -N, .., N , we have

∀t, h ℓ (û ℓ )(t) = h ℓ (û ℓ )(0) (3.12)
which provides 2N + 1 invariant quantities (we notice that H = ∑ N ℓ=-N h ℓ ). 3.2. Illustration of the stability and instabilities of the plain parareal in time algorithm. In all what follows, the discrete solutions will always, for any time t, belong to S N and will result from a Galerkin process. We will thus denote by u n the approximation (in time) of the value of u N (T n ) (hence the full approximation of u(T n )). In order to test the parareal in time algorithm on this simple second order equation, we further discretize the semi discrete problem (3.5) by adding a time discretization, that is based, both for the coarse G ∆T and the fine F ∆T propagators, on the velocity Verlet scheme which is a second order in time, stable and consistent with the preservation of invariants (the wave equation being an Hamiltonian system, the symplectic nature of the velocity Verlet scheme makes it a method of choice).

The plain parareal scheme reads:

{ u 0 n+1 = G ∆T (u 0 n ) u k+1 n+1 = G ∆T (u k+1 n ) + F ∆T (u k n ) -G ∆T (u k n ). (3.13)
We perform simulations with the discretization parameters equal to 1. T = 100, ∆T = 2.10 -1 , dT = 4.10 -3 , δt = 2.10 -6 . 2. N = 30. for different initial conditions of different regularities:

u 0 (x) = ∑ ℓ∈Z ûℓ e iℓ 2π T x with ûℓ = { 1 |ℓ| p , if ℓ ̸ = 0 0, if ℓ = 0,
with p p = 1, 5 or 10. The cases with p = 1 and p = 5 diverges while the case with p = 10 converges.

3.3.

Results with the plain parareal in time algorithm on a realistic test problem. In this example, the excitation is provided by a Ricker pulse, which is analytically defined as follows 1. T = 5000 m, c = 2000 m/s 2. and the initial conditions are

f (x) = 0, g(x) = (1-2(f s π x-xs c ) 2 ) exp(-(f s π x-xs c ) 2
). Here, f s = 2.5 Hz, is the frequency and x s = 2500 m is the position of the vibration source.

The numerical results are obtained with a coarse propagator G and fine propagator F that differ only with time step dT and δt, respectively :

1. T = 100, ∆T = 2.10 -1 , dT = 4.10 -3 , δt = 2.10 -6 . 2. N = 30. The plot showing the evolution in time of the solution obtained after 15 iterations at different positions in the periodic space [0, , T] on figure 3.1 reveals an instability visible starting at time t = 80. The situation is even worse after 25 iterations. Actually this is one case that illustrates the bad behavior of the parareal in time algorithm for this hyperbolic equation. Note however that for many other choices of the discretization parameters δt or dT or even N but the very same initial condition, the simulation is perfectly fine, stable and convergent. The instability is thus not systematic, and changing slightly the discretization parameters, is often enough to make the scheme work. The plain parareal in time algorithm does lack of robustness, and this is what we want to correct in what follows. The exact solution shown on the next plot (figure 3.2) illustrates the differences. The plots of the error between the exact solution and the solution obtained after k iterations, with k = 1, .., 15 (see figures 3.3 and 3.4) illustrate clearly that the solution obtained by the parareal in time algorithm 1. converges for short time, e.g. only 7 iterations are sufficient to recover the solution obtained with the sole sequential fine propagator at time equal to 10. 2. does not converge for longer propagations, here for a time equal to 100 3. During the first 5 to 6 iterations the error of the parareal in time algorithm is dominated by the error of the coarse scheme, but after k = 6, the error increases largely to values that are order 1 and bigger 4. After 15 iterations of the parareal in time algorithm the parareal solution is computed with an acceptable error until time t ≃ 30 only. Taking into account the variants of the parareal in time algorithm introduced in [START_REF] Dai | parareal algorithms for Hamiltonian systems[END_REF] to simulate Hamiltonian systems we infer that the bad behavior of the parareal scheme comes from the non preservation of the invariant quantities. Let us for example observe the evolution of the value of the quantity H defined in (3.9) The figures 3.5 and 3.6 illustrate the lack of conservation, that is coherent with the divergence of the scheme in this case.

Finally we want to notice that we have performed the same simulation with a Ricker pulse of varying force (by modifying the value of f s ). The instability starts much earlier for larger f s confirming the relation between the instability and the singularity of the solution.

Parareal algorithm with projection.

Let us define the manifold M corresponding to the set of all functions v such that H(v) = H(u(0)). The Standard Projection Method on M (see chapter 4 of [START_REF] Hairer | Geometric numerical integration: structurepreserving algorithms for ordinary differential equations[END_REF]) can be written as follows: For a given ũ, solve the system in (u, λ)

u = u + λ∇H( u), where λ ∈ R is such that H(u) = H(u(0)), i.e. H( u + λ∇H( u)) = H(u(0)). (3.14)
We then set

u = π M ( u) (3.15)
In the general case, equation (3.14) can only be solved by some iterative method, for example, Newton iterations, but in the current case, using the fact that H is quadratic, we notice and take into account that equation (3.14) is a second order polynomial in λ that can be solved exactly.

Then, combining the plain parareal scheme (3.13) with the projection method above, we obtain the parareal method with projection to M:

   u 0 n+1 = G ∆T (u 0 n ) ũk+1 n+1 = G ∆T (u k+1 n ) + F ∆T (u k n ) -G ∆T (u k n ), u k+1 n+1 = π M (ũ k+1 n+1 ).
(3.16)

We now show the results obtained by parareal method with projection. We first note on figure 3.7 that the global solution after k = 15 iterations does not seem to present the instabilities of the plain parareal in time algorithm. The same conclusion holds after 25 iterations. This is confirmed by the figures 3.8 and 3.9 where we note an improvement in the error between the exact solution and the discrete one as the iterations evolve. The error is always at worse of the same order as the error of the coarse propagator.

Of course this is far from being what is expected, but the conservation of H does improve the behavior of the algorithm. The following figures 3.10 and 3.11 illustrate the good preservation of H along the trajectory for the parareal method with projection. Although the solution does not diverge as in the plain parareal case, it does not converge as we would like.

In the results above, when doing the projection, we view all the Fourier modes as a whole system, and treat them altogether, with a unique Lagrange multiplier. Remember that each Fourier coefficient individually has an energy that is preserved as is explained in Remark 2. Instead of projecting on the manifold M, where only one invariant is taken into account, we can have a manifold defined with more invariant properties. We have defined only two invariants based on two groups in the Fourier coefficients, those corresponding to high wave numbers and those corresponding to low wave numbers. That is, we partition all the Fourier modes into two groups, one contains the modes |k| ≥ K 0 , the other one contains all the modes |k| < K 0 , we then perform a projection for each group separately to make the energy of each group conserved.

The motivation comes from the fact that, when we discretize the solution by Fourier series, the contribution of the high frequency modes is much smaller than that of the low frequency modes. If we use the same Lagrange multiplier λ to correct all the modes, the contribution of the high frequency modes may be to magnify them and pollute the solution.

So, let us introduce two Lagrange multipliers λ 1 and λ 2 , one is for the high frequency modes, the other is for the low frequency modes only. The method is not more complex to implement. The numerical simulation still with the same N = 30 is based on K 0 = 20; the results are much improved as can be seen on figures 3.12, 3.13 and 3.14.

The solution itself, after 15 iterations, is in excellent accordance with the very fine approximation as can be seen on figure 3.15.

For the linear second order wave equation, the preservation of the invariant quantities based on the decomposition of the coefficients representing the solution into two groups: the low modes and high frequency modes, is thus a sufficiently mild ingredi- ent to improve the behavior of the parareal in time algorithm. In the next section we shall extend this variant of the parareal in time algorithm to a nonlinear PDE.

Note that, if when the method converges, i.e. ∀n = 1, .., N , u k n -→ u ∞ n , this limit solution satisfies thanks to (3.16)

u ∞ n+1 = π M (F ∆T (u ∞ n ))
that, here, due to the fact that the propagator F ∆T is symplectic and thus preserves the Hamiltonian, reads

u ∞ n+1 = F ∆T (u ∞ n
) that proves that the method converges towards the fine solution, as expected.

∂u ∂t -ν ∂ 2 u ∂x 2 + 1 2 ∂(u 2 ) ∂x = 0, (4.1) u(x, 0) = u 0 (x), u(x + T, t) = u(x, t),
The viscosity ν > 0 is planned to be very small so that the solution will exhibit very sharp gradients.

Taking into account that the periodic frame of the problem allows to use Fourier spectral approximation, we thus consider a semi discretization based on the Galerkin Fourier approximation in S N

( ∂u N ∂t , v N ) + ν( ∂u N ∂x , ∂v N ∂x ) + 1 2 ( ∂(u N ) 2 ∂x , v N ) = 0, ∀v N ∈ S N . (4.2)
We then choose to use a second order symmetrized splitting method based on the splitting between the diffusion effects (treated in the Fourier space exactly) and the nonlinear convection effects (treated in the physical space). The only difference between the coarse propagator G and fine propagator F are relative to the time steps dT and fine time step δt, respectively. Denoting by δ either dT or δt, the splitting schemes at each iteration proceeds as follows:

1. Solve the following problem in Fourier space in the interval [τ n = nδ, τ

n+ 1 2 = τ n + δ 2 ] ( ∂u N ∂t , v N ) + ν( ∂u N ∂x , ∂v N ∂x ) = 0, ∀v N ∈ S N (4.3) u N (τ n ) = u N,n , (4.4)
and denote ũN,n+ 1 2 = u N (τ n+ 1 2 ). Note that, as in the discretization of the wave equation, by choosing iteratively v N = e ik 2π T x , k ∈ [-N, N ], this results in a system of differential equations in the Fourier coefficients.

Solve the following problem in

[τ n , τ n+1 ] ( ∂u N ∂t , v N ) + 1 2 ( ∂(u N ) 2 ∂x , v N ) = 0, ∀v N ∈ S N (4.5) u N (τ n ) = ũN,n+ 1 2 , (4.6)
and denote u N (τ n+1 ) as u N,n+ 1 2 . Note that here, the solution procedure cannot be exact since the problem is non linear: the discretization in time is based on a 3 rd order Runge Kutta method and an exact evaluation of the Fourier coefficients of ∂(uN ) 2 ∂x by a Discrete Fourier Method based on S 2N .

Solve the following problem in Fourier space in [τ

n+ 1 2 , τ n+1 ] ( ∂u N ∂t , v N ) + ν( ∂u N ∂x , ∂v N ∂x ) = 0, ∀v N ∈ S N (4.7) u N (τ n+ 1 2 ) = u N,n+ 1 2 . (4.8)
4.2. Results for the plain parareal in time algorithm. We have performed 3 tests for the plain algorithm. First we note that, provided that the viscosity is large enough, the parareal in time algorithm is working fine. Indeed in case referred later as Case 1, we have chosen Case 1:

1. µ = 1e -2, 2. N c = 2 8 , N f = 2 8 , 3. T = 2, ∆T = 2.10 -2 , dT = 1.10 -2 , δt = 1.10 -4
The evolution is presented on figure 4.1, the relative error in the solution with respect to a solution computed with a much smaller time step is presented on figure 4.2 showing that after 3 parareal iterations, the same convergence as for the fine sequential simulation is achieved. For the second test, we have diminished the value of the viscosity. This case, referred to as Case 2, corresponds to the choices Case 2: For this case, the plain parareal method is not stable, as can be seen on the figure 4.3, after 3 iterations instabilities occur soon after time 1 where the solution of the inviscid Burgers equation becomes discontinuous.

1. µ = 1e -3; 2. N c = 2 9 , N f = 2 9 , 3. T = 2, ∆T = 2.
It should be noted that, as explained in [START_REF] Maday | The "parareal in time" algorithm[END_REF], that using a coarser spacial representation helps in the stabilization of the process. Indeed, Case 3 below provides a stable parareal approximation Case 3:

1. µ = 1e -3; 2. N c = 2 8 , N f = 2 9 , 3. T = 2, ∆T = 2.10 -2 , dT = 1.10 -2 , δt = 1.10 -5 .
As can be seen on the series of figures 4.4, 4.5 and 4.6, the method is stable, converges but very slowly. In the next sections, we shall improve the two last simulations by extending the method presented in the previous section.

4.3.

Results for the plain parareal in time algorithm with projection alone. To start with, let us notice that, due to dissipation, there is no quantity that is preserved along a trajectory, nevertheless, by using u N as a test function in (4.2), we derive that

∥u N (t)∥ 2 L 2 (0,T) + 2ν ∫ t 0 [ ∂u N ∂x ] 2 = ∥u N (0)∥ 2 L 2 (0,T) ,
and actually for any two times

T n+1 > T n ≥ 0 ∥u N (T n+1 )∥ 2 L 2 (0,T) + 2ν ∫ Tn+1 Tn [ ∂u N ∂x ] 2 = ∥u N (T n )∥ 2 L 2 (0,T) . (4.9)
In the parareal in time algorithm, starting from u k+1 n (remember that this represents the Galerkin approximation that should be denoted as u k+1 N,n ), we expect that Since we do not want to compute F ∆T [u k+1 n ] at this stage, we approximate ∥F ∆T [u k+1 n ]∥ L 2 (0,T) by

∥u k+1 n+1 ∥ L 2 (0,T) = ∥F ∆T [u k+1 n ]∥ L 2 (0,T) (4
∥F ∆T [u k+1 n ]∥ L 2 (0,T) ≃ ∥F ∆T [u k n ]∥ L 2 (0,T) ∥u k n ∥ L 2 (0,T) ∥u k+1 n ∥ L 2 (0,T) . (4.11)
What can thus be done is to consider the manifold M k+1 n corresponding to the set of all functions v such that ∥v∥ L 2 (0,T) =

∥F∆T [u k n ]∥ L 2 (0,T) ∥u k n ∥ L 2 (0,T)
∥u k+1 n ∥ L 2 (0,T) Second, we have noticed in the analysis of the wave equation that treating all the frequencies with a unique Lagrange multiplier is not effective enough. We thus partition the frequency space of functions in S N into two (or more if necessary) groups G i , i = 1, 2 (or i = 1, .., I) defined by

G i = span{e iℓ 2π T x , N - i ≤ |ℓ| < N + i }, (4.12)
where

L - 1 = 0, 0 < L + 1 = L - 2 -1, L + 2 = N + 1 (or L - 1 = 0, 0 < L + 1 = L - 2 -1, L - i < L + i = L - i+1 + 1, L + I = N + 1 if
more groups are around). The L 2 -projection on group G i is denoted as P i . In this case we consider the manifold [M I ] k+1 n defined as being the intersection

[M I ] k+1 n = ∩ I i=1 [m i ] k+1 n (4.13)
where the group manifolds [m i ] k+1 n are defined as follows .

[

m i ] k+1 n = {v ∈ S N , ∥P i (v)∥ L 2 (0,T) = ∥P i [F ∆T [u k n ]]∥ L 2 (0,T) ∥u k n ∥ L 2 (0,T) ∥u k+1 n ∥ L 2 (0,T) } (4.14)
To any given ṽ ∈ S N , we associate the element π [M I ] k+1 n (ṽ) defined by the Standard Projection Method on [M I ] k+1 n , that involves the determination of I Lagrange multipliers.

With this manifold the parareal scheme with projection reads

   u 0 n+1 = G ∆T (u 0 n ) ũk+1 n+1 = G ∆T (u k+1 n ) + F ∆T (u k n ) -G ∆T (u k n ), u k+1 n+1 = π [M I ] k+1 n (ũ k+1 n+1 ).
(4.15)

The results are as follows (see figures 4.7, 4.8, 4.9). We notice that the parareal in time algorithm is stabilized, but it does not converge much. Increasing the number of groups does not improve the behavior of this variant of the parareal in time algorithm.

Remark 3. If we would know the exact value of the energy of each groups, the projection method would converge. The problem of course is that due to dissipation, and also to nonlinear effect that mixes up the low and high frequencies, the value of the exact energy in each groups is not known. This is the reason why we get into these difficulties.

4.4.

Results for the plain parareal in time algorithm with projection and damping. We notice that, for this non linear problem, the projection allows to stabilize the parareal in time algorithm but a little more is required. We have chosen to damp the high frequencies at each coarse step before proceeding to the next parareal iteration. This means that we introduce a damping function Φ that, to any element v N in S N provides an element w

N = Φ(v N ) in S N , such that ∀ℓ ∈ G 1 , ŵℓ = vℓ (4.16) ∀ℓ ∈ G i , i = 2, ..I, ŵℓ = f (i)v ℓ (4.17)
where f (i) is chosen such that high frequency parts are more and more filtered. In all the following results, we choose f (i) = 1 i 2 +1.0 , except for 2 groups, were we have chosen scheme with projection and damping reads

       u 0 n+1 = G ∆T (u 0 n ) ũk+1 n+1 = G ∆T (ū k+1 n ) + F ∆T (u k n ) -G ∆T (ū k n ), u k+1 n+1 = π [M I ] k+1 n (ũ k+1 n+1 ) ūk+1 n+1 = Φ(u k+1 n+1 ). (4.18)
The following figures 4.10, 4.11 and 4.11 present the results of this improved parareal method. We notice that the parareal method is now converging.

Remark 4. Note that the strategy used here consisting in adding a little bit of dissipation on the high modes to the parareal scheme goes in the direction suggested in [START_REF]On the convergence and the stability of the parareal algorithm to solve partial differential equations[END_REF] and is coherent with the spectral viscosity proposed in [START_REF] Tadmor | Convergence of Spectral Methods for Nonlinear Conservation Laws[END_REF]. Note also that the solution should be measured after the projection and before the damping. Note finally that the fine propagator acts on the undamped solution so that the dissipativity that is added does not affect the quality of the solution (e.g. as any parareal scheme, after N th iterations, the solution is the same as the one of the plain fine propagator) but improves the quality of the scheme.

The projection with more than two groups performs slightly better, see e.g. figure 4.13 where only the results for iterations between 7 and 15 are presented. It may be interesting to keep this in mind even though the separation in more than two groups is not completely straightforward in the non periodic case where finite difference, finite differences or finite volume techniques will be used.

Preliminary analysis.

In this section, we provide some simple elements of analysis that prove the stabilty of the previous approximation with the parareal method with projection (and damping).

Let us first assume that the fine propagator is L 2 stable in the sense that there This proves the stability by induction ∥u k+1 n+1 ∥ L 2 (0,T) ≤ (1 + α∆T ) n+1 ∥u k+1 0 = u 0 ∥ L 2 (0,T) ≤ e αT ∥u 0 ∥ L 2 (0,T)

(5.3)

6. Concluding remarks. In this paper we have analyzed and proposed a cure to the lack of robustness of the parareal in time algorithm applied to hyperbolic systems or convection diffusion problems with small diffusion. The lack of robustness, leading sometimes to instabilities of the algorithm before convergence is a consequence of the lack of regularity that the solution develops when time runs. The cure is inherited from the extension of the parareal in time algorithm for Hamiltonian system and consists in projecting the solution over an energy manifold defined on the fly, and improved iteratively, where the exact solution lies. This simple procedure allows to stabilizes the algorithm that now converges after some iterations. It should be noted that, similarly as for other types of equations, the convergence of the parareal in time algorithm allows to speed up the solution procedure but is not optimal in terms of parallel efficiency since, roughly speaking, a convergence that is achieved after 10 or 15 iterations leads to a parallel efficientcy of 1 10 or 1 15 at best, which is far from optimal. It is not the purpose of this paper to provide a full optimal implementation of the approach. Let us just notice that the combination of the parareal in time algorithm with standard domain decomposition methods is a way to improve the efficiency and use more of the processors that are available on current architectures than the one that a plain domain decomposition technique is able to optimally use. We refer to [START_REF] Guetat | Mthode de parallélisation en temps : Application aux méthodes de décomposition de domaine[END_REF] for more on this matter.

2. 1 .

 1 Basics on the parareal in time algorithm. Let us consider a possibly non-linear time dependent partial differential equation of the form ∂y ∂t + A(t; y) = 0, y(0) = y 0 (2.1)

Fig. 3 . 1 .

 31 Fig. 3.1. Vibration at different positions obtained by the parareal method(δt = 2.10 -6 , dT = 4.10 -3 , ∆T = 2.10 -1 ) revealing the instabilities.

Fig. 3 . 2 .Fig. 3 . 3 .

 3233 Fig. 3.2. Vibration at different positions obtained by the fine sequential method(δt = 2.10 -6 )

Fig. 3 . 4 .

 34 Fig. 3.4. Error on the solution obtained by the parareal method(δt = 2.10 -6 , dT = 4.10 -3 , ∆T = 2.10 -1 )

Fig. 3 . 5 .Fig. 3 . 6 .

 3536 Fig. 3.5. Error on the energy obtained by the parareal method(δt = 2.10 -6 , dT = 4.10 -3 , ∆T = 2.10 -1 )

Fig. 3 . 7 .

 37 Fig. 3.7. Vibration at different positions obtained by the parareal method with projection to the manifold M (δt = 2.10 -6 , dT = 4.10 -3 , ∆T = 2.10 -1 )

Fig. 3 . 8 .Fig. 3 . 9 .

 3839 Fig. 3.8. Error on the solution obtained by the parareal method with projection to the manifold M (δt = 2.10 -6 , dT = 4.10 -3 , ∆T = 2.10 -1 )

Fig. 3 . 10 .Fig. 3 . 11 .

 310311 Fig. 3.10. Error on the energy obtained by parareal method with projection to the manifold M (δt = 2.10 -6 , dT = 4.10 -3 , ∆T = 2.10 -1 )

Fig. 3 . 12 .Fig. 3 . 13 .

 312313 Fig.3.12. Error on the solution obtained by the parareal with projection on two groups (δt = 2.10 -6 , dT = 4.10 -3 , ∆T = 2.10 -1 )

Fig. 3 . 14 .Fig. 3 . 15 .

 314315 Fig. 3.14. Error on the solution obtained by the parareal with projection on two groups (δt = 2.10 -6 , dT = 4.10 -3 , ∆T = 2.10 -1 )

4 .

 4 Burgers equation. 4.1. The total Fourier/time splitting discretization. Let us consider in this section the following periodic viscous Burgers' equation in one dimension

Fig. 4 . 1 .

 41 Fig. 4.1. Case 1: solution obtained by the parareal method(δt = 1e -4, dT = 1e -2, DT = 2e -2)

10 - 2 Fig. 4 . 2 .Fig. 4 . 3 .

 1024243 Fig. 4.2. Case 1: relative error on the solution obtained by the parareal method(δt = 1e-4, dT = 1e -2, DT = 2e -2)

Fig. 4 . 4 .

 44 Fig. 4.4. Case 2-1: relative error of the solution obtained by the parareal method(δt = 1e -5, dT = 1e -2, DT = 2e -2)

Fig. 4 . 5 .Fig. 4 . 6 .

 4546 Fig. 4.5. Case 2-1: relative error of the solution obtained by the parareal method(δt = 1e -5, dT = 1e -2, DT = 2e -2)

Fig. 4 . 7 .Fig. 4 . 8 .Fig. 4 . 9 .

 474849 Fig. 4.7. Case 2-2: relative error of the solution obtained by the method parareal with projection to 2 groups(δt = 1e -5, dT = 1e -2, DT = 2e -2)

Fig. 4 . 10 .Fig. 4 . 11 .Fig. 4 . 12 .Fig. 4 . 13 .Fig. 4 . 14 .≤ ( 1 + 2 L 2 (

 410411412413414122 Fig. 4.10. Case 2-2: relative error of the solution obtained by the method parareal with projection to 2 groups with damping(δt = 1e -5, dT = 1e -2, DT = 2e -2)

Acknowledgments. This work has been supported in part by the Agence Nationale de la Recherche, under the grant ANR-06-CIS6-007 (PITAC).

* This work was supported by ANR-06-CIS6-007.


A Heusler
email: a.heusler@mpi-hd.mpg.de

T Faestermann

R Hertenberger

R Krücken

H-F Wirth

P Von Brentano

Structure of states in 208 Pb with major components of the neutron particle-hole configuration ν(d +1 5/2 p -1 3/2 )

Keywords: numbers: 21, 10, -k 21, 10, Hw 24, 30, Gd 27, 80, +w

Using the 208 Pb(p, p ) reaction via isobaric analog resonances in 209 Bi, a dozen states in the doubly magic nucleus 208 Pb with major strength of the neutron particle-hole configuration ν(d 5/2 +1 p 3/2 -1 ) are investigated. Among them, two doublets at Ex = 5.81, 6.01 MeV are resolved by comparing the mean cross sections on the d 5/2 resonance. The 5.81 MeV doublet consists of two states with major configuration ν(d 5/2 +1 p 3/2 -1 ) and spins 2 -and 3 -; the 2 - state is newly identified. The distance between the two states is 0.4 ± 0.2 keV. The 6.01 MeV doublet consists of a 4 -state with dominant configuration ν(d 5/2 +1 p 3/2 -1 ) and a 3 -state with a small admixture of ν(d 5/2 +1 p 3/2 -1 ). The distance between the two states is 1.8 ± 0.3 keV. The 5886 state has been newly assigned the spin of 4 -. The spins of seven more states are verified. Almost the complete ν(d 5/2 +1 p 3/2 -1 ) strength is identified in these twelve states. The center of gravity for the configuration ν(d 5/2 +1 p 1/2 -1 ) and ν(d 5/2 +1 p 3/2 -1 ) is in line with the predictions by the shell model without residual interaction to within 10 keV.

Introduction

States in the doubly magic nucleus 208 Pb are of several different types, collective states, pairing and octupole vibrations, as well as one-particle one-hole and two-particle twohole states. Most states at E x < 6. 5 MeV are described as one-particle one-hole states. Studying them deserves special interest.

Seven low-lying states in 209 Pb are known to be rather pure single-particle states (g 9/2 , i 11/2 , j 15/2 , d 5/2 , s 1/2 , g 7/2 , and d 3/2 ). By coupling a neutron hole (p 1/2 , p 3/2 , f 5/2 , f 7/2 , h 9/2 , i 13/2 ) with these states, particle-hole states in 208 Pb are created. The neutron pickup reaction on 209 Pb * is capable of exciting these particle-hole states. An equivalent reaction is the inelastic proton scattering via isobaric analog resonances (IAR). One isobaric analog state in 209 Bi corresponds to each particle state in 209 Pb.

The analog state can be created by adjusting the proton energy in the 208 Pb(p, p ) reaction. The opportunity to determine from the angular distributions spectroscopic amplitudes of neutron holes coupled with the neutron particle is an advantage of the (p, p ) reaction via IAR [1].

In this paper, in particular the configuration ν(d 5/2 +1 p 3/2 -1 ) located at E x ≈ 5.9 MeV in 208 Pb is discussed. Only a few of the other 42 neutron particle-hole configurations admix to the d 5/2 p 3/2 states (namely g 9/2 f 7/2 , g 7/2 p 1/2 , d 3/2 p 1/2 , and s 1/2 f 5/2 ), and admixtures of proton particle-hole configurations are negligible.

The shell model without residual interaction (SSM) describes the structure of many states in 208 Pb by particle-hole configurations. Ref. [2], henceforth called Paper S, describes the SSM in detail. The SSM predicts four states with the configuration d 5/2 p 3/2 . However, there are ten more states with the same spins (1 -, 2 -, 3 -, 4 -) predicted to have similar excitation energies. Hence, at least fourteen states are expected to share the major d 5/2 p 3/2 fractions in the region at E x ≈ 5.9 MeV in 208 Pb.

In 1970, by essentially using only data from the 208 Pb(p, p ) reaction via the lowest IAR in 209 Bi, the structure of twenty negative parity states at E x < 4.7 MeV in 208 Pb was determined in a rather complete manner. The residual interaction among particle-hole configurations was deduced with the result that the mean value of an offdiagonal matrix element for unnatural parity is about 50 keV while for natural parity it is about 150 keV [3].

However, it is crucial to correctly identify the states, as well as to obtain firm spin and parity assignments. For the states at E x < 4.5 MeV in 208 Pb this was realized in 1982 [4] and is now confirmed [5]. The complete identification, and spin and parity assignment for states at E x < 4.8 MeV was carried out thirty years later [6], when the Q3D magnetic spectrograph at the Maier-Leibnitz-Laboratorium (Munich) [7,8,9] had matured to yield high linearity [START_REF] Heusler | Maier-Leibnitz Laboratorium[END_REF] and high particle resolution (Paper S).

More than a hundred states are observed in the range of 4.8 < E x < 6.1 MeV [5]. Therefore, many doublets with distances between two states of less than 4 keV are known [2,5,6,[START_REF] Heusler | [END_REF]12,13,14,15]. The main data discussed in this paper (twelve states at 5.75 < E x < 6.05 MeV) has been presented in Refs. [14,15], henceforth called Paper P; Table 2 therein is called Table T.

The theoretical description of the (p, p ) reaction via IAR is briefly described in Appendix A where further references are given. In this paper we essentially rely on relative mean cross sections derived from Table T. Determining spectroscopic amplitudes from cross sections depends strongly on the energy of the scattered protons and their angular momenta (the so-called penetrability). In addition, special assumptions are made. Therefore, required formulas are given in Sec. 3.2.3. An important experimental tool is the computer code GASPAN [16], see also Paper S. It allows spectra of the (p, p ) and (d, p) reactions to be analyzed, using many different methods. Eventually, excitation energies for a limited energy region can be determined with a precision of 0.1 keV.

In this paper, two doublets at E x = 5.81 and 6.01 MeV with distances of 0.4 keV and 1.8 keV are disentangled (Secs. 3.3, 4.4, and 4.6). The distribution of almost the complete d 5/2 p 3/2 strength across twelve states has been identified (Sec. 4). The 5812 2 -state has been newly identified, spin of 4 -is newly assigned to the 5886 state. The identification of the 6012 state [5] and its tentative assignment of spin 4 -by Kulleck et al [17] is confirmed.

Model configurations

In 208 Pb at excitation energies E x < 6.4 MeV, a dozen doublets with distances of less than 2 keV are known. Since the determination of excitation energies develops over time, each state identified is given a unique four-digit energy label Ẽx ≈ E x (Paper S).

Most states in 208 Pb with spin and parity I π and excitation energy Ẽx are described by particle-hole configurations in the SSM,

| Ẽx I π = LJ,lj c Ẽx,I π LJ lj |LJ lj , with real amplitudes c, -1 < c Ẽx,I π LJ lj < +1. ( 1
)
The SSM is explained in detail in Paper S. The particles LJ and the holes lj move in orbitals with orbital momenta L, l and spins J, j, respectively; both the neutron particle-hole configurations ν(LJ +1 lj -1 ) and the proton configurations π(LJ +1 lj -1 ) are abbreviated to LJ lj. For brevity, the magnetic quantum numbers of the state and configuration are omitted.

In this paper, we discuss states containing major fractions of the configuration d 5/2 p 3/2 . They are predicted by the SSM at the excitation energy E SSM x = 5.896 MeV (see Table 1 in Paper P). The d 5/2 p 3/2 states with spins 2 -and 4 -are not yet identified and some spin assignments listed in the Nuclear Data Sheets [5] (henceforth called NDS2007) are doubted. In this work, the fragmentation of the d 5/2 p 3/2 strength over twelve states with spins 1 -to 4 -is determined. In these states, components of additional particle-hole configurations are observed, but not discussed in detail.

Experiments

We discuss measurements performed with the Q3D spectrograph at the 14 MV tandem accelerator in Munich for the 208 Pb(p, p ) reaction on all known IAR and especially on top of the d 5/2 IAR; data for the 207 Pb(d, p) reaction is also considered (Paper P).

At excitation energies 4.8 < E x < 6.1 MeV in 208 Pb, the mean distance between states is about 10 keV. The resolution of about 3 keV is hardly sufficient to resolve all states. In general, it is difficult to identify states in a spectrum.

The term "peak" is used to define a visible hump in the raw spectrum of some reaction such as 208 Pb(p, p ).

The term "level" defines the result of the peak fitting using computer codes such as GASPAN [16]. The result of a fit may depend on details such as the handling of secondary peaks produced by contamination reactions such as 14 N(p, p ). The electron knockout reaction, 208 82 Pb (p, p + z e -) 208 82 Pb z+ , z= 0, • • • , 82, enlarges each peak and produces secondary peaks [START_REF] Heusler | In Changing Facets of Nuclear Structure[END_REF]. The term "state" defines a physical state in a nucleus described by some model configurations such as Eq. (1). The correspondence between a level and a single state or an ensemble of states is not straightforward.

Therefore the term "doublet" defines an unresolved or incompletely resolved peak in the spectrum of the reaction considered, produced by a pair of states (similarly "triplet" etc). Finally, a previously assumed doublet may be recognized as a single state.

The term "multiplet" defines a group of states with a certain configuration LJ lj and spins |J -j| ≤ I ≤ J + j.

Data evaluation

For some Ẽx I π state, the measured angular distribution for the 208 Pb(p, p ) reaction via LJ IARs in 209 Bi is fitted by a series of even order Legendre polynomials P K ,

dσ dΩ ( Ẽx , I π , Θ, E p , LJ) = Kmax K=0,2,••• A K ( Ẽx , I π , E p , LJ)P K (Θ). (2)
For natural parity states, a considerable direct (p, p ) component often overwhelms the resonant part for Θ < ∼ 100 • . Hence, if a direct (p, p ) component is present, Eq. ( 2) should be used for scattering angles Θ > ∼ 100 • only. More precise angular distributions (often with K max = 6) are obtained through experiments with semiconductor detectors at scattering angles up to Θ = 173 • [17,[START_REF] Glöckner | [END_REF]; however, not all levels are observed and the resolution was only 18 and 15 keV, respectively.

The isotropic component A 0 defines the mean cross section of the state Ẽx I π on the LJ IAR. The dependence on the beam energy E p is eliminated by applying the inverse Lorentzian factor a Lz [Eq. (A.10)]. Hence, the angle-integrated (mean) cross section on top of the LJ IAR is determined by

A mean 0 ( Ẽx , I π , LJ) = 1 N N i=1 A 0 ( Ẽx , I π , E i p , LJ) a Lz (E i p , LJ) . (3)
The values shown in Tables 1 and2 . The basic assumption of an isolated IAR is taken for granted for the states under discussion. The distance between the d 5/2 and the s 1/2 IARs is 465 keV while the widths of the IARs is about 300 keV [22]. The mean cross section of a state with a d 5/2 lj component is calculated to decrease by a factor of about 5 and 10, respectively, in relation to the proton energy corresponding to the s 1/2 IAR and the g 7/2 +d 3/2 doublet IAR. For the 5.81 and 6.01 MeV doublets, the available excitation functions (Fig. 7 of Ref. [22]) agree with the calculations. An excitation on the d 3/2 IAR is present (see also Table T); however, no significant enhancement on the s 1/2 IAR is observed. Spectra for 5.4 < E x < 6.2 MeV on all IARs are shown in Fig. 1 of Paper P. The 5.81 MeV level is excited on the g 9/2 , d 5/2 and g 7/2 + d 3/2 IARs. Spectra taken on the g 9/2 and d 5/2 IARs for 5.75 < E x < 5.90 MeV are shown in Fig. 4 of Ref. [12]. The 5778 2 -level is strongly excited on the d 5/2 IAR only (Table T).

Mean cross sections.

Should configurations LJ lj with different holes lj contribute to the state, then the strength of each configuration can only be determined with difficulty. The determination of the spin of the states, the configuration strengths in the states and the s.p. widths is highly correlated (Appendix A). However, for states with energies E x ≈ 5.0 and E x ≈ 5.9 MeV strongly excited on the d 5/2 IAR, the configurations d 5/2 p 1/2 and d 5/2 p 3/2 are assumed to contribute in each case mainly, in

Ẽx I π Ref. A corr 0 c 2 LJ lj a ptra c d,p LJ 2 G λ A 2 /A 0
Eq.( 4) Eq.( 5) Eq.(A.5) Eq.( 10) [START_REF] Valnion | [END_REF]21] Eq.( 2) comparison to the SSM (Table 1 in Paper P). We define a cross section corrected for the change in the penetrability of the dominant configuration LJ lj [Eq. (A.5)],

Ref

A corr 0 ( Ẽx , I π , LJ lj) = 1 N N i=1 A 0 ( Ẽx , I π , E i p , LJ) a ptra LJ lj E i p , Ẽx 2 a Lz (E i p , LJ) . (4)
Tables 1 and2 show penetrability-corrected cross sections for E x ≈ 5.9 and 5.0 MeV, respectively. They correspond to Table T where data for states in the wider range 4.6 < E x < 6.4 MeV is listed more completely: The excitation energies, the mean cross section A mean 0 [Eq. ( 3)] on all IARs, as well as spectroscopic factors for the 207 Pb(d, p) reaction [START_REF] Valnion | [END_REF] and the 209 Bi(d, 3 He) reaction [24] are given.

The s.p. widths Γ s.p. LJ for particles LJ and the s.p. widths Γ s.p. lj for holes lj are not well known (see Table 2 in Paper S). The agreement of some products Γ s.p. LJ Γ s.p. lj with earlier results is discussed in Sec. 4.1.2. However, in this paper we discuss only ratios of the cross section A corr 0 together with further simplifying assumptions (Secs. 4.4, 4.5, 4.6, and 4.7).

States with one major configuration.

A weak admixture of a secondary configuration LJ l j does not have much effect upon the mean cross section A calc 0 [Eq. (A.13)] provided that Γ s.p. l j Γ s.p. lj . In this paper we discuss states with major d 5/2 p 3/2 components. Here, weak admixtures of the configuration d 5/2 f 5/2 may be neglected since Γ s.p.

f 5/2 = 0.25 Γ s.p. p 3/2 [6,22,23]; admixtures of d 5/2 p 1/2 with Γ s.p. p 1/2 = 1.8 Γ s.p. p 3/2 are small [14,[START_REF] Valnion | [END_REF]21], see Table T. (For states with major d 5/2 p 1/2 components, weak admixtures of the configurations d 5/2 f 5/2 , d 5/2 p 3/2 may be similarly neglected.) The discussion in Sec. 4 is based on two assumptions.

(i) By assuming a single neutron particle-hole configuration LJ lj to contribute to the cross section of a certain state Ẽx I π , the calculated mean cross section [Eq. (A.13)] is described by a single term. Hence the penetrability-corrected cross section blue for the fit of the smooth background, green for each peak, red for the whole spectrum). By fitting the 6.01 MeV peak with a single level (left), the residuum spectra exceed the ±2 σ limits (dashed lines) and the excitation energies determined from the 208 Pb(p, p ) spectra on the two IAR are different in contrast to a fit of the 6.01 MeV peak with two levels (right). The peaks at Ex ≈ 5.986 and 6.024 MeV are satellites from neighbouring peaks because L-electrons carry away about 14 keV from the proton energy [START_REF] Heusler | In Changing Facets of Nuclear Structure[END_REF] (see Sec. 3).

[Eq. ( 4)] may be used to determine the d 5/2 p 3/2 strength in states at E x ≈ 5.9 MeV and their distribution across the four different spins 1 -to 4 -. (ii) In the discussion (Sec. 4), another important assumption is the completeness of the total strength for the configuration LJ lj in a certain group of states with spins |J -j| ≤ I ≤ J + j. Then, by introducing relative cross sections

σ rel ( Ẽx , I π , LJ lj) = A corr 0 ( Ẽx , I π , LJ lj) Ẽx I A corr 0 ( Ẽx , I π , LJ lj) , (5)
the complete LJ lj strength for all states Ẽx I π with spin I yields

Ẽx σ rel ( Ẽx , I π , LJ lj) = (2I + 1), (6)
meaning that the relative cross section σ rel [Eq. 5)] is proportional to the configuration strength |c LJ lj | 2 [Eq. ( 1)]. The total sum of all relative strengths yields

I Ẽx σ rel ( Ẽx , I π , LJ lj) = (2J + 1)(2j + 1), (7)
since J+j I=|J-j|

(2I + 1) = (2J + 1)(2j + 1). (8)
In the selected group of states at E x ≈ 5.9 MeV, the spin I of each state can be determined by observing the strength of each configuration LJ lj, and especially d 5/2 p 3/2 , to be in the range 0 < c 2 < 1 within the experimental uncertainties. Two states distinguished by different energy labels Ẽx may build an unresolved doublet with the same excitation energy E x . In such a doublet, the sum of the strengths c 2 weighted by the spin factors 2I + 1 can only be determined.

Determination of excitation energies

The line shape for 208 Pb(p, p ) is asymmetric. The instrumental resolution is 3 keV, the energy loss in the target produces a tail to each peak, see Paper S. In addition, satellite peaks from the 208 82 Pb(p, p + z e -) 208 82 Pb z+ reaction appear [START_REF] Heusler | In Changing Facets of Nuclear Structure[END_REF]. Some satellite peaks are comparable to peaks from physical states. Hence, determining excitation energies depends on several circumstances of the fit.

Determining energies using GASPAN.

Excitation energies are determined by using GASPAN [16] and additional computer programs described in Paper S. The calibration of the excitation energies has been verified as correct to within better than 1 keV for 3.1 < E x < 8.0 MeV by comparison with NDS2007.

Fig. 1 shows 208 Pb(p, p ) spectra of only 70 keV length taken on the d 5/2 and g 7/2 IARs for 3.1 < E x < 8.0 MeV on a logarithmic scale. The peak-to-valley ratio is 500 : 1. Besides the strong peaks at E x = 5.97 and 6.01 MeV, several minor peaks are fitted. Some peaks are satellites from the electron knockout reaction [START_REF] Heusler | In Changing Facets of Nuclear Structure[END_REF].

On the d 5/2 IAR, by assuming a single level for the 6.01 MeV peak, the residuum spectra deviates by up to ±4σ in a systematic manner (bottom left frame). Also, the excitation energies differ by 1.2 keV for the spectra taken on the d 5/2 and g 7/2 IARs (left frames). Clearly, the 6.01 MeV peak can be fitted only by assuming two levels having a distance of about 1.5 keV (right frames).

Mean excitation energy and variance.

The excitation energy

E GASP AN x ± δE GASP AN
x is determined from the fit by GASPAN for several runs of typically 30 minutes beam time. Typical individual uncertainties are between 0.04 keV for levels with high statistics and 0.30 keV for levels with low cross sections and low peak-tovalley ratios (e.g. 2 µb/sr and 3 : 1). Fig. 2 shows the values E GASP AN x and their uncertainties for the 5778 2 -state and the E x = 5.81 MeV doublet.

A study for the 6.01 MeV doublet reveals similar shifts of the values E GASP AN x with scattering angle and proton energy but with a spread of excitation energies within a range of about 2 keV. The fit by GASPAN is more difficult than for the 5.81 MeV doublet because of the presence of satellite peaks from L-electrons and several other close-lying physical states (Fig. 1). The mean value does not depend on the proton energy or scattering angle for the 5778 2 -state (Fig. 2, left frames). However for the E x = 5.81 MeV doublet, the mean value on the g 9/2 IAR is different from that on the d 5/2 IAR (right frames). On the g 9/2 IAR, the values are less precise because of the much lower mean cross section (15 µb/sr on the g 9/2 IAR vs. 270 µb/sr on the d 5/2 IAR, see Table T).

E GASP AN x = n i=1 [E GASP AN x (i)/δE GASP AN x (i)] 2 n i=1 [1/δE GASP AN x (i)] 2 (9)

Data from the 207 Pb(d, p) reaction

The strength of the configuration LJ p 1/2 is derived from the spectroscopic factors G λ LJ determined by the 207 Pb(d, p) reaction [START_REF] Valnion | [END_REF]21], In this sense, for spins 1 -and 2 -, the gaps between the respective configurations d 3/2 p 1/2 and s 1/2 p 3/2 with E SSM x = 6.033 and 6.371 MeV, and for spins 3 -and 4 -, the gaps between the respective configurations g 7/2 p 1/2 and j 15/2 i 13/2 with E SSM x = 5.922 and 6.487 MeV are considered to be sufficiently large, see Table 1 in Paper P. For unnatural parity states the mean matrix element of the residual interaction is about 50 keV, for natural parity about 150 keV [3]. Hence at least for unnatural parity states, the unitarity relations may be assumed to hold for the states at E x < ∼ 6.3 MeV. At 5.7 < E x < 6.3 MeV (Table 1 of Paper P) a strong mixing of the d 5/2 p 3/2 configuration with other configurations can be expected, since the SSM predicts three, four, four and three respective configurations for spins 1 -, 2 -, 3 -and 4 -within 300 keV.

c d,p LJ ( Ẽx , I π ) 2 = 2 2I + 1 G λ LJ ( Ẽx , I π ). ( 10
)
Fig. 4 in Ref. [22] shows an excitation function of the 6087 2 -state. It is strongly excited on the d 3/2 IAR and weakly on the s 1/2 IAR. The cross section at E p = 16.45 MeV is explained by the Lorentzian tails. A vanishing excitation on the d 5/2 IAR is assumed (Table 1). The 6191 3 -state [5] is weakly excited on all IARs (see Figs. 1 and2 in Ref. [13]). The 6087 2 -and 6191 3 -states are not discussed any further below, but shown in Fig. 3.

For states containing only d 5/2 p 1/2 components, the shape of the angular distribution on the d 5/2 IAR should be isotropic, K max = 0. The non-vanishing anisotropy coefficients A 2 /A 0 of the angular distribution indicate admixtures of d 5/2 f 5/2 and d 5/2 p 3/2 (Table 2). However, as is known [25], even small admixtures produce large, non-vanishing values A K /A 0 , K > 0 [Eq. (A.11)] through the interference with the ‡ NDS2007 adds the following footnote to the 5813.27 level observed by 208 Pb(n, n γ). 2005YaZW report a transition with Eγ =2338.765 14 with Iγ =1.30 6 placed from the 4953 level. In (d,pγ), there is a 2338γ placed from the 5813 level. From Iγ (2338γ)/Iγ (3199γ)=0.85 14, one expects Iγ =0.133 23 for a component of the 2338γ placed from the 5813 level in (n,n'γ). This leaves Iγ=1.17 6 for placement of the 2338γ from the 4953 level. Eγ for placement from the 5813 level is taken from the level energies.

§ NDS2007 adds the following footnote to the 6009.77 level observed by 208 Pb(n, n γ). 2005YaZW report levels at 6010, deexciting via Eγ =1924 and 2535, and at 6011 deexciting via Eγ =1687 and 3397.

For the 6011 level, the two transitions give discrepant energies, the 1687 γ giving E(level)=6010.72 15 and the 3397 γ giving 6011.64 6. For the 6010 level, the two transitions give consistent energies of 6009.78 11 and 6009.78 3, respectively. The 1687 γ thus gives an E(level) lying 1 keV below that given by the 3398 γ and 1 keV above that given by the 1924 and 2535 γ's. In (d,p γ), the 3397 γ is not reported and a 1686 γ is placed from the 6010 level. The branchings of the 1686, 1924, and 2535 γ's in (d,p γ) and (n,n' γ) agree well. Thus, in spite of the energy discrepancy, the (d,p γ) scheme seems correct and is adopted by the evaluator. The 1687 γ is thus included with the 6010 level and removed from the 6011 level. Eγ is taken from the level energies.

major configuration. Hence, a dominant d 5/2 p 1/2 strength may be still assumed in these states.

For states containing only d 5/2 p 3/2 components, the shape of the angular distribution on the d 5/2 IAR is described by Legendre polynomials with K max = 2, see Eq. (A.12). For a d 5/2 f 5/2 admixture, K max = 4 is needed to describe the shape [Eq. (A.11)]. Kulleck at al [17] found the anisotropy coefficients A K /A 0 for K = 4, 6 to be small for all levels with 4.6 < E x < 6.3 MeV which do not contain a natural parity state. Hence, a dominant d 5/2 p 3/2 strength may be still assumed in these states.

Table 3 shows the values for LJ = d 5/2 and lj = p 1/2 , p 3/2 . They differ from previous values by 10 -30% [6,26]. The exact value of Γ s.p. d 5/2 Γ s.p. p 3/2 is not important for the following discussion (Secs. 4.2, 4.4, 4.5, 4.6, and 4.7).

Center of gravity.

The center of gravity of states assumed to contain the complete LJ lj strength for a certain spin is calculated from the penetrability-corrected cross sections A corr 0 [Eq. ( 4)],

E Cof G x (LJ lj, I π ) = N i=1 Ẽi x A corr 0 ( Ẽi x , I π , LJ lj) N i=1 A corr 0 ( Ẽi x , I π , LJ lj) . ( 12
)
The center of gravity for all states assumed to contain the complete LJ lj strength is given by

E Cof G x (LJ lj) = I (2I + 1)E Cof G x (LJ lj, I π ) I (2I + 1) . ( 13
)
For the states with major d 5/2 p 3/2 fractions, the common center of gravity agrees well with the SSM prediction (Table 4).

The difference between the center of gravity for all states assumed to contain the complete LJ lj strength [Eq. (13] and that for a certain spin [Eq. (12] defines the multiplet splitting. For the d 5/2 p 1/2 and d 5/2 p 3/2 configurations, Table 4 shows the center of gravity for each spin E Cof G x [Eq. ( 12)], the center of gravity for the complete configuration E Cof G x [Eq. ( 13)], and the multiplet splitting

E Cof G x -E Cof G x
. Fig. 3 displays the values for the d 5/2 p 3/2 configuration.

In a similar way to the states with dominant d 5/2 p 3/2 strength, the center of gravity for all states with major d 5/2 p 1/2 fractions (Sec. 4.1.2) is in line with the SSM prediction (Table 4). The multiplet splitting is much smaller than for d 5/2 p 3/2 .

LJ lj

Reaction E SSM

x

I π E Cof G x E Cof G x E Cof G x -E Cof G x [2]
[Eq. ( 13)] [Eq. ( 12 From the SSM, most of the d 5/2 p 3/2 strength is expected in the region E x ≈ 5.9 MeV. By assuming the spin assignments of NDS2007 and by assuming that the complete d 5/2 p 3/2 strength [Eq. ( 7)] has been identified, no apparent combination of the measured cross sections A corr 0 can reproduce the expected values (2I + 1) [Eq. ( 6)]. With any assumed spin for the 6012 state, the sum of the d 5/2 p 3/2 strength largely differs from the expected values; even the sum for some combinations of spins disagrees by factors two to five, see Table 5 under assumptions A-D.

In Secs. 4.3 and 4.7, the 5886 state is assigned the spin of 4 -. In Sec. 4.4 the 6012 state is shown to have the spin of 4 -. In Secs. 4.5 and 4.6.1, the missing 1 -and 2 - strengths are located. In Secs. 4.6.1, the 5.81 MeV level is shown to be a doublet of two states with spins 2 -and 3 -. The energy spacing is less than 1 keV (Sec. 4.6.3).

Tentatively, Kulleck at al [17] already assigned a spin of 4 -to the 6.01 MeV level. Yet the doublet was not resolved, and the even more strongly excited 5.81 MeV level was not analyzed because of contaminations from the 208 Pb(p, d) reaction.

Figure 3:

Excitation energy for states with major fragments of the configuration ν(d 5/2 +1 p 3/2 -1 )I π with spins I π = 1 -, 2 -, 3 -, and 4 -(Table 1). The thin bar shows the full strength c 2 = 1 while the thick bar shows the d 5/2 p 3/2 strength determined by Eq. ( 5). The global center of gravity [dashed line continued across all spins, Eq. ( 13)] is in line with with the SSM value shown on the left. The center of gravity for spins 1 -and 4 -is higher than the global center of gravity, for spins 2 -and 3 -it is lower (diamonds on the right).

The 5886 state

According to several arguments, the 5886 state is assigned a spin of 4 -; only three are mentioned here. In Sec. 4.4 an additional argument is given. (In Fig. 1 of Paper P, the 5886 state is prematurely already assigned a spin of 4 -.)

Among over a hundred identified states below E x = 6.1 MeV, the lowest distance of any pair of states with equal spin and parity is determined as being larger than 25 keV [2,6,12,13,14]. Therefore, the 5874 and 5886 states have different spins.

The shape of an angular distribution for the pure d 5/2 p 3/2 states with spins 3 -and 4 -differs greatly, A K /A 0 = -0.629 and +0.289, respectively [Eq. (A.12)]. If they are assigned spins 3 -and 4 -, respectively, the angular distributions of the 5874 and 5886 states have a shape corresponding to the configuration d 5/2 p 3/2 , but not in the case of the reverse assignment.

For the 5874 state, the angular distribution of the 207 Pb(d, p) reaction is well described by LJ = d 5/2 [21] whereas for the 5886 state two components are needed. Valnion [21] proposed a mixture of LJ = d 5/2 and g 7/2 with J = L + 1 2 and L-1 2 , respectively, to explain the angular distribution for the polarization asymmetry. However, an equally good fit is obtained with LJ = g 7/2 and g 9/2 and spectroscopic factors G λ g 7/2 ≈ 0.10, G λ g 9/2 ≈ 0.03 [Eq. [START_REF] Heusler | Maier-Leibnitz Laboratorium[END_REF]] signifying a spin of 4 -.

, I π B , • • •) is evaluated below the respective heading A, B, • • •.
If the spins are correctly assigned (as finally achieved by assumption H), the sum should correspond to 2I + 1 [Eq. ( 6)]. The question mark hints at assumed assignments, the exclamation mark new assignments. Assumptions are ruled out if the sum rule for any spin or any pair of spins deviates significantly from the expected value of 2I + 1. Values which are too low are printed in italics, values which are too high in bold face.

Ẽx σ

rel I π A I π B I π C I π D I π E σ rel I π F I π G σ rel I π H Eq. (5) a Eq. (5) b Eq. (5) c I π 2I +1
Configuration strength 11.9 (c) Relative cross sections of levels investigated by assumption H.

P I P Ẽx σ rel [Eq. (6)] for assumption A-H A B C D E F G H (accepted) 1 - 3 7

Resolving the 6.01 MeV doublet

The 6012 state is essentially excited only on the d 5/2 IAR while the 6010 state is excited on all IAR (at proton energies 14.8 < E p < 18.2 MeV). The 6.01 MeV doublet is resolved by the 208 Pb(p, p ) reaction (Sec. 3.3.1 and Fig. 1) yielding E GASP AN x = 6009.6±0.2 and 6011.4±0.3keV from about 60 runs (Table T). The 207 Pb(d, p) reaction excites only the 6010 3 -state yielding E GASP AN x = 6009.5 ± 0.1 keV [Table T, Eq. ( 9)]. The 6010 3 -state contains 29% g 7/2 p 1/2 strength [START_REF] Valnion | [END_REF]. The difference between the excitation energies for the 6010 and 6012 states is determined as 1.8 ± 0.3 keV. NDS2007 declares a distance of 1.9 ± 0.1 keV.

On the d 5/2 IAR for scattering angles of Θ > 100 • , the 6010 3 -state is much more weakly excited than the 6012 state, on all other IAR much more strongly. The mean cross section A mean 0 [Eq. ( 3)] of the 6.01 MeV doublet agrees with the values from Ref. [17] (Table T).

The 6012 state is assigned the spin of 4 -by observing the mean cross sections A corr 0 [Table 1, Eq. ( 4)]. Namely, the large cross section for the 6012 state excludes spins of 1 -, 2 -, and 3 -(see assumptions A-C in Table 5).

Assumption D excludes spin 3 -for the 5886 state and the spin of 4 -is assigned by assumption E. The assignment deduced in Sec. 4.3 is thus confirmed. The 5886 4 -state contains the remaining d 5/2 p 3/2 strength, while the 5969 4 -state is only weakly excited on the d 5/2 IAR (Table 1). The sum of the d 5/2 p 3/2 strength in the 5886 4 -, 5969 4 -, and 6012 4 -states agrees with the expectation within 15% while for spins 3 -, 2 -, and 1 -there are still discrepancies of 50% and factors two and five, respectively.

The shape of the angular distribution of the 6012 4 -state is well described by assuming a pure d 5/2 p 3/2 component (Eq. (A.12) and Table 1). From experiments with semiconductor detectors [17], the minimum at Θ ≈ 90 • predicted for the d 5/2 p 3/2 angular distribution is verified by our data for the 6.01 MeV doublet. As suggested by Kulleck at al [17], the direct (p, p ) contribution to the 6010 3 -state (starting already at Θ ≈ 120 • ) partially explains the expressive A 4 /A 0 contribution near Θ = 90 • . Yet in addition, a weak d 5/2 f 5/2 admixture to the 6012 state may also explain the large value A 4 /A 0 [Eq. (A.11)].

4.5. The d 5/2 p 3/2 1 -states Only one 1 -state is located within 5.7 < E x < 6.1 MeV; much d 5/2 p 3/2 strength for the spin of 1 -is missing. Indeed, the 6264 1 -and 6314 1 -states contain the missing fractions (Table 1), but the penetrability is low, namely a ptra d 5/2 p 3/2 = 0.60 and 0.32, respectively. Fig. 2 of Ref. [13]) shows spectra with the 6264 1 -and 6314 1 -states.

An open question is the spin and structure assignment to the five states at 5.638 < E x < 5.650 MeV where the highest member is a (9) + state (Paper S); one state may have a spin of 1 - [START_REF] Heusler | [END_REF]. Apparently, the configuration mixing among the 1 - states is large. Therefore the uncertainty of the center of gravity [Eq. (13)] may be even greater than shown in Fig. 3 and Table 4. A large fraction of the configuration g 7/2 p 1/2 is derived from the 207 Pb(d, p) reaction (Table T). The cross sections A corr 0 of the 5.81 MeV level and the 6012 4 -state are similar (Table 1). Since the 5813, 5874 and 6010 states contain some fraction of the d 5/2 p 3/2 3 -strength and the 5886 and 6012 states most of the d 5/2 p 3/2 4 -strength, another state with a spin of either 1 -or 2 -must be present. We are assigning the energy labels 5812 and 5813 to the members of the doublet.

Spins of 1 -and 4 -are excluded since the sum of the cross sections for all states with spins of 1 -and 4 -in the region 5.7 < E x < 6.4 MeV would exceed expectations (see assumptions F and G in Table 5). Hence the 5812 state has a spin of 2 -.

The distribution of the d 5/2 p 3/2 strength across the four spins (1 -, 2 -, 3 -, and 4 -) is in full agreement with this finding. The cross section of the 5.81MeV doublet can be divided up between the 5812 2 -state and the 5813 3 -state by about 45 ± 5 : 55 ∓ 5%, see assumption H in Table 5.

4.6.2.

Particle-hole configurations in the 5812 and 5813 states. The 5.81 MeV doublet consists of two states with spins of 2 -and 3 -(Sec. 4.6.1). In the Q3D experiments, the doublet at E x = 5.81 MeV is not resolved, see however Secs. 3.3.2 and 4.6.3. By assuming the cross section to be shared equally by both the 5812 2 -state and the 5813 3 -state, the angular distribution of the 5.81 MeV doublet is well described. A superposition of d 5/2 p 3/2 components with spins of 2 -and 3 -yields an angular distribution with a high maximum near Θ = 90 • . It is derived from the 3 -state, while the angular distribution for the 2 -state is rather flat [Eq. (A.12)].

A 12% admixture of the configuration g 7/2 p 1/2 to the 3 -state in the 5.81 MeV level is determined from the 207 Pb(d, p) reaction [START_REF] Valnion | [END_REF]; a d 3/2 p 1/2 admixture to the 2 - state is estimated to be much less.

4.6.3.

Excitation energies of the 5812 and 5813 states.

The excitation energy derived from 207 Pb(d, p) is higher than the mean value derived from the 208 Pb(p, p ) reaction (Table T). Since only the 3 -state is excited by the 207 Pb(d, p) reaction (LJ = g 7/2 [START_REF] Valnion | [END_REF]21]), the 5813 3 -state has a higher excitation energy than the 5812 2 -state.

From the SSM, another major component is assumed to be g 9/2 f 7/2 ; the weak cross section on the g 9/2 IAR is in line with such an interpretation (Table T). The mean excitation energy determined on the g 9/2 IAR is lower than on all other IARs (Fig. 2),

E GASP AN x (E p = 14.92 MeV) = 5812.96 ± 0.27 keV, E GASP AN x (15.7 < E p < 18.1 MeV) = 5813.23 ± 0.17 keV. ( 14
)
The difference is explained by the different angular distributions of the g 9/2 f 7/2 and d 5/2 p 3/2 components on the respective IARs. Similar as for the 5686, 5694, 5836 states with spins of 6 -, 7 -, and 8 - [12], respectively, the 5812, 5813 may be described by pure g 9/2 f 7/2 configurations. The excitation of the 5812 2 -is enhanced on the g 9/2 IAR at scattering angles |π -Θ| ≈ 90 • and, hence, the mean excitation energy is lower (large symbols in Fig. 2). As a common result from the 208 Pb(p, p ) and the 207 Pb(d, p) reactions, we obtain The sum of the d 5/2 p 3/2 strength [Eq. ( 7)] located in twelve states at 5.7 < E x < 6.4 MeV is complete within about 10% (assumption H in Table 5).

E x (5812 2 -) = 5812.8 ± 0.2 keV, E x (5813 3 -) = 5813.2 ± 0.2 keV, E x (5813 3 -) -E x (5812 2 -) = 0.4 ± 0.2 keV. ( 15
For spin 1 -a strength of 2.8±0.7 is determined while 2I +1 = 3 [Eq. ( 6)] is expected, for spin 2 -a strength of 4.8 ± 0.5 ± 0.3 compared to 2I + 1 = 5, for spin 3 -a strength of 6.9 ± 0.7 ∓ 0.3 compared to 2I + 1 = 7, for spin 4 -a strength of 9.4 ± 1.5 compared to 2I + 1 = 9. Fig. 3 displays the distribution of the d 5/2 p 3/2 strength. The 5812 2 -and 6012 4 -states contain about 60% and 70% of the d 5/2 p 3/2 strength, respectively. Each of the three 3 -states consists of at least three configurations having less than 60% strength. While the 5947 1 -state contains 90% of the d 3/2 p 1/2 strength [START_REF] Valnion | [END_REF] and the 6314 1 -state 60 -80% of the s 1/2 p 3/2 strength [13], the 6264 1 -state has a complex structure.

The finding of the d 5/2 p 3/2 strength to be complete within 10% verifies the assumption made in Sec. 4.1.1 that only minor d 5/2 p 3/2 fractions contribute to the states at E x < 5.7 MeV.

Summary

In the doubly magic nucleus 208 Pb, the SSM predicts eighty negative parity particle-hole states for E x < 6.1 MeV. Only the structure of all negative parity states at E x < ∼ 4.8 MeV and a few more higher-lying states is known up to now. In order to determine the residual interaction among particle-hole configurations, the structure of all eighty states must be known, however.

Twelve states with spins of 1 -, 2 -, 3 -, and 4 -are studied in detail at excitation energies 5.75 < E x < 6.35 MeV. They are identified as containing the complete strength of the neutron particle-hole configuration d 5/2 p 3/2 within at least 90% for each spin. The center of gravity for the configuration d 5/2 p 1/2 and d 5/2 p 3/2 is in line with the predictions by the shell model without residual interaction (SSM) within 10 keV. The centroid energies of the d 5/2 p 3/2 configuration for spins of 1 -to 4 -follow a parabolic shape with deviations from the global centroid between -90 keV (2 -) and +280 keV (1 -).

The doublet at E x = 6.01 MeV is resolved into two states having a distance of 1.8 keV, E x = 6010 and 6012 keV. In line with the tentative assignment by Kulleck et al [17], the spin and parity of the 6012 state is assigned as 4 -; the excitation energy tentatively determined by NDS2007 is confirmed. The 5886 state has been newly assigned the spin of 4 -. The d 5/2 p 3/2 4 -strength is almost completely located in the 5886 and 6012 states. The 6010 3 -state has a complex structure similar to the 5813 3 -state.

The doublet at E x = 5.81 MeV is resolved into two states having a distance of 0.4 keV, E x = 5812 and 5813 keV. The 5812 state has been newly identified; it has the spin of 2 -and consists of the two configurations d 5/2 p 3/2 and g 9/2 f 7/2 with a strength of about 50% each. The 5813 3 -state contains about 10% g 7/2 p 1/2 strength besides about half of the d 5/2 p 3/2 strength.

The mean cross section of many states with major d 5/2 f 5/2 strength in the region 5.3 < E x < 5.7 MeV is greatly increased by weak d 5/2 p 1/2 , d 5/2 p 3/2 admixtures. However, the knowledge of the distribution of almost the complete d 5/2 p 1/2 and d 5/2 p 3/2 strengths presented in this paper restricts such admixtures to be less than 10%. By this means, the structure of most negative parity particle-hole states up to E x ≈ 6.1 MeV can be determined in future. States with dominant proton particle-hole configurations may also be studied in more detail. Ultimately, the residual interaction among all eighty particle-hole configurations with negative parity at E x < 6.1 MeV in 208 Pb can be determined with high reliability. may be viewed as a two-stage process 1 → r → 2. Starting from the ground state E 0

x 0 + in 208 Pb, in the initial stage a LJ IAR is created. In the subsequent proton decay of the IAR, neutron particle-hole states in 208 Pb with excitation energy E I x , and final spin and parity I π are excited and a proton lj is emitted,

E 0 x 0 + ⊗ LJ → E p LJ → E I x I π ⊗ lj. (A.1)
In the E I x I π state, the particle LJ and the hole lj couple to the final spin I. The decay amplitude γ(E p LJ → E I

x I π ⊗ lj) is normalized in such a way that the absolute square is the partial width (see Eq. (3F-1) in Ref. [START_REF] Bohr | Nuclear Structure[END_REF])

Γ(E p LJ → E I x I π ⊗ lj) = |γ(E p LJ → E I x I π ⊗ lj)| 2 . (A.2)
This gives the probability per unit time (multiplied by ) for the decay into the channel considered. The time reversal invariance implies Γ(E p LJ → E I x I π ⊗ lj) = Γ(E I

x I π ⊗ lj → E p LJ). Since the decay involves the transmission through a large barrier associated with centrifugal and Coulomb forces, the decay amplitudes strongly depend on the energy of the proton.

In the 208 Pb(p, p ) reaction via IAR, the protons enter with the beam energy E p and leave the two-stage process with the energy given by the difference between the beam energy E p and the excitation energy of the E I

x I π state, E p (E I

x , E p ) = E p -E I x .

(A.3)

In the SSM, the energy for populating a model configuration does not depend on the spin I,

E SSM p (LJ lj) = E p -E SSM x (LJ lj). (A.4)
To state the case briefly, the decay width Γ(E p LJ → E I x I π ⊗ lj) [Eq. (A.

2)] is written as Γ lj (E p ) where E p is the energy of the scattered proton [Eq. (A.3)]. The decay width Γ lj depends strongly on the proton energy [START_REF] Clarkson R G, Brentano P Von | [END_REF], in fact on the penetrability of the outgoing proton through the Coulomb barrier and the additional centrifugal barrier. We define the relative penetrability as depends only slightly on the orbital momentum l. The ratio of the penetrabilities for the spin j = l + 1 2 to j = l -1 2 increases, the higher the angular momentum l becomes. Typically, the relative penetrability amplitude a ptra LJ lj changes by a factor of two for a 1 MeV difference in the proton energy [START_REF] Clarkson R G, Brentano P Von | [END_REF].

The decay width Γ LJ for the entrance channel is calculated to vary little in the vicinity of the IAR, provided that the beam energy does not change by more than a few percent [START_REF] Clarkson R G, Brentano P Von | [END_REF]. This holds for the data considered in this work.

For convenience, single particle (s.p.) widths are introduced by defining Γ

( a )

 a Instead of the uncertainty of the cross sections A corr 0 uncertainty of the strength c 2 is shown. (b) Secs. 4.6.1, 4.6.3. (c) For the unresolved 5812 + 5813 doublet. (d) See assumption H in Table 5 discussed in Sec. 4.6.1. (e) Sec. 4.3. (f ) Sec. 4.4. (g) Included for completeness (Sec. 4.1.1).

Figure 2 :

 2 Figure 2: Excitation energies determined by using GASPAN (Sec. 3.3.1), (left) for the 5778 level, (right) for 5813 level. For each run, E GASP AN x is shown by a dot, the uncertainty δE GASP AN x by a horizontal bar. The runs are sorted according (i) to the scattering angle, (ii) to the proton energy. The vertical lines show the mean value E GASP AN x [Eq. (9)] and the uncertainty determined from the values E GASP AN x for all runs, (top left) 5777.86 ± 0.20, (bottom left) 5777.92±0.15 keV, (top right) 5813.23±0.17, (bottom right) 5812.96±0.27 keV. On the g 9/2 , the excitation energies for 80 • ≤ |180 • -Θ| ≤ 90 • (large symbols) are lower than for Θ < 80 • .

  (a) Relative cross sections of levels investigated by assumptions A-E. (b) Relative cross sections of levels investigated by assumptions F-G.

4. 6 .

 6 Resolving the 5.81 MeV doublet 4.6.1. Spin of the 5812 state. The 5.81 MeV level certainly contains a 3 -state [5].

) 4 . 7 .

 47 The complete d 5/2 p 3/2 strength

2 =

 2 a ptra LJ lj E p , E I x Γ lj E p (E I x , E p ) Γ lj E SSM p (LJ lj) = Γ lj E p -E I x Γ lj [E p -E SSM x (LJlj)] . (A.5) This becomes unity on top of the LJ IAR (E p = E res LJ ) at the excitation energy E SSM x (LJ lj) of the SSM configuration [Eq. (A.4)]. The relative penetrability a ptra LJ lj

  are determined from data taken on the d 5/2 IAR shown in TableT. Here the factor a Lz varies by 0.6 -1 for proton energies 16.4 < E p < 16.6 MeV. The anisotropy coefficient A K /A 0 [Eq. (2)] is given for K = 2.

	3.2. Determination of s.p. widths, spins, configuration strengths, and excitation
	energies
	3.2.1.	Basic assumptions.	The basic assumption of an overwhelming resonant
	reaction is well founded for states with unnatural parity at scattering angles Θ > ∼ 40 • and
	for states with natural parity at scattering angles Θ > ∼ 100 • [2, 6, 12, 13, 14, 17, 22, 23]

Table 1 :

 1 For states in208 Pb with excitation energies 5.7 < Ex < 6.4 MeV, strongly excited by the 208 Pb(p, p ) reaction via the d 5/2 IAR with A mean 0 > 10 µb/sr, the energy label, the spin I π , the cross sections A mean

				0	[Eq. (3)], A corr 0	, the fraction c 2 of the d 5/2 p 3/2 strength,
	the calculated penetrability, and the anisotropy AK/A0 are shown. The state 5812 has been
	newly identified, the spins for the 5886 and 6012 states are new (printed in bold face).
					LJ lj = d 5/2 p 3/2	on d 5/2 IAR
	Ẽx	I π		Ref.	A corr 0	c 2 LJ lj	a ptra	A mean 0	A 2 /A 0
		[5]	this		Eq.(4)	Eq.(5)	Eq.(A.5)	[14]	Eq.(2)
			work		µb/sr	×100 a		µb/sr
	5778 2 -, 3 -	2 -	[13]	68 33± 2		0 . 9 5 3	7 5 + 0 .2± 0.1
	5812		2 -b	c	54± 7±6 d	0.973	c	c
	5813	3 -	3 -	[5]	256 c 50± 5∓5 d	0.973	270 c -0.7± 0.2
	5874	3 -	3 -	[5]	96 33± 5		1 . 0 0 8	9 5 -0.2± 0.3
	5886	3 -	4 -e	108 29± 5		1 . 0 1 5	1 0 5 + 0 .3± 0.2
	5924	2 -	2 -	[5]	16	7± 3		1 . 0 2 0	1 5 -0.5± 0.3
	5947	1 -	1 -	[5]	21 17± 4		1 . 0 2 7	2 0 + 0 .7± 0.3
	5969	4 -	4 -	[5]	26	7± 2		1 . 0 3 4	2 5 + 0 .3± 0.2
	6010	3 -	3 -	[5]	48 16± 8		1 . 0 9 5	4 0 -0.0± 0.5
	6012		4 -f	252 68±16	1.096	210 +0.3± 0.2
	6087	1 -	2-[13]	< 13	< 6		1 . 1 4 8	< 10 g
	6191	3 -	3 -	[5]	< 8	< 3		1 . 2 2 8	< 5 g
	6264	1 -	1 -	[5]	33 26±13	1.286	20 +0.9± 0.3
	6314	1 -	1 -	[5]	63 51±25	1.773	20 -0.0± 0.3

Table 2 :

 2 Data similar to Table 1 for excitation energies 4.2 < Ex < 5.3 MeV. In addition, the strength c 2 d 5/2 p 1/2 determined from the 207 Pb(d, p) reaction [with an estimated uncertainty of 10%, Eq. (10)] is shown. The non-vanishing value of A2/A0 indicates admixtures of d 5/2 f 5/2 and d 5/2 p 3/2 .

	LJ lj = d 5/2 p 1/2	on d 5/2 IAR
	208 Pb(p, p )

207

Pb(d, p)

Table 2

 2 NDS2007 discusses the 5.81 and 6.01 MeV doublets. The 5.81 MeV level is attributed to a single state ‡. The 6.01 MeV doublet is recognized as possibly consisting of two states §.4.1.1. Assumption of completeness.The spin assignments are mainly based on the assumption that the configuration strengths for configurations in a certain subspace is complete. Then, as shown for the lowest negative parity states with E x < 4.7 MeV[3], a sufficiently large gap in the model configuration space may yield unitarity relations for a group of states and configurations.

	shows values G λ LJ and c d,p LJ ( Ẽx , I π )	2	, see also Table T.

Table 3 :

 3 Product of s.p. widths Γ s.p. The SSM predicts four states at E x = 5.896 MeV with the d 5/2 p 3/2 configuration and two states at E x = 4.998 MeV with the d 5/2 p 1/2 configuration. Fourteen states at E x ≈ 5.9 MeV may be assumed to contain most of the d 5/2 p 3/2 strength (Table1). Six states in the region 4.2 < E x < 5.3 MeV contain almost the complete d 5/2 p 1/2 strength (Table2); both the 208 Pb(p, p ) reaction and the 207 Pb(d, p) reaction yield similar results for the distribution of the d 5/2 p 1/2 strength.With the assumptions (i) and (ii) declared in Sec. 3.2.3, the product of the decay widths can be determined from the sum of the penetrability-corrected cross sections

	LJ Γ s.p. lj .

π , LJ lj).

Table 4 :

 4 Center of gravity for states in208 Pb with major parts of the particle-hole configu- rations d 5/2 p 1/2 and d 5/2 p 3/2 .

Table 5 :

 5 Discussion of spin assignments to states in 208 Pb with excitation energies 5.7 < Ex < 6.4 MeV (Secs. 4.2, 4.4, 4.6.1, and 4.7). For each assumption A, B, • • •, the sum of the relative cross sections σ rel for different spin assignments (I π A

  , p 3/2 , f 5/2 , f 7/2 , h 9/2 , i 13/2 . (A.7)

	s.p. LJ = Γ[LJ → E 0 x 0 + ]	a t t h e r e s o n a n c e e n e r g y E res LJ ,	(A.6)
	Γ s.p. lj = Γ[g 9/2 → E SSM x	(g 9/2 lj)],	lj= p 1/2

Acknowledgments

We thank T. von Egidy, G. Graw, H.-L. Harney, J. Jolie, N. Pietralla, F. Riess, C. Scholl, V. Werner for discussions. This work has been supported by DFG Br799/20-1.

In this paper we mainly discuss angle integrated (mean) cross sections. For details, we refer the reader to more complete descriptions of the reaction mechanism in Sec. III of Ref. [23], Sec. 2.2 of Ref. [1], and Appendix 3F of Ref. [START_REF] Bohr | Nuclear Structure[END_REF].

The differential cross section is given by [1] dσ dΩ

The normalization constant is given by

2 .

(A.9)

Near the resonance, the shape of the excitation function is Lorentzian,

In general however, the energy dependence of the total width Γ tot LJ produces an asymmetry of the excitation function.

The geometrical anisotropy coefficients a K /a 0 defines the shape of the angular distribution. The shape is described by Legendre polynomials of maximum order K max ≤ min(2L, 2J, max(l, l ), max(j, j )).

(