Geneviève Dusson
email: genevieve.dusson@math.cnrs.fr.

Israel Michael Sigal
email: im.sigal@utoronto.ca.

Benjamin Stamm

The Feshbach-Schur map and perturbation theory

Keywords: Mathematics Subject Classification 2020. Primary 47A55, 35P15, 81Q15; Secondary 47A75, 35J10 Keywords. Perturbation theory, Spectrum, Feshbach-Schur map, Schroedinger operator, Atomic systems, Helium-type ions, Ground state

This paper deals with perturbation theory for discrete spectra of linear operators. To simplify exposition we consider here self-adjoint operators. This theory is based on the Feshbach-Schur map and it has advantages with respect to the standard perturbation theory in three aspects: (a) it readily produces rigorous estimates on eigenvalues and eigenfunctions with explicit constants; (b) it is compact and elementary (it uses properties of norms and the fundamental theorem of algebra about solutions of polynomial equations); and (c) it is based on a self-contained formulation of a fixed point problem for the eigenvalues and eigenfunctions, allowing for easy iterations. We apply our abstract results to obtain rigorous bounds on the ground states of Helium-type ions.

a robust mathematical field in the works of Kato and Rellich. It was extended to quantum resonances by Simon, see [START_REF] Cycon | Schrödinger Operators (with Applications to Quantum Mechanics and Global Geometry[END_REF][START_REF] Hislop | Introduction to Spectral Theory[END_REF][START_REF] Kato | Perturbation theory for linear operators[END_REF][START_REF] Reed | Fourier Analysis and Self-Adjointness[END_REF][START_REF] Reed | Methods of Modern Mathematical Physics IV: Analysis of Operators[END_REF][START_REF] Rellich | Perturbation Theory of Eigenvalue Problems[END_REF][START_REF] Simon | A Comprehensive Course in Analysis, Part 4: Operator Theory[END_REF][START_REF] Simon | Tosio Kato's work on non-relativistic quantum mechanics: Part 1 and 2[END_REF] for books and a book-size review.

A different approach to the eigenvalue perturbation problem going back to works of Feshbach and based on the Feshbach-Schur map was introduced in [START_REF] Bach | Quantum electrodynamics of confined non-relativistic particles[END_REF] and extended in [START_REF] Gustafson | Mathematical Concepts of Quantum Mechanics[END_REF][START_REF] Dusson | Analysis of the Feshbach-Schur method for the planewave discretizations of Schrödinger operators[END_REF].

In this paper, we develop further this approach proposing a self-contained theory in a form of a fixed point problem for the eigenvalues and eigenfunctions. It is more compact and direct than the traditional one and, as we show elsewhere, extends to the nonlinear eigenvalue problem.

We show that this approach leads naturally to bounds on the eigenvalues and eigenfunctions with explicit constants, which we use in an estimation of the ground state energies of the Helium-type ions.

The approach can handle tougher perturbations, non-isolated eigenvalues (see [START_REF] Bach | Quantum electrodynamics of confined non-relativistic particles[END_REF][START_REF] Gustafson | Mathematical Concepts of Quantum Mechanics[END_REF]) and continuous spectra as well as discrete ones. In this paper, we restrict ourselves to the latter. Namely, we address the eigenvalue perturbation problem for operators on a Hilbert space of the form

𝐻 = 𝐻 0 + 𝑊, (1.1)
where 𝐻 0 is an operator with some isolated eigenvalues and 𝑊 is an operator, small relative to 𝐻 0 in an appropriate norm. The goal is to show that 𝐻 has eigenvalues near those of 𝐻 0 and estimate these eigenvalues. Specifically, with • standing for the vector and operator norms in the underlying Hilbert space, we assume that (A) 𝐻 0 is a self-adjoint, non-negative operator (𝐻 0 ≥ 𝛽 > 0); (B) 𝑊 is symmetric and form-bounded with respect to 𝐻 0 , in the sense that

𝐻 -1 2 0 𝑊 𝐻 -1 2 0 < ∞.
(C) 𝐻 0 has an isolated eigenvalue 𝜆 0 > 0 of a finite multiplicity, 𝑚.

Here • is the operator norm and 𝐻 -𝑠 0 , 0 < 𝑠 < 1, is defined either by the spectral theory or by the explicit formula

𝐻 -𝑠 0 := 𝑐 ∫ ∞ 0 (𝐻 0 + 𝜔) -1 𝑑𝜔 𝜔 𝑠 ,
where 𝑐 :=

∫ ∞ 0 (1 + 𝜔) -1 𝑑 𝜔 𝜔 𝑠 -1
. It turns out to be useful in the proofs below to use the following form-norm

𝑊 𝐻 0 := 𝐻 -1 2 0 𝑊 𝐻 -1 2 0 ,
Let 𝑃 be the orthogonal projection onto the span of the eigenfunctions of 𝐻 0 corresponding to the eigenvalue 𝜆 0 and let 𝑃 ⊥ := 1 -𝑃. Let 𝛾 0 be the distance of 𝜆 0 to the rest of the spectrum of 𝐻 0 , and 𝜆 * := 𝜆 0 + 𝛾 0 . In what follows, we often deal with the expression

Φ(𝑊) := 𝜆 0 𝜆 * 𝛾 0 𝑃 ⊥ 𝑊 𝑃 2 𝐻 0 . (1.
2)

The following theorem proven in Section 2 is the main result of this paper.

Main Theorem 1.1. Let Assumptions (A)-(C) be satisfied and assume that, for some 0 <𝑏 < 1 and 0 <𝑎 < 1 -𝑏,

𝑃 ⊥ 𝑊 𝑃 ⊥ 𝐻 0 ≤ 𝑏𝛾 0 𝜆 * , (1.3)
𝑃𝑊 𝑃 + 𝑘Φ(𝑊) < 𝑎𝛾 0 , (1.4)

𝑘Φ(𝑊) < 1 2 (𝑎𝛾 0 -𝑃𝑊 𝑃 ), (1.5)
where 𝑘 := 1 1-𝑎-𝑏 . Then the spectrum of the operator 𝐻 near 𝜆 0 consists of isolated eigenvalues 𝜆 𝑖 of the total multiplicity 𝑚 satisfying, together with their normalized eigenfunctions 𝜑 𝑖 , the following estimates

|𝜆 𝑖 -𝜆 0 | ≤ 𝑃𝑊 𝑃 + 𝑘Φ(𝑊), (1.6)
𝜑 𝑖 -𝜑 0𝑖 ≤ 𝑘 √ Φ(𝑊) 𝛾 0 , (1.7)
where 𝜑 0𝑖 are appropriate eigenfunctions of 𝐻 0 corresponding to the eigenvalue 𝜆 0 .

Remark 1 (Comparison with [START_REF] Dusson | Analysis of the Feshbach-Schur method for the planewave discretizations of Schrödinger operators[END_REF]). A similar result was already proven in [START_REF] Dusson | Analysis of the Feshbach-Schur method for the planewave discretizations of Schrödinger operators[END_REF]. Here, the theory is made self-contained and formulated as the fixed point problem and the bounds are tightened.

Remark 2 (Conditions on 𝑊)

The fine tuning the conditions (1.3)-(1.5) on 𝑊 is used in application of this theorem to atomic systems in Section 3. Note that because of the elementary estimate

Φ(𝑊) ≤ 1 𝛾 0 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 ,
see (3.11) below, the computation of 𝑃𝑊 𝑃 and Φ(𝑊) reduces to computing the largest eigenvalues of the simple 𝑚 × 𝑚-matrices 𝑃𝑊 𝑃, -𝑃𝑊 𝑃 and 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃.

Remark 3 (Higher-order estimates and non-degenerate 𝜆 0 ). In fact, one can estimate 𝜆 𝑖 (as well as the eigenfunctions 𝜑 𝑖 ) to an arbitrary order in 𝑃 ⊥ 𝑊 𝑃 𝐻 0 /𝛾 0 . As a demonstration, we derive (after the proof of Theorem 1.1) the second-order estimate of the eigenvalue in the rank-one (𝑚 = 1) case

|𝜆 1 -𝜆 0 -𝑊 | ≤ 𝑘Φ(𝑊), (1.8)
where 𝜑 0 is the eigenfunction of 𝐻 0 corresponding to 𝜆 0 , 𝑊 := 𝜑 0 , 𝑊 𝜑 0 .

For degenerate 𝜆 0 , we would like to prove a similar bound on |𝜆 𝑖 -𝜆 0 -𝜇 𝑖 |, where 𝜇 𝑖 is the corresponding eigenvalue of the 𝑚 × 𝑚-matrix 𝑃𝑊 𝑃. Here, we have a partial result (proven at the end of the next section) for the lowest eigenvalue 𝜆 1 of 𝐻:

𝜆 0 + 𝜇 min -𝑘Φ(𝑊) ≤ 𝜆 1 ≤𝜆 0 + 𝜇 min , (1.9)
where 𝜇 min denotes the smallest eigenvalue of the matrix 𝑃𝑊 𝑃.

Remark 4 (Non-self-adjoint 𝐻). With the sacrifice of the explicit constants in (1.6) and (1.7) (mostly coming from (2.5)), the self-adjointness assumption on 𝐻 can be removed. However, for the problem of quantum resonances, one can still obtain explicit estimates.

In the rest of this section 𝐻 is an abstract operator not necessarily self-adjoint or of the form (1.1). Our approach is grounded in the following theorem (see [START_REF] Gustafson | Mathematical Concepts of Quantum Mechanics[END_REF], Theorem 11.1). Theorem 1.2. Let 𝐻 be an operator on a Hilbert space and 𝑃 and 𝑃 ⊥ , a pair of projections such that 𝑃 + 𝑃 ⊥ = 1. Assume 𝐻 ⊥ := 𝑃 ⊥ 𝐻𝑃 ⊥ is invertible on Ran 𝑃 ⊥ and the expression

𝐹 𝑃 (𝐻) := 𝑃(𝐻 -𝐻𝑅 ⊥ 𝐻)𝑃, (1.10)
where 𝑅 ⊥ := 𝑃 ⊥ (𝐻 ⊥ ) -1 𝑃 ⊥ , defines a bounded operator. Then 𝐹 𝑃 , considered as a map on the space of operators, is isospectral in the following sense:

(a) 𝜆 ∈ 𝜎(𝐻) if and only if 0 ∈ 𝜎(𝐹 𝑃 (𝐻 -𝜆));

(b) 𝐻𝜓 = 𝜆𝜓 if and only if 𝐹 𝑃 (𝐻 -𝜆) 𝜑 = 0; (c) dim 𝑁𝑢𝑙𝑙 (𝐻 -𝜆) = dim 𝑁𝑢𝑙𝑙𝐹 𝑃 (𝐻 -𝜆).
(d) 𝜓 and 𝜑 in (b) are related as 𝜑 = 𝑃𝜓 and 𝜓 = 𝑄 𝑃 (𝜆)𝜑, where

𝑄 𝑃 (𝜆) := 𝑃 -𝑃 ⊥ (𝐻 ⊥ -𝜆) -1 𝑃 ⊥ 𝐻𝑃. (1.11)
Finally, 𝐹 𝑃 (𝐻) * = 𝐹 𝑃 * (𝐻 * ) and therefore, if 𝐻 and 𝑃 are self-adjoint, then so is 𝐹 𝑃 (𝐻).

A proof of this theorem is elementary and short; it can be found in [START_REF] Bach | Quantum electrodynamics of confined non-relativistic particles[END_REF], Section IV.1, pp 346-348, and [START_REF] Gustafson | Mathematical Concepts of Quantum Mechanics[END_REF], Appendix 11.6, pp 123-125.

The map 𝐹 𝑃 on the space of operators, is called the Feshbach-Schur map. The relation 𝜓 = 𝑄 𝑃 (𝜆)𝜑 allows us to reconstruct the full eigenfunction from the projected one. By statement (a), we have (1.12)

Concentrating on the eigenvalue problem, Corollary 1.3 shows that the original problem

𝐻𝜓 = 𝜆𝜓, (1.13)
is mapped into the equivalent eigenvalue problem,

𝐻 (𝜆)𝜑 = 𝜆𝜑, (1.14)
nonlinear in the spectral parameter 𝜆, but on the smaller space Ran 𝑃.

Since the projection 𝑃 is of a finite rank, the original eigenvalue problem (1.13) is mapped into an equivalent lower-dimensional/finite-dimensional one, (1.14). Of course, we have to pay a price for this: at one step we have to solve a one-dimensional fixed point problem that can be equivalently seen as a non-linear eigenvalue problem and invert an operator in Ran 𝑃 ⊥ .

We call this approach the Feshbach-Schur map method, or FSM method, for short. It is rather compact, as one easily see skimming through this paper and entirely elementary.

We call 𝐻 (𝜆) the effective Hamiltonian (matrix) and write as

𝐻 (𝜆) = 𝑃𝐻𝑃 + 𝑈 (𝜆), (1.15)
with the self-adjoint effective interaction, or a Schur complement, 𝑈 (𝜆), defined as

𝑈 (𝜆) := -𝑃𝐻𝑃 ⊥ (𝐻 ⊥ -𝜆) -1 𝑃 ⊥ 𝐻𝑃. (1.16)
It is shown in Lemma 2.1 below that (1.16) defines a bounded operator family.

We mention here some additional properties discussed in [START_REF] Dusson | Analysis of the Feshbach-Schur method for the planewave discretizations of Schrödinger operators[END_REF].

Proposition 1.4. Let 𝐻 be self-adjoint, let Λ be the same as in Corollary 1.3, let 𝑚 be the rank of 𝑃 and let the eigenvalues of 𝐻 (𝜆), 𝜆 ∈ Λ ∩ R, be labeled in the order of their increase counting their multiplicities so that 𝜈 1 (𝜆) ≤ . . . ≤ 𝜈 𝑚 (𝜆).

(1.17)

Then we have that (a) a solution of the equation 𝜈 𝑖 (𝜆) = 𝜆, for 𝜆 ∈ Λ ∩ R, is the 𝑖-th eigenvalue, 𝜆 𝑖 , of 𝐻 and vice versa; (b) 𝜈 𝑖 is differentiable in 𝜆 and 𝜈 𝑖 (𝜆)≤0, for 𝜆 ∈ Λ ∩ R.

Proof. (a) By (1.17), 𝜆 𝑖 = 𝜈 𝑖 (𝜆 𝑖 ) < 𝜈 𝑖+1 (𝜆) = 𝜆 𝑖+1 (except for the level crossings), which proves the result. For (b), by the explicit formula 𝑈 (𝜆) := -𝑃𝐻𝑃 ⊥ (𝐻 ⊥ -𝜆) -2 𝑃 ⊥ 𝐻𝑃 ≤ 0, we have 𝑈 (𝜆) ≤ 0, which by the Hellmann-Feynman theorem (cf. (2.19)) implies 𝜈 𝑖 (𝜆) ≤ 0. Remark 5 (Perturbation expansion). In the context of Hamiltonians of form (1.1) satisfying Assumptions (A)-(C), the FSM leads to a perturbation expansion to an arbitrary order. Indeed, in this case, 𝑃 is the orthogonal projection onto Null(𝐻 0 -𝜆 0 ), 𝑃 ⊥ := 1 -𝑃 and 𝑈 (𝜆) can be written as

𝑈 (𝜆) := -𝑃𝑊 𝑃 ⊥ (𝐻 ⊥ -𝜆) -1 𝑃 ⊥ 𝑊 𝑃.
Now, using the notation 𝐴 ⊥ := 𝑃 ⊥ 𝐴𝑃 ⊥ | Ran 𝑃 ⊥ and expanding

(𝐻 ⊥ -𝜆) -1 = (𝐻 ⊥ 0 + 𝑊 ⊥ -𝜆) -1
in 𝑊 ⊥ and at the same time iterating fixed point equation (1.12), we generate a perturbation expansion for eigenvalues 𝐻 to an arbitrary order (see also Remark 2 above). The paper is organized as follows. In Section 2, we present the proof of the main result, Theorem 1.1. Section 3 uses this result to obtain bounds of the ground-state energy of Helium-type ions.

Perturbation estimates

We want to use Theorem 1.2(b, c) to reduce the original eigenvalue problem to a simpler one. In this section, we assume that Conditions (A)-(C) of Section 1 are satisfied. Recall that 𝛾 0 is the distance of 𝜆 0 to the rest of the spectrum of 𝐻 0 and the expression Φ(𝑊) is defined in (1.2). First, we prove that the operator 𝐹 𝑃 (𝐻 -𝜆) is well-defined for 𝜆 ∈ Ω, where, with 𝑎 the same as in Theorem 1.1,

Ω := {𝑧 ∈ C : |𝑧 -𝜆 0 | ≤ 𝑎𝛾 0 }. (2.1)
Recall that 𝑃 denotes the orthogonal projection onto Null(𝐻 0 -𝜆 0 ) and 𝑃 ⊥ := 1 -𝑃. Denote

𝐻 ⊥ := 𝑃 ⊥ 𝐻𝑃 ⊥ | Ran 𝑃 ⊥ and 𝑅 ⊥ (𝜆) := 𝑃 ⊥ (𝐻 ⊥ -𝜆) -1 𝑃 ⊥
and recall 𝑘 = 1 1-𝑎-𝑏 . We have Lemma 2.1. Recall 𝜆 * := 𝜆 0 + 𝛾 0 and assume (1.3). Then, for 𝜆 ∈ Ω, the following statements hold

(a) The operator 𝐻 ⊥ -𝜆 is invertible on Ran 𝑃 ⊥ ; (b) The inverse 𝑅 ⊥ (𝜆) := 𝑃 ⊥ (𝐻 ⊥ -𝜆) -1 𝑃 ⊥ defines a bounded, analytic operator- family; (c) The expression 𝑈 (𝜆) := -𝑃𝐻𝑅 ⊥ (𝜆)𝐻𝑃 (2.2)
defines a finite-rank, analytic operator-family, bounded as

𝑈 (𝜆) ≤ 𝑘Φ(𝑊). (2.3) (d) 𝑈 (𝜆) is symmetric for any 𝜆 ∈ Ω ∩ R and therefore 𝐻 (𝜆) is self-adjoint as well.
Proof. With the notation 𝐴 ⊥ := 𝑃 ⊥ 𝐴𝑃 ⊥ | Ran 𝑃 ⊥ , we write

𝐻 ⊥ = 𝐻 ⊥ 0 + 𝑊 ⊥ .
To prove (a), we let 𝐻 𝜆 := 𝐻 ⊥ 0 -𝜆 and use the factorization

𝐻 𝜆 = |𝐻 𝜆 | 1 2 𝑉 𝜆 |𝐻 𝜆 | 1 2
, where 𝑉 𝜆 is a unitary operator, and use that, for 𝜆 ∈ Ω, the operator 𝐻 𝜆 is invertible and therefore we have the identity

𝐻 ⊥ -𝜆 = |𝐻 𝜆 | 1 2 [𝑉 𝜆 + 𝐾 𝜆 ] |𝐻 𝜆 | 1 2 , (2.4)
where

𝐾 𝜆 := |𝐻 𝜆 | -1 2 𝑊 ⊥ |𝐻 𝜆 | -1 2 . Next, for 𝜆 ∈ Ω, we have 𝐻 ⊥ 0 |𝐻 𝜆 | -1 𝑃 ⊥ = 𝑓 (𝜆), where
𝑓 (𝜆) := sup 𝑠 𝑠 -𝜆 : 𝑠 ≥ 0, |𝑠 -𝜆 0 | ≥ 𝛾 0 .
Assuming |𝜆 -𝜆 0 | ≤ 𝑎𝛾 0 and using

𝑠 𝑠 -𝜆 ≤ 1 + 𝜆 𝑠 -𝜆 ≤ 1 + 𝜆 0 + 𝑎𝛾 0 (1 -𝑎)𝛾 0 = 𝜆 * (1 -𝑎)𝛾 0 , we obtain 𝑓 (𝜆) ≤ 𝜆 * (1 -𝑎)𝛾 0 .
Since 𝐻

1 2 0 |𝐻 𝜆 | -1 2 𝑃 ⊥ = (𝐻 0 |𝐻 𝜆 | -1 | Ran 𝑃 ⊥ ) 1 2 𝑃 ⊥ , we have for 𝜆 ∈ Ω, 𝐻 1 2 0 |𝐻 𝜆 | -1 2 𝑃 ⊥ ≤ 𝜆 * (1 -𝑎)𝛾 0 1 2 , ( 2.5)
which implies in particular that 𝐾 𝜆 ≤ 𝜆 * (1-𝑎) 𝛾 0 𝑊 ⊥ 𝐻 0 . By the assumption (1.3), i.e., 𝑊 ⊥ 𝐻 0 ≤ 𝑏𝛾 0 𝜆 * , we have

𝐾 𝜆 ≤ 𝑏 1 -𝑎 .
(2.6)

Since 1 -𝑎 > 𝑏, by (2.4), the operator 𝐻 ⊥ -𝜆 is invertible and its inverse is analytic in 𝜆 ∈ Ω, which proves (a) and (b). We show that statement (c) is also satisfied. Since 𝑃𝐻 0 = 𝐻 0 𝑃 and 𝑃𝑃 ⊥ = 0, we have

𝑃𝐻𝑃 ⊥ = 𝑃𝑊 𝑃 ⊥ , 𝑃 ⊥ 𝐻𝑃 = 𝑃 ⊥ 𝑊 𝑃.
These relations and definition (2.2) yield

𝑈 (𝜆) = -𝑃𝑊 𝑅 ⊥ (𝜆)𝑊 𝑃.
(2.7)

Inverting (2.4) on Ran 𝑃 ⊥ and recalling the notation

𝑅 ⊥ (𝜆) := 𝑃 ⊥ (𝐻 ⊥ -𝜆) -1 𝑃 ⊥ gives 𝑅 ⊥ (𝜆) = |𝐻 𝜆 | -1 2 𝑃 ⊥ [𝑉 𝜆 + 𝐾 𝜆 ] -1 𝑃 ⊥ |𝐻 𝜆 | -1 2 .
(2.8)

Now, using identity (2.8), estimate (2.5) and (2.6), we find, for 𝜆 ∈ Ω,

𝐻 1 2 0 𝑅 ⊥ (𝜆)𝐻 1 2 0 ≤ 𝑘𝜆 * 𝛾 0 . (2.9)
Furthermore, by the eigen-equation 𝐻 0 𝑃 = 𝜆 0 𝑃, we have

𝐻 1 2 0 𝑃 2 = 𝜆 0 . (2.10)
Using expression (2.7) and estimates (2.9) and (2.10), we arrive at inequality (2.3). The analyticity follows from (2.7) and the analyticity of 𝑅 ⊥ (𝜆). For (d), since 𝐻 0 , 𝑊 and 𝑃 are self-adjoint, then so are 𝑈 (𝜆), for any 𝜆 ∈ Ω ∩ R, and, since 𝑈 (𝜆) is bounded, 𝐻 (𝜆) is self-adjoint as well.

Proof of Theorem 1.1. Let Ω be given by equation (2.1). Recall that, by Lemma 2.1 and Theorem 1.2, the 𝑚 × 𝑚 matrix-family 𝐻 (𝜆) := 𝐹 𝑃 (𝐻 -𝜆) + 𝜆𝑃, with 𝐹 𝑃 given in (1.10), is well defined, for each 𝜆 ∈ Ω, and can be written as (1.15). Since 𝑃𝐻𝑃 = 𝜆 0 𝑃 + 𝑃𝑊 𝑃, Eq. (1.15) can be rewritten as

𝐻 (𝜆) = 𝜆 0 𝑃 + 𝑃𝑊 𝑃 + 𝑈 (𝜆).
(2.11)

Eqs (2.
3) and (2.11) imply the inequality

𝐻 (𝜆) -𝜆 0 𝑃 -𝑃𝑊 𝑃 ≤ 𝑘Φ(𝑊).
(2.12)

By a fact from Linear Algebra, for each 𝜆 ∈ Ω ∩ R, the total multiplicity of the eigenvalues of the 𝑚 × 𝑚 self-adjoint matrix 𝐻 (𝜆) is 𝑚.

Denote by 𝜈 𝑖 (𝜆), 𝑖 = 1, 2, . . . , 𝑚, the eigenvalues of 𝐻 (𝜆), repeated according to their multiplicities. Eq (2.12) yields

|𝜈 𝑖 (𝜆) -𝜆 0 | ≤ 𝑃𝑊 𝑃 + 𝑘Φ(𝑊).
(2.13) Indeed, let 𝑃 𝑖 (𝜆) be the orthogonal projection onto Null(𝐻 (𝜆) -𝜈 𝑖 (𝜆)). Then

𝐻 (𝜆)𝑃 𝑖 (𝜆) = 𝜈 𝑖 (𝜆)𝑃 𝑖 (𝜆),
which, due to (2.11) and 𝑃 𝑖 (𝜆)𝑃 = 𝑃 𝑖 (𝜆), can be rewritten as

(𝜈 𝑖 (𝜆) -𝜆 0 )𝑃 𝑖 (𝜆) = (𝑃𝑊 𝑃 + 𝑈 (𝜆))𝑃 𝑖 (𝜆).
Equating the operator norms of both sides of this equation and using (2.12) and 𝑃 𝑖 (𝜆) = 1 gives (2.13). By Corollary 1.3, the eigenvalues of 𝐻 in the interval Ω ∩ R are in one-to-one correspondence with the solutions of the equations

𝜈 𝑖 (𝜆) = 𝜆 (2.14)
in Ω ∩ R. If this equation has a solution, then, due to (2.13), this solution would satisfy (1.6). Thus, we address (2.14). Let

Ω := {𝑧 ∈ R : |𝑧 -𝜆 0 | ≤ 𝑟𝛾 0 }, (2.15)
with

𝑟 := 1 𝛾 0 𝑃𝑊 𝑃 + 𝑘Φ(𝑊) .
By our assumption (1.4), 𝑟 < 𝑎 < 1 -𝑏.

Recall that by the definition, a branch point is a point at which the multiplicity of one of the eigenvalue families (branches) changes. One could think on a branch point as a point where two or more distinct eigenvalue branches intersect. Our next result shows that the eigenvalue branches of 𝐻 (𝜆) could be chosen in a differentiable way and estimates their derivatives.

Proposition 2.2. The following statements hold.

i) The eigenvalues 𝜈 𝑖 (𝜆) of 𝐻 (𝜆) and the corresponding eigenfunctions can be chosen to be differentiable for 𝜆 ∈ Ω .

ii) The derivatives 𝜈 𝑖 (𝜆), 𝜆 ∈ Ω , are bounded as

|𝜈 𝑖 (𝜆)| ≤ 𝑘 𝑎 -𝑟 Φ(𝑊) 𝛾 0 .
(2.16)

iii) 𝜈 𝑖 (𝜆) maps the interval Ω into itself.

Consequently, since by (1.5) the r.h.s. of (2.16) is < 1, the equations 𝜈 𝑖 (𝜆) = 𝜆 have unique solutions in Ω .

Proof. Proof of (i) for simple 𝜆 0 . For a simple eigenvalue 𝜆 0 , 𝑃 is a rank-one projection on the space spanned by the eigenvector 𝜑 0 of 𝐻 0 corresponding to the eigenvalue 𝜆 0 and therefore Eq. (2.11) implies that 𝐻 (𝜆) = 𝜈(𝜆)𝑃, with

𝜈(𝜆) := 𝜆 0 + 𝜑 0 , (𝑊 + 𝑈 (𝜆))𝜑 0 .
This and Lemma 2.1 show that the eigenvalue 𝜈(𝜆) is analytic. Proof of (i) for degenerate 𝜆 0 . We pick an arbitrary point 𝜇 in Ω and let 𝑃 𝜇𝑖 be the orthogonal projections onto the eigenspaces of 𝐻 (𝜇) corresponding to the eigenvalues 𝜈 𝑖 (𝜇), i.e. for a fixed 𝑖, 𝑃 𝜇𝑖 projects on the spans of all eigenvectors with the same eigenvalue 𝜈 𝑖 (𝜇). We now show that, in a neighbourhood of 𝜇, the eigenfunctions 𝜒 𝑖 (𝜆) can be chosen in a differentiable way. We introduce the system of equations for the eigenvalues 𝜈 𝑖 (𝜆) and corresponding eigenfunctions 𝜒 𝑖 (𝜆):

𝜈 𝑖 (𝜆) = 𝜒 𝑖 (𝜆), 𝐻 (𝜆) 𝜒 𝑖 (𝜆) 𝜒 𝑖 (𝜆) 2
, (2.17)

𝜒 𝑖 (𝜆) = 𝜒 𝑖 (𝜇) -𝑅 ⊥ (𝜈 𝑖 (𝜆), 𝜆)𝑊 𝜆 𝜒 𝑖 (𝜇), (2.18)
where

𝑅 ⊥ (𝜈, 𝜆) := 𝑃 ⊥ 𝜇𝑖 (𝑃 ⊥ 𝜇𝑖 𝐻 (𝜆)𝑃 ⊥ 𝜇𝑖 -𝜈) -1 𝑃 ⊥ 𝜇𝑖 , 𝑊 𝜆 := 𝑈 (𝜆) -𝑈 (𝜇).
The expression on the right side of (2.18) is the 𝑄-operator for 𝐻 (𝜆) and 𝑃 𝜇𝑖 defined according to (1.11) For each 𝑖, we can reduce system (2.17)-(2.18) to the single equation for 𝜈 𝑖 (𝜆), by treating 𝜒 𝑖 (𝜆) in (2.17) as given by (2.18). This leads to the fixed point problem for the functions 𝜈 𝑖 (𝜆):

𝜈 = 𝐹 𝑖 (𝜈, 𝜆),
where

𝐹 𝑖 (𝜈, 𝜆) := 𝜒 𝑖 (𝜈, 𝜆), 𝐻 (𝜆) 𝜒 𝑖 (𝜈, 𝜆) 𝜒 𝑖 (𝜈, 𝜆) 2 , 𝜒 𝑖 (𝜈, 𝜆) := 𝜒 𝑖 (𝜇) -𝑅 ⊥ (𝜈, 𝜆)𝑊 𝜆 𝜒 𝑖 (𝜇).
Notice that since 𝑃 ⊥ 𝜇𝑖 𝐻 (𝜆)𝑃 ⊥ 𝜇𝑖 are self-adjoint, the resolvent 𝑅 ⊥ (𝜈, 𝜆) and its derivatives in 𝜈 are uniformly bounded in a neighbourhood of (𝜈 𝑖 (𝜇), 𝜇), which does not contain branch points, except, possibly, for 𝜇.

To be more specific, the resolvent 𝑅 ⊥ (𝜈, 𝜆) is uniformly bounded in a neighbourhood Ω 𝜇𝑖 of (𝜈 𝑖 (𝜇), 𝜇) whether 𝜇 is a branch point or not, as long as Ω 𝜇𝑖 does not contain any other branch point than 𝜇.

Hence 𝐹 𝑖 (𝜈, 𝜆) is differentiable in 𝜈 and 𝜆 and |𝜕 𝜈 𝐹 𝑖 (𝜈, 𝜆)| → 0, as 𝜆 → 𝜇 (since 𝜕 𝜈 𝜒 𝑖 (𝜈, 𝜆) = 𝑅 ⊥ (𝜈, 𝜆) 2 𝑊 𝜆 𝜒 𝑖 (𝜇) and therefore

𝜕 𝜈 𝜒 𝑖 (𝜈, 𝜆) ≤ 𝑅 ⊥ (𝜈, 𝜆) 2 (𝐻 0 + 𝛼) 1 2 𝑊 𝜆 𝐻 0 , 𝛼 (𝐻 0 + 𝛼) 1 2 𝜒 𝑖 (𝜇) → 0,
as 𝜆 → 𝜇). Moreover, 𝜈 𝑖 (𝜇) -𝐹 𝑖 (𝜈 𝑖 (𝜇), 𝜇) = 0. Let 𝑓 𝑖 (𝜈, 𝜆) := 𝜈 -𝐹 𝑖 (𝜈, 𝜆). Then, by the above, 𝑓 𝑖 (𝜈 𝑖 (𝜇), 𝜇) = 0 and |𝜕 𝜈 𝑓 𝑖 (𝜈, 𝜆) -1| → 0, as 𝜆 → 𝜇. Thus, the implicit function theorem is applicable to the equation 𝑓 𝑖 (𝜈, 𝜆) = 0 and shows that there is a unique solution 𝜈 𝑖 (𝜆) in a neighbourhood of (𝜈 𝑖 (𝜇), 𝜇) and that this solution is differentiable in 𝜆.

Next (𝜆) . Writing equation (2.17) as

𝜈 𝑖 (𝜆) = χ𝑖 (𝜆), 𝐻 (𝜆) χ𝑖 (𝜆)
and differentiating this with respect to 𝜆, we obtain

𝜈 𝑖 (𝜆) = χ𝑖 (𝜆), 𝐻 (𝜆) χ𝑖 (𝜆) + 𝜂(𝜆),
where 𝜂(𝜆) := χ 𝑖 (𝜆), 𝐻 (𝜆) χ𝑖 (𝜆) + χ𝑖 (𝜆), 𝐻 (𝜆) χ 𝑖 (𝜆) . Now, moving 𝐻 (𝜆) in the last term to the l.h.s. and using that 𝐻 (𝜆) χ𝑖 (𝜆) = 𝜈 𝑖 (𝜆) χ𝑖 (𝜆), we find

𝜂(𝜆) = 𝜈 𝑖 (𝜆) [ χ 𝑖 (𝜆), χ𝑖 (𝜆) + χ𝑖 (𝜆), χ 𝑖 (𝜆) ] = 𝜈 𝑖 (𝜆)𝜕 𝜆 χ𝑖 (𝜆), χ𝑖 (𝜆) = 0,
|𝜈 𝑖 (𝜆) -𝜆 0 | ≤ 𝑟𝛾 0 .
By (2.15), this shows that 𝜈 𝑖 (𝜆) maps the interval Ω into itself which proves (iii). Hence, under condition (1.4), the fixed point equations 𝜆 = 𝜈 𝑖 (𝜆) have unique solutions on the interval Ω , proving Proposition 2.2. By Corollary 1.3, the eigenvalues of 𝐻 in Ω are in one-to-one correspondence with the solutions of the equations 𝜈 𝑖 (𝜆) = 𝜆. By Proposition 2.2, these equations have unique solutions, say, 𝜆 𝑖 . Then, estimate (2.13) implies inequality (1.6).

To obtain estimate (1.7), we recall from Theorem 1.2 that 𝑄(𝜆 𝑗 )𝜑 0 𝑗 = 𝜑 𝑗 , where the operator 𝑄(𝜆) is given by 𝑄(𝜆) := 1 -𝑅 ⊥ (𝜆)𝑃 ⊥ 𝑊 𝑃, 𝜆 𝑗 are the eigenvalues of 𝐻 and 𝜑 0 𝑗 are eigenfunctions of 𝐻 0 corresponding to 𝜆 0 . This gives

𝜑 0 𝑗 -𝜑 𝑗 = 𝜑 0 𝑗 -𝑄(𝜆 𝑗 )𝜑 0 𝑗 = 𝑅 ⊥ (𝜆 𝑗 )𝑃 ⊥ 𝑊 𝑃. (2.22)
Now, as in the derivation of (2.9), using identity (2.8), estimates (2.5) and 𝜒 𝑖 (𝜆) = -𝑅 ⊥ (𝜈 𝑖 (𝜆), 𝜆)𝑈 (𝜆) 𝜒 𝑖 (𝜆).

|𝐻 𝜆 | -
Finally we prove the relations (1.8) and (1.9) stated in the introduction.

Proof of (1.8). Eq. (2.11) gives

𝐻 (𝜆) = (𝜆 0 + 𝑊 + 𝑈 (𝜆) )𝑃,
where we use the notation 𝐴 := 𝜑 0 , 𝐴𝜑 0 . Since 𝑃 is a rank one projection, it follows that 𝐻 (𝜆) has only one eigenvalue and this eigenvalue is

𝜈 1 (𝜆) = 𝜆 0 + 𝑊 + 𝑈 (𝜆) .
This formula and estimate (2.3) show that 𝜈 1 (𝜆) obeys the estimate

|𝜈 1 (𝜆) -𝜆 0 -𝑊 | ≤ 𝑘Φ(𝑊).
By Proposition 2.2, the eigenvalue 𝜆 1 of 𝐻 satisfies the equation 𝜆 1 = 𝜈 1 (𝜆 1 ). Then the estimate above implies inequality (1.8).

Proof of Eq (1.9). Let 𝑊 be symmetric and denote by 𝜇 min the smallest eigenvalue of the matrix 𝑃𝑊 𝑃. Then 𝑈 (𝜆) ≤ 0 and (2.11) implies 𝐻 (𝜆) ≤ 𝜆 0 𝑃 + 𝑃𝑊 𝑃. This, together with inf 𝐴 ≤ inf 𝐵 for 𝐴 ≤ 𝐵, gives the upper bound 𝜈 1 (𝜆) ≤ 𝜆 0 + 𝜇 min on the smallest eigenvalue-branch 𝜈 1 (𝜆) of 𝐻 (𝜆).

For the lower bound, Eqs (2.3) and (2.11) imply the following estimate

𝜆 0 + 𝜇 min -𝑘Φ(𝑊) ≤ 𝜈 1 (𝜆).
The last two estimates and the equation 𝜈 𝑖 (𝜆) = 𝜆 (see (1.12)) yield (1.9).

Application: The ground state energy of the Heliumtype ions

In this section, we will use the inequalities obtained above to estimate the ground state energy of the Helium-type ions, which is the simplest not completely solvable atomic quantum system. For simplicity, we assume that the nucleus is infinitely heavy, but we allow for a general nuclear charge 𝑧𝑒, 𝑧 ≥ 1. Then the corresponding Schrödinger operator (describing 2 electrons of mass 𝑚 and charge -𝑒, and the nucleus of infinite mass and charge 𝑧𝑒) is given by ).

𝐻 𝑧,2 = 2 ∑ 𝑗=1 - ℏ 2 2𝑚 Δ 𝑥 𝑗 - 𝜅𝑧𝑒 2 |𝑥 𝑗 | + 𝜅𝑒 2 |𝑥 1 -𝑥 2 | , ( 3
For computations of the relativistic and QED contributions, see [START_REF] Pachucki | Testing fundamental interactions on the helium atom[END_REF][START_REF] Turbiner | Few-electron atomic ions in non-relativistic QED: The ground state[END_REF][START_REF] Yerokhin | Theoretical energies of low-lying states of light helium-like ions[END_REF]. Now, we derive some consequences of estimates (3.2). Let 𝑧 * be the smallest 𝑧 for which the error bound in the symmetric case is less than or equal to the smallest explicit (subleading) term 𝑤 1 . Inequality (3.2) shows that the latter bound is satisfied for 𝑧 such that 𝑤 1 𝑧 ≥ 10𝑤 2 𝛾 0 , which shows that 𝑧 * = 10𝑤 2 𝑤 1 𝛾 0 and consequently, 𝑧 * ≈ 12. According to equation (3.2), the symmetric ground state energy is lower than the anti-symmetric one, if -𝑧 2 + 𝑤 1 𝑧 < - 5 8 𝑧 2 + 𝑤 as 1 𝑧 -160 𝑤 as 2 . Using the values 𝑤 1 ≈ 0.6, 𝑤 2 ≈ 0.3, 𝑤 as 1 ≈ 0.20 and 𝑤 as 2 ≈ 0.01, we find that, based on (3.2), the ground state is symmetric and therefore its spin is 0 for

𝑧 ≥ 4 3 2 5 + √ 4 25 + 96 5 = 8 3 ≈ 2.6.
We conjecture that the symmetric ground state energy is lower than the anti-symmetric one for all 𝑧's.

Finally, note that on symmetric functions, the eigenvalue of the unperturbed operator is simple and on anti-symmetric ones, has the degeneracy 4, see below.

Proof of Proposition 3.2. First, we rescale the Hamiltonian (3.1) as 𝑥 → 𝑥/𝜇, with ) , where 𝑈 𝜇 : 𝜓(𝑥) → 𝜇 3 𝜓(𝜇𝑥) and 𝐻 (𝑧) is the rescaled Hamiltonian given by

𝜇 := 𝑧𝜅𝑒 2 𝑚 ℏ 2 = 𝑧 the Bohr radius , to obtain 𝑈 -1 𝜇 𝐻 𝑧,2 𝑈 𝜇 = 𝑧 2 𝜅 2 𝑒 4 𝑚 ℏ 2 𝐻 (𝑧
𝐻 (𝑧) = 2 ∑ 𝑗=1 - 1 2 Δ 𝑥 𝑗 - 1 |𝑥 𝑗 | + 1/𝑧 |𝑥 1 -𝑥 2 | .
Thus, it suffices to estimate the ground state energy of 𝐻 (𝑧) . We consider

𝑊 = 1/𝑧 |𝑥 1 -𝑥 2 |
as the perturbation and 2 𝑗=1 (-

1 2 Δ 𝑥 𝑗 -1 |𝑥 𝑖 |
) as the unperturbed operator. First, we consider the rescaled Hamiltonian 𝐻 (𝑧) on the symmetric functions subspace. On symmetric functions, the ground state energy of 2 𝑗=1 (-

1 2 Δ 𝑥 𝑗 -1 |𝑥 𝑗 |
) is -1 (see below), so we shift the operator 𝐻 (𝑧) by 1 + 𝛽, for some 𝛽 > 0, so that

𝐻 (𝑧) + 1 + 𝛽 = 𝐻 0 + 𝑊, with 𝐻 0 := 2 ∑ 𝑗=1 - 1 2 Δ 𝑥 𝑗 - 1 |𝑥 𝑗 | + 1 + 𝛽 and 𝑊 = 1/𝑧 |𝑥 1 -𝑥 2 | .
Now, the ground state energy of 𝐻 0 is 𝛽 and we can use inequality (1.8) and Proposition 1.4(c) to estimate the ground state energy of 𝐻 (𝑧) + 1 + 𝛽. By the HVZ theorem, the spectrum of 𝐻 0 on symmetric functions consists of the continuum [𝑒 1 + 1 + 𝛽, ∞) and the eigenvalues 𝑒 𝑚 + 𝑒 𝑛 + 1 + 𝛽, 𝑚, 𝑛 ≥ 1, where 𝑒 𝑛 , 𝑛 = 1, 2, . . . , denote the discrete eigenvalues of the Hydrogen Hamiltonian -1 2 Δ 𝑥 -1 |𝑥 | . The eigenvalues 𝑒 𝑛 , 𝑛 = 1, 2, . . . , are known explicitly:

𝑒 𝑛 = - 1 2𝑛 2 , (3.3)
with the multiplicities of 𝑛 2 and the corresponding eigenfunctions are given by Eq (A.1) below.

Then the ground state energy of 𝐻 0 is 𝜆 0 = 2𝑒 1 + 1 + 𝛽 = 𝛽, as claimed, and the gap,

𝛾 0 = 𝑒 1 + 𝑒 2 -2𝑒 1 = 𝑒 2 -𝑒 1 = 3/8.
Next, we show that the condition (1.3) is satisfied for 𝑧 ≥ 31.25. We begin with the really rough estimate:

1 |𝑥 1 -𝑥 2 | ≤ ℎ 𝑥 1 + ℎ 𝑥 2 + 10, (3.4)
where

ℎ 𝑥 := -1 2 Δ 𝑥 -1 |𝑥 | . First, we use Hardy's inequality, -Δ ≥ 1 4|𝑥 | 2 , and the estimate 1 4𝑚|𝑥 | 2 -1 |𝑥 | ≥ -𝑚 to obtain - 1 𝑚 Δ - 1 |𝑥| ≥ -𝑚. (3.5)
Next, passing to the relative and centre-of-mass coordinates, 𝑥 -𝑦 and 1 2 (𝑥 + 𝑦), we find

-Δ 𝑥 -Δ 𝑦 = -2Δ 𝑥-𝑦 - 1 2 Δ 1 2 ( 𝑥+𝑦) ≥ -2Δ 𝑥-𝑦 .
which, together with (3.5), yields

1 |𝑥 1 -𝑥 2 | ≤ - 1 2 Δ 𝑥 1 -𝑥 2 + 2 ≤ - 1 4 (Δ 𝑥 1 + Δ 𝑥 2 ) + 2. Now, we use -1 2 Δ 𝑥 = ℎ 𝑥 + 1 |𝑥 | ≤ ℎ 𝑥 -1 4 Δ 𝑥 + 4 to obtain - 1 4 Δ 𝑥 ≤ ℎ 𝑥 + 4.
The last two inequalities yield (3.4). Eq. (3.4), together with the relation

𝐻 0 = ℎ 𝑥 1 + ℎ 𝑥 2 + 1 + 𝛽,
implies the estimate:

1 |𝑥 1 -𝑥 2 | ≤ 𝐻 0 + 9 -𝛽. (3.6)
Now, we check the condition (1.3):

𝑃 ⊥ 𝑊 𝑃 ⊥ 𝐻 0 ≤ 𝑏𝛾 0 𝜆 * .
Using the last relation, we estimate

𝑧(𝐻 ⊥ 0 ) -1 2 𝑊 (𝐻 ⊥ 0 ) -1 2 ≤ (𝐻 ⊥ 0 ) -1 (𝐻 ⊥ 0 + 9 -𝛽) = 1 + (9 -𝛽) (𝐻 ⊥ 0 ) -1 .
Recalling that 𝜆 0 = 𝛽, we see that 𝐻 ⊥ 0 ≥ 𝛽 + 𝛾 0 , so that the last inequality gives

0 ≤ (𝐻 ⊥ 0 ) -1 2 𝑊 (𝐻 ⊥ 0 ) -1 2 ≤ 1 𝑧 9 + 𝛾 0 𝛽 + 𝛾 0 , provided 𝛽 < 9. This implies 𝑃 ⊥ 𝑊 𝑃 ⊥ 𝐻 0 ≤ 1 𝑧 9 + 𝛾 0 𝛽 + 𝛾 0 . (3.7)
Since

𝜆 * = 𝛽 + 𝛾 0 , condition (1.3) is satisfied, if 9+𝛾 0

𝑧

≤ 𝑏𝛾 0 , which gives 𝑧 ≥ 9+𝛾 0 𝑏𝛾 0 . Since 𝛾 0 = 3/8 and 𝑏 = 0.8, this implies 𝑧 ≥ 31.25. We will need the following lemma, whose proof is given after the proof of the proposition. Lemma 3.3. Recall that 𝑃 is the orthogonal projection onto the eigenspace of 𝐻 0 corresponding to the lowest eigenvalue 𝛽. We have

𝑊 = 𝑃𝑊 𝑃 = 𝑤 1 𝑧 and Φ(𝑊) ≤ 𝑤 2 𝑧 2 𝛾 0 , (3.8)
with, recall, 𝑊 := 𝜑 0 , 𝑊 𝜑 0 , where 𝜑 0 is the eigenfunction of 𝐻 0 corresponding to 𝜆 0 .

We check now the second condition of Theorem 1.1, (1.4). Recall the definition

𝑟 := 𝑃𝑊 𝑃 + 𝑘Φ(𝑊) 𝛾 0 .
Then condition (1.4) can be written as 𝑟 < 𝑎. By Lemma 3.3,

𝑟 ≤ 1 𝛾 0 𝑤 1 𝑧 + 𝑘𝑤 2 𝑧 2 𝛾 0 .
Recall the condition 𝑧 ≥ 31.25 and 𝑘 = 10, 𝑤 1 ≈ 0.6 and 𝑤 2 ≈ 0.3. Thus

𝑘𝑤 2 𝑧 2 𝛾 0 ≤ 𝑤 1 𝑧 𝑘 𝑤 2 31.25 𝛾 0 𝑤 1 ≤ 0.5 𝑤 1 𝑧 .
Then we find 𝑟 ≤ 1.5 𝑤 1 𝑧𝛾 0 and therefore condition (1.4) is satisfied if 𝑎 > 1.5 𝑤 1 𝑧𝛾 0 , or (for 𝑎 = 0.1)

𝑧 > 3 2 8 3 0.6 𝑎 = 24.
This is less restrictive than 𝑧 ≥ 31.25.

For the third condition (1.5), recalling that 𝜆 0 = 𝛽, so that 𝜆 * = 𝛾 0 + 𝛽 (= 𝜆 1 , the second eigenvalue of 𝐻 0 ) and using both relations in (3.8), we obtain that condition

(1.5) is satisfied if 𝑘𝑤 2 /(𝑧 2 𝛾 0 ) < 1 2 (𝑎𝛾 0 -𝑤 1 /𝑧), which is equivalent to 𝑧 > 1 2𝑎𝛾 0 𝑤 1 + √ 𝑤 2 1 + 8𝑎𝑘𝑤 2 .
Using the values 𝛾 0 = 3 8 , 𝑤 1 ≈ 0.6 and 𝑤 2 ≈ 0.3, and 𝑏 = 0.8 and 𝑎 = 0.1, giving 𝑘 = 10, this shows that the latter inequality, and therefore (1.5), holds if 𝑧 > 30, which is less restrictive than 𝑧 ≥ 31.25. Thus, all three conditions are satisfied for 𝑧 ≥ 31.25.

Next, we use (1.8) to estimate the ground state energy, 𝐸 (𝑧) sym of 𝐻 (𝑧) . The first and second relations in (3.8), together with (1.8) and the fact that 𝑈 (𝜆) ≤ 0, gives the following bounds

- 𝑘𝑤 2 𝛾 0 𝑧 2 ≤ 𝐸 (𝑧) sym + 1 - 𝑤 1 𝑧 ≤ 0.
which after rescaling gives (3.2). Now, we consider the ground state of 𝐻 (𝑧) on anti-symmetric functions. In this case, the ground state energy of 2 𝑗=1 (- Since 𝛾 as 0 = 5 72 and 𝑏 = 0.8, this implies 𝑧 ≥ 170 for the anti-symmetric case. We skip the verification of conditions (1.4) and (1.5), which are simpler and similar to that of the symmetric case.

1 2 Δ 𝑥 𝑗 -1 |𝑥 𝑖 | ) is 𝑒 1 + 𝑒 2 = -
To estimate the ground state energy, 𝐸 (𝑧) as , of 𝐻 (𝑧) on anti-symmetric functions, we use bound (1.9). First, we claim that Φ(𝑊) defined in (1.2) satisfies (cf. Remark 2)

Φ(𝑊) ≤ 𝜇 max 𝛾 as 0 , (3.10)
where 𝜇 max is the largest eigenvalue of the positive, 4 × 4-matrix 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃. Indeed, recall that

Φ(𝑊) := 𝜆 0 𝜆 * 𝑃 ⊥ 𝑊 𝑃 2
𝐴𝐴 * = 𝐻 -1 2 0 𝑃𝑊 𝑃 ⊥ 𝐻 -1 0 𝑊 𝑃𝐻 -1 2 0 , 𝐻 -1 2 0 𝑃 = 𝜆 -1 2 0 𝑃 (since 𝐻 0 𝑃 = 𝜆 0 𝑃),

𝑧

and

𝜇 max = 𝑤 as 2 𝑧 2 , (3.12)
for some positive 𝑧-independent constants 𝑤 as 1 and 𝑤 as 2 . These constants are computed in Appendix B. Now, we see that (1.9) and (3. Proof of Lemma 3.3. The ground state energy, 𝑒 1 , of -1 2 Δ 𝑥 -1 |𝑥 | is non-degenerate, with the normalized ground state (known as the 1s wavefunction) given by (see e.g. [START_REF] Messiah | Quantum Mechanics[END_REF], or Wiki, with 𝑎 0 replaced by 1 due to the rescaling above)

𝜓 1 (𝑥) = 1 √ 𝜋 𝑒 -| 𝑥 | .
Hence, the ground state energy, 0, of 𝐻 0 is also non-degenerate, with the normalized ground state 𝜓 1 (𝑥 1 )𝜓 1 (𝑥 2 ) (in the symmetric subspace). First, we compute 𝑊 and 𝑃𝑊 𝑃 . Let (𝜓 ⊗ 𝜙) (𝑥, 𝑦) := 𝜓(𝑥) 𝜙(𝑦) and, for an operator 𝐴 on 𝐿 2 (R 6 ), we denote 𝐴 := 𝜓 1 ⊗ 𝜓 1 , 𝐴𝜓 1 ⊗ 𝜓 1 . In the symmetric case, 𝑃𝑊 𝑃 = 𝑊 𝑃, where

𝑊 := 𝜓 1 ⊗ 𝜓 1 , 𝑊𝜓 1 ⊗ 𝜓 1 = 1 𝑧 ∫ R 6 |𝜓 1 (𝑥)𝜓 1 (𝑦)| 2 |𝑥 -𝑦| 𝑑𝑥 𝑑𝑦 =: 1 𝑧 𝑤 1 . (3.13)
Hence 𝑃𝑊 𝑃 = 𝑤 1 𝑧 . This gives the first relation in (3.8). Now, using that 𝑃 ⊥ = 1 -𝑃, we compute

𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 = 𝑃𝑊 2 𝑃 -𝑃𝑊 𝑃𝑊 𝑃 = 𝑤 2 𝑧 2 𝑃, (3.14)
and

𝑤 2 := 𝑧 2 𝑊 2 -𝑊 2 = ∫ R 6 |𝜓 1 (𝑥)𝜓 1 (𝑦)| 2 |𝑥 -𝑦| 2 𝑑𝑥 𝑑𝑦 - ∫ R 6 |𝜓 1 (𝑥)𝜓 1 (𝑦)| 2 |𝑥 -𝑦| 𝑑𝑥 𝑑𝑦 2 . ( 3

A The eigenfunctions of the Hydrogen-like Hamiltonian

The eigenfunctions 𝜓 𝑛ℓ 𝑘 (𝑥) are given by (see Wiki, "Hydrogen atom" with 𝑎 0 replaced by 1)

𝜓 𝑛ℓ 𝑘 (𝑥) = √ 2𝑧 𝑛 3 (𝑛 -ℓ -1)! 2𝑛(𝑛 + ℓ)! 𝑒 -𝜌/2 𝜌 ℓ 𝐿 2ℓ+1 𝑛-ℓ-1 (
1 := ∫ R 6 |𝜓 1 (𝑥)𝜓 1 (𝑦)| 2 |𝑥 -𝑦| 𝑑𝑦 𝑑𝑥,
Quoting from Wikipedia (https://en.wikipedia.org/wiki/Hydrogen_atom, 30.10.2020): "... the generalized Laguerre polynomials are defined differently by different authors. The usage here is consistent with the definitions used by [START_REF] Messiah | Quantum Mechanics[END_REF], p. 1136, and Wolfram Mathematica. In other places, the Laguerre polynomial includes a factor of (𝑛 + ℓ)!. ... or the generalized Laguerre polynomial appearing in the hydrogen wave function is 𝐿 2ℓ+1 𝑛+ℓ (𝜌) instead."

𝑤 2 := ∫ R 6 |𝜓 1 (𝑥)𝜓 1 (𝑦)| 2 |𝑥 -𝑦| 2 𝑑𝑦 𝑑𝑥 -𝑤 2 1 ,
it becomes clear that we need to compute terms of the form

𝐶 𝛼 = ∫ R 6 |𝜓 1 (𝑥)𝜓 1 (𝑦)| 2 |𝑥 -𝑦| 𝛼 𝑑𝑦 𝑑𝑥 = ∫ R 3 |𝜓 1 (𝑥)| 2 ∫ R 3 |𝜓 1 (𝑥 -𝑧)| 2 |𝑧| 𝛼 𝑑𝑧 𝑑𝑥,
with 𝛼 = 1, 2. We thus introduce a numerical quadrature in order to approximate the values of these integrals. We do not aim within this work to obtain the most efficient implementation. We define a numerical quadrature grid in terms of spherical coordinates with integration points 𝑥 𝑖,𝑛 and weights 𝜔 𝑖,𝑛 given by 𝑥 𝑖,𝑛 = 𝑟 𝑖 𝑠 𝑛 , 𝜔 𝑖,𝑛 = ℎ 𝑟 Then, we approximate 𝑤 1 ≈ 𝐶 1 and 𝑊 2 ≈ 𝐶 2 -𝐶 2 1 . We now shed our attention to the constants 𝑤 as 1 and 𝑤 as 2 and we first introduce the bi-electronic integrals, for (ℓ, 𝑘), (ℓ , 𝑘 ) = (0, 0), (1, -1), (1, 0), (1, 1), by Then, using the definition (3.9) and symmetry properties, we derive the following expressions for the matrix elements 2.

We then employ the same quadrature rule as above to approximate the matrix elements We then again, approximate the different integrals with the above introduced quadrature rule and find the maximal eigenvalue of Q in order to approximate 𝑤 as 2 . In Figure 1, we plot the the approximate value of 𝑤 1 , 𝑤 2 , 𝑤 as 1 and 𝑤 as 2 for different values of the parameters 𝑅 max , 𝑁 𝑟 and 𝑁 leb given in Table 2.

We

Corollary 1 . 3 .

 13 Assume there is an open set Λ ⊂ C such that 𝐻 -𝜆, 𝜆 ∈ Λ, is in the domain of the map 𝐹 𝑃 , i.e. 𝐹 𝑃 (𝐻 -𝜆) is well defined. Define the operator-family 𝐻 (𝜆) := 𝐹 𝑃 (𝐻 -𝜆) + 𝜆𝑃, and let 𝜈 𝑖 (𝜆), for 𝜆 in Λ, denote its eigenvalues counted with multiplicities. Then the eigenvalues of 𝐻 in Λ are in one-to-one correspondence with the solutions of the equations 𝜈 𝑖 (𝜆) = 𝜆.

with 1 20

 1 𝜆 0 = 𝛽 and observe that 𝑃𝑊 𝑃 ⊥ 𝐻 0 = 𝐴 , with 𝐴 := 𝐻 -𝑃𝑊 𝑃 ⊥ 𝐻 -

  and setting 𝑘 = 10 gives (3.2) in the antisymmetric case.

3 𝜓 3 |𝜓 1 ( 3 𝜓 1 ( 3 𝜓 1 (

 3313131 2ℓ 𝑘 (𝑥)𝜓 2ℓ 𝑘 (𝑥) ∫ R 𝑥)𝜓 2ℓ 𝑘 (𝑥) ∫ R 𝑦)𝜓 2ℓ 𝑘 (𝑦)|𝑥 -𝑦| 𝑑𝑦 𝑑𝑥.

M

  ℓ ,𝑘 ℓ,𝑘 = 𝑧 𝜓 2ℓ 𝑘 |𝑊 |𝜓 2ℓ 𝑘 = A ℓ ,𝑘 ℓ,𝑘 -B ℓ ,𝑘 ℓ,𝑘 .

Figure 1 :

 1 Figure 1: Value of the approximate coefficients 𝑤 1 , 𝑤 2 , 𝑤 as 1 and 𝑤 as 2 depending on the quadrature rule reported in Table2.

3 𝜓 3 |𝜓 1 ( 3 𝜓 1 ( 3 𝜓 1 (

 3313131 ℓ,𝑘 , B ℓ ,𝑘 ℓ,𝑘 and compute the smallest eigenvalue of M as an approximation of 𝑤 as 1 . Finally, the approximation of 𝑤 as 2 involves the matrix𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 = 𝑃𝑊 2 𝑃 -𝑃𝑊 𝑃𝑊 𝑃,whose elements are given byQ = N -M 2 with N ℓ ,𝑘 ℓ,𝑘 = 𝑧 2 𝜓 2ℓ 𝑘 |𝑊 2 |𝜓 2ℓ 𝑘 = C ℓ ,𝑘 ℓ,𝑘 -D ℓ 2ℓ 𝑘 (𝑥)𝜓 2ℓ 𝑘 (𝑥) ∫ R 𝑦)| 2 |𝑥 -𝑦| 2 𝑑𝑦 𝑑𝑥, 𝑥)𝜓 2ℓ 𝑘 (𝑥) ∫ R 𝑦)𝜓 2ℓ 𝑘 (𝑦) |𝑥 -𝑦| 2𝑑𝑦 𝑑𝑥.

  To bound the r.h.s. of (2.20), we use the analyticity of 𝑈 (𝜆) in Ω and estimate (2.3). Ω . Combining Eqs. (2.20) and (2.21), we arrive at(2.16).Due to the definition of 𝑟 after (2.15), we can rewrite estimate (2.13) as

	for 𝜆 ∈
	which implies (2.19).
	To prove Eq. (2.16), we use formula (2.19) above and the normalization of χ𝑖 (𝜆)
	to estimate 𝜈	𝑖 (𝜆) as
		|𝜈 𝑖 (𝜆)| ≤ 𝑈 (𝜆) .	(2.20)
	Indeed, by the Cauchy integral formula, we have
		𝑈 (𝜆) ≤	1	sup	𝑈 (𝜆 ) ,
			𝑅𝛾 0	|𝜆 -𝜆|=𝑅𝛾 0
				𝑘 𝑎 -𝑟	𝛾 0 Φ(𝑊)	,	(2.21)

with 𝑅 ≤ 𝑎 -𝑟, so that {𝜆 ∈ C : |𝜆 -𝜆| < 𝑅𝛾 0 } ⊂ Ω, for 𝜆 ∈ Ω . This, together with (2.3) and under the conditions of Lemma 2.1, gives the estimate 𝑈 (𝜆) ≤

  .1) acting on the space 𝐿 2 𝑠𝑦𝑚 (R 6 ) of 𝐿 2 -functions symmetric (or antisymmetric) w.r.t. the permutation of 𝑥 1 and 𝑥 2 . Here 𝜅 =1 4 𝜋 𝜀 0 is Coulomb's constant and 𝜀 0 is the vacuum By the Pauli principle, the product of coordinate and spin wave functions should antisymmetric w.r.t. permutation of the particle coordinates and spins. Hence, in the two particle case, after separation of spin variables, coordinate wave functions could be either antisymmetric or symmetric.𝛾 0 ) and the relative difference defined by err 𝑧 := |𝐸 exact -𝐸

	𝑧							10	20	30	40	50
	-𝐸 exact (from [16, 17, 18])		93.9	387.7	881.4	1575.2 2468.9
	main part -𝐸	sym,lead 𝑧,2	in (3.2)		94	388	882	1576	2470
	relative difference err 𝑧		0.11% 0.077% 0.068% 0.051% 0.045%
	relative error term Δ 𝑧 in (3.2) 8.51% 2.06%	0.91%	0.51%	0.32%
	Table 1: Comparison of non-rigorous computations (see [8, 9, 16, 17, 18]) with the
	main part -𝐸	sym,lead 𝑧,2	:= 𝑧 2 -𝑤 1 𝑧 in equation (3.2), its relative contribution of the
	error estimation Δ 𝑧 := 10𝑤 2 /(-𝐸	sym,lead 𝑧,2
			sym,lead 𝑧,2	|/(-𝐸	sym,lead 𝑧,2

  )𝜓 2ℓ 𝑘 (𝑥 2 ) -𝜓 2ℓ 𝑘 (𝑥 1 )𝜓 100 (𝑥 2 )),(3.9) for (ℓ, 𝑘) = (0, 0), (1, -1), (1, 0),[START_REF] Bach | Quantum electrodynamics of confined non-relativistic particles[END_REF][START_REF] Bach | Quantum electrodynamics of confined non-relativistic particles[END_REF], where 𝜓 𝑛ℓ 𝑘 (𝑥) for ℓ = 0, 1, 2, . . . , 𝑛 -1, 𝑘 = -ℓ, -ℓ+1, . . . , ℓ, 𝑛 = 1, 2, . . ., are the eigenfunctions of the Hydrogen-like Hamiltonian -1 2 Δ 𝑥 -1 | 𝑥 | corresponding to eigenvalues 𝑒 𝑛 , so that 𝜓 100 (𝑥) = 𝜓 1 (𝑥), see Appendix A. Now, we define the unperturbed operator as 𝑒 1 + 𝑒 3 -(𝑒 1 + 𝑒 2 ) = 𝑒 3 -𝑒 2 giving 𝛾 as 0 = 5 72 . For condition (1.3), note first that since now 𝐻 0 = ℎ 𝑥 1 + ℎ 𝑥 2 + 5

	5 8 of multiplicity 4, + 𝛽(< 𝑒 2 + 𝑒 2 + + 5 8 + 𝛽, 5 8 + 𝛽), 0 = 8 + 𝛽, Eq. (3.4) with the ground states 𝜙 2ℓ 𝑘 (𝑥 1 , 𝑥 2 ) := 1 √ 2 2 ∑ 𝑗=1 -1 2 Δ 𝑥 𝑗 -1 |𝑥 𝑗 | to have 𝛽 as the lowest eigenvalue of 𝐻 0 . For the gap, since the first two eigenvalues are 𝑒 1 + 𝑒 2 + 5 8 + 𝛽 = 𝛽 and 𝑒 1 + 𝑒 3 + 5 8 we have 𝛾 as implies 1 |𝑥 1 -𝑥 2 | ≤ 𝐻 0 + 9 + 3 8 -𝛽, instead of (3.6), which shows that (3.7) should be modified as 𝑃 ⊥ 𝑊 𝑃 ⊥ 𝐻 0 ≤ 1 𝑧 9 + 3 8 + 𝛾 as 0 𝛽 + 𝛾 as 0 , (𝜓 100 (𝑥 1 𝐻 0 := in the anti-symmetric case. Since 𝜆 as * = 𝛽 + 𝛾 as 0 , condition (1.3) is satisfied, if 1 𝑧 (9 + 3 8 + 𝛾 as 0 ) ≤ 𝑏𝛾 as 0 , which gives 𝑧 ≥ 9 + 3 8 + 𝛾 as 0 0 𝑏𝛾 as .

  and 𝐻 -1 0 𝑃 ⊥ ≤ 𝜆 -1 * 𝑃 ⊥ , we find𝐴𝐴 * ≤ (𝜆 0 𝜆 * ) -1 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃,which, together with 𝐴 = 𝐴 * = 𝐴𝐴 * 1 2 , gives𝑃𝑊 𝑃 ⊥ 2 𝐻 0 ≤ (𝜆 0 𝜆 * ) -1 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 .(3.11)Since 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 ≥ 0, we have 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 = 𝜇 max , which gives (3.10). Hence, we have to compute the smallest eigenvalue 𝜇 min of the 4 × 4 matrix 𝑃𝑊 𝑃 and the largest eigenvalue, 𝜇 max , of the positive, 4 × 4-matrix 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃, see(1.9).By scaling, it is easy to see that the matrices 𝑃𝑊 𝑃 and 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 are of the form 𝑃𝑊 𝑃 = 𝑧 -1 𝐴 and 𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 = 𝑧 -2 𝐵, where 𝐴 and 𝐵 are 𝑧-independent 4 × 4 matrices. Hence

	𝜇 min =	𝑤 as 1

  The second author is grateful to Volker Bach and Jürg Fröhlich for enjoyable collaboration. The authors are indebted to the anonymous referee for many useful suggestions and remarks.

			.15)
	By the Schwarz inequality and the normalization of 𝜓 1 , we have 𝑤 2 ≥ 0, which,
	together with (3.14), gives
	𝑃𝑊 𝑃 ⊥ 𝑊 𝑃 =	𝑤 2 𝑧 2 .	(3.16)
	Eqs. (3.11) and (3.16) imply the second relation in (3.8).

Acknowledgments.

  𝑁 𝑟 -1 and where {𝑠 𝑛 , 𝜔 leb 𝑛 } denote Lebedev integration points on the unit sphere with 𝑁 leb points. Thus, 𝑁 𝑟 denotes the number of points along the 𝑟-coordinate, 𝑅 max the radius of the furthest points away from the origin.We thus approximate 𝐶 𝛼 by𝐶 𝛼 ≈ |𝜓 1 (𝑥 𝑖,𝑛 -𝑥 𝑗,𝑚 )| 2 |𝑥 𝑗,𝑚 | 𝛼 .

	𝑁 𝑟 ∑	𝑁 leb ∑
	𝑖=1
	𝑁 𝑟 ∑	𝑁 leb ∑	𝜔 𝑗,𝑚
	𝑗=1	𝑚=1

3 𝑖 𝜔 leb 𝑛 with 𝑟 𝑖 = exp(-ln(𝑅 max ) + (𝑖 -1)ℎ), 𝑖 = 1, . . . , 𝑁 𝑟 , ℎ = 2 ln(𝑅 max ) 𝑛=1 𝜔 𝑖,𝑛 |𝜓 1 (𝑥 𝑖,𝑛 )| 2 Φ 𝛼 (𝑥 𝑖,𝑛 ) with Φ 𝛼 (𝑥 𝑖,𝑛 ) =

Table 2 :

 2 obtain approximate values𝑤 1 ≈ 0.625, 𝑤 2 ≈ 0.2628, 𝑤 as 1 ≈ 0.3, 𝑤 as 2 ≈ 0.01where the last reported digit is stable over the last three quadrature rules. Parameters for the different quadrature rules.

	index	1	2	3	4	5	6	7	8	9	10	11	12
	𝑅 max	16	18	20	22	24	26	28	30	32	34	34	34
	𝑁 𝑟	10	12	14	16	18	20	22	24	26	28	30	32
	𝑁 leb	194	266	350	590	974	1454	2030	2702	3470	4334	5294	5810

Supported in part by the French "Investissements d'Avenir" program, project ISITE-BFC (contract ANR-15-IDEX-0003).

permittivity. For 𝑧 = 2, 𝐻 𝑧,2 describes the Helium atom, for 𝑧 = 1, the negative ion of the Hydrogen and for 2 < 𝑧 ≤ 94 (or ≤ 103, depending on what one counts as stable elements), Helium-type positive ion. (We can call (3.1) with 𝑧 > 103 a 𝑧-ion.)

It is well known that 𝐻 𝑧,2 has eigenvalues below its continuum. Variational techniques give excellent upper bounds on the eigenvalues of 𝐻 𝑧,2 , but good lower bounds are hard to come by. Thus, we formulate Problem 3.1. Estimate the ground state energy of 𝐻 𝑧,2 .

The most difficult case is of 𝑧 = 1, the negative ion of the hydrogen, and the problem simplifies as 𝑧 increases.

Here we present fairly precise bounds on the ground state energy of 𝐻 𝑧,2 implied by our actual estimates. However, the conditions under which these estimates are valid impose rather sever restrictions in 𝑧. We introduce the reference energy

where, for symmetric functions, 𝑐 = 1, 𝛾 0 = 3 8 and 𝑤 1 ≈ 0.6 and 𝑤 2 ≈ 0.27, defined in equations (3.13) and (3.15), and, for anti-symmetric functions, 𝑐 as = 5 8 , 𝛾 as 0 = 5 72 and 𝑤 as 1 ≈ 0.20 and 𝑤 as 2 ≈ 0.01, defined in (3.12). The approximate values of 𝑤 1 , 𝑤 2 , 𝑤 as 1 and 𝑤 as 2 are computed numerically in Appendix B. Here we report the stable digits of our computations.

The inequalities 𝑧 ≥ 31.25 and 𝑧 ≥ 170 come from condition (1.3), while estimates (3.2), from (1.8), with 𝑏 = 0.8 and 𝑎 = 0.1 (which give 𝑘 = 10).

Table 1 compares the result for symmetric functions with computations in quantum chemistry. (We did not find results for the antisymmetric space.)

We observe that, except for 𝑧 = 50, the results of [START_REF] Turbiner | Helium-like and Lithium-like ions: Ground state energy[END_REF][START_REF] Turbiner | Ultra-compact accurate wave functions for He-like and Li-like iso?electronic sequences and variational calculus: I. Ground state[END_REF][START_REF] Turbiner | Few-electron atomic ions in non-relativistic QED: The ground state[END_REF] lie in the interval provided by the estimation (3.2).


Ph.D Marie-Odile Berger

M.D René Anxionnat
email: r.anxionnat@chu-nancy.fr-l.picard@chu-nancy.fr

Ph.D Erwan Kerrien
email: erwan.kerrien@loria.fr

M.D Luc Picard

M.D Michael Söderman
email: michael.soderman@karolinska.se

Projet Magrit

A methodology for validating a 3D imaging modality for brain AVM delineation: Application to 3DRA

Keywords: Angiography, Brain Arteriovenous Malformation, Radiosurgery, Validation, Principal Component Analysis

A general methodology is described to validate a 3D imaging modality with respect to 2D digital subtracted angiography (DSA) for brain AVMs (BAVM) delineation. It relies on the assessment of the statistical compatibility of the radiosurgical target delineated in 3D with its delineations in 2D. This methodology is demonstrated through a preliminary evaluation of rotational 3D angiography (3DRA). Generally speaking, BAVM delineation cannot be performed on 3DRA alone. However, in our study, 3DRA showed similar performances to DSA for rather easy cases, and even better for 3 patients. Conversely, 3 problematic cases are identified and discussed.

Introduction

Brain arteriovenous malformations (BAVMs) are complex pial or subpial arteriovenous shunts which can pose therapeutic problems. The major threat from a BAVM is brain hemorrhage. The natural course is still not very well known, but it seems that the risk of hemorrhage increases with BAVM volume, patient age, central location and prior hemorrhage. Most patients with BAVM are therefore subjected to treatment with surgery, embolization or radiosurgery (RS), unless the risk posed by the remedy itself is considered to be unacceptable.

RS is unique in that the outcome can be predicted from mathematical models [START_REF] Söderman | Volume determination and predictive models in the management of cerebral arteriovenous malformations[END_REF]. The chances for BAVM obliteration after RS depend on the dose of radiation given to the periphery of the target [START_REF] Karlsson | Can the probability for obliteration after radiosurgery for arteriovenous malformations be accurately predicted?[END_REF][START_REF] Bollet | Efficacy and morbidity of arc-therapy radiosurgery for cerebral arteriovenous malformations: a comparison with the natural history[END_REF]. Adverse radiation effects (ARE), in the worst case radionecrosis, constitute the main potential drawback. The risk for ARE depends on the location of the lesion and total amount of energy delivered, i.e. BAVM volume and radiation dose [START_REF] Karlsson | Can the probability for obliteration after radiosurgery for arteriovenous malformations be accurately predicted?[END_REF]. There is also a risk of hemorrhage during the latency period before obliteration, a risk that varies according to treatment parameters, AVM volume and patient age [START_REF] Karlsson | Risk for hemorrhage during the 2-year latency period following gamma knife radiosurgery for arteriovenous malformations[END_REF].

RS may be carried out with proton beams, the Gamma Knife (Elekta, Stockholm, Sweden) or a linear accelerator. In the latter case multileaf collimators are required to irradiate complex-shaped targets, such as BAVMs. The definition of the target volume is an essential stage of the treatment regardless of equipment. The totality of the nidus (which is the active part of the AVM) and as little as possible of the surrounding brain tissue should be included within the prescription isodose, in order to limit the risks of ARE. Indeed, the reasons for failed AVM obliteration are various but partial volume irradiation is put forward as a major cause by many authors [START_REF] Flickinger | A multi-institutional analysis of complication outcomes after arteriovenous malformation radiosurgery[END_REF][START_REF] Gallina | Failure in radiosurgery treatment of cerebral arteriovenous malformations[END_REF][START_REF] Pollock | Factors associated with successful arteriovenous malformations radiosurgery[END_REF][START_REF] Colombo | Three-dimensional angiography for radiosurgical treatment planning for arteriovenous malformations[END_REF][START_REF] Zhang | Clinical significance of 3D reconstruction of arteriovenous malformation using digital subtraction angiography and its modification with CT information in stereotactic radiosurgery[END_REF].

Whereas in the case of tumours this volume is obtained from computed tomography (CT) or magnetic resonance (MRI) imaging, digital subtracted angiography (DSA) is still the gold standard for determining the target volume in the case of BAVMs [START_REF] Zhang | Clinical significance of 3D reconstruction of arteriovenous malformation using digital subtraction angiography and its modification with CT information in stereotactic radiosurgery[END_REF][START_REF] Guo | Application of MR in stereotactic radiosurgery[END_REF][START_REF] Gauvrit | Three-dimensional dynamic magnetic resonance angiography for the evaluation of radiosurgically treated cerebral arteriovenous malformations[END_REF]. The AVM nidus consists of a complex tangle of abnormal vessels with arteriovenous shunting. Only the detailed analysis of different times of opacification of an angiographic series allows us to precisely determine the contours of this nidus and to differentiate it from the related arteries and the draining veins. The spatial and especially temporal resolutions of DSA are better than those of CT and MRI and thus allow for a more precise target definition.

Using DSA, the target volume is reconstructed from the contours drawn independently on two or more perpendicular or oblique views (see figure 1), in stereotactic conditions. Due to the nature of DSA, where the image is projected onto two dimensions, the threedimensional (3D) target volume cannot be reconstructed with complete accuracy [START_REF] Laurentini | How far 3D shapes can be understood from 2D silhouettes[END_REF]. This limitation of DSA has long been acknowledged, especially in the case of complexly shaped BAVMs [START_REF] Bova | Stereotactic angiography: an inadequate database for radiosurgery?[END_REF]. Methods have been developed that aim at better estimating the nidus shape from two perpendicular angiographic projections [START_REF] Söderman | Volume measurement of cerebral arteriovenous malformations from angiography[END_REF][START_REF] Coste | 3D reconstruction of the encapsulating contour of arteriovenous malformations for radiosurgery using digital subtraction angiography[END_REF]. These techniques basically assume the target presents an ellipse shape in every axial cut plane. Though they were shown to improve the accuracy of the target volume computation (in the sense of containment), their ability to provide an accurate target shape has not been demonstrated.

The aim of this paper is twofold: i. In section 2.1, we propose a method to validate AVM targets determined from a 3D imaging modality. ii. In section 2.2, we describe a clinical pilot study investigating the use of ("3D") rotational angiography (3DRA) to delineate the AVM radiosurgical target. A preliminary account of these findings has been presented in [START_REF] Berger | A methodology for validating a new imaging modality with respect to a gold standard imagery: Example of the use of 3DRA and MRI for AVM delineation[END_REF].

Material and methods

Statistical methods

Basic assumptions and outline of the method

The general method for validating a given 3D modality in the framework of RS was designed assuming the following experimental conditions. A set of patients was involved in the experiment. For each patient, a set of DSA views and one 3D image (e.g. MRI, CT or 3DRA) were available. The DSA views were corrected for geometric distorsions and the 3D image was spatially registered with the DSA views.

Two segmentation processes were to be compared: i. A test process, applied on the 3D images. This process could be either manual, with a new delineation protocol, or automatic, thanks to a new image processing algorithm. ii. A reference process, applied on the 2D images. In our case, this process was the manual contouring of the AVMs in the DSA views by experts. Each process generates 3D targets. Consequently, for each patient, a set of 3D targets was to be compared to a set of reference 3D targets by experts.

For each patient, all the test 3D targets and all the expert 3D targets were reprojected onto all DSA views. A set of test 2D curves together with a set of expert 2D curves were thus provided for each patient and each DSA view.

The general strategy of our method consisted in first designing a statistical test for a given test 2D curve to be accepted as valid, given a reference curve distribution. An acceptance rate was then computed for the test curves and compared to an expected acceptance rate estimated from the reference data.

Building a statistical shape model

A statistical shape model was built as introduced by Cootes [17]. This mathematical tool makes it possible to learn what variations in shapes are plausible, based on a set of expert reference curves. One model was built per DSA view, by analyzing shape variations over the training set of M 2D expert curves { C i } 1≤i≤M

Cootes' method [START_REF] Cootes | Statistical grey-level models for object location and identification[END_REF] requires that all curves are resampled on the same number N of points so that they are set into a one-to-one correspondence. This matching step must be done carefully for the shape model to make sense [START_REF] Davies | A minimum description length approach to statistical shape modeling[END_REF]. For our particular clinical application, dense point-to-point correspondences between curves were established using a semi-interactive interface (see figure 2). Segments were drawn by the user and linear interpolation was performed between these cuts. In practice, cuts were chosen so as to delimit areas where the expert curves were in good agreement from areas which presented high variability.

Each curve C can therefore be represented as a 2N dimensional vector, concatenating the x and y coordinates of the N sample points. The average curve is defined as the average vector:

 C= 1 M ∑ i C i (1)
with the covariance of the curve data:

Σ = 1 M ∑ i  C i - C  C i - C  t (2)
Still following Cootes [START_REF] Cootes | Statistical grey-level models for object location and identification[END_REF], a linear eigen-model was built by applying principal component analysis (PCA). The eigenvectors of Σ corresponding to the t largest eigenvalues represent the t principal modes of variation. Let P be the 2N×t matrix gathering these eigenvectors. A particular curve C is then approximated by:

C ' ≈  CPb
(3) where b=P t C-̄C  . The squared error of this approximation is:

r 2 C =∥C-C ' ∥ 2 =C- C  t C- C -b t b
(4)

Acceptance test for one curve

Once the model had been built, an acceptance test had to be designed to decide how a test curve matched the model. A first measure of the quality of fit is classically deduced from the PCA:

F m C =∑ k=1 t b k 2 λ k (5)
where λ k is the k -th largest eigenvalue in the PCA and b k the k -th component of b .

This criterion compares, for each mode k , the observed variation b k 2 on curve C to the expected variation λ k . This value only measures the distance by which the shape deviates from the mean along the modes of the model. It does not consider the distance r 2  C  of the test curve from its linear approximation (see equation ( 4)). Ideally, this distance should be small, meaning that the linear model given in equation ( 3) is able to closely model any curve bound to be an acceptable AVM contour.

In order to define what "small" means, an expected residual distance was evaluated for each linear model using leave-one-out experiments: for each expert curve C i , a reduced PCA model was built from all the expert curves except C i . The residual of this reduced PCA is denoted r i 2 . The expected residual was given by:

r e 2 = 1 M ∑ i=1 M r i 2  C i  (6)
Besides, the residual observed on a new curve C in this experiment was:

r o 2 = 1 M ∑ i=1 M r i 2 C  (7)
The compatibility of the new curve C was given by:

F r C = r o 2 r e 2 = ∑ i=1 M r i 2 C  ∑ i=1 M r i 2  C i  (8)
The entire measure of fit was then defined as the sum

F C =F m C F r C  (9)

Study of the distribution of F

F theoretically obeys a chi-square law. The actual distribution could be estimated through leave-one-out experiments on the expert curves as follows. For each expert curve C i , a reduced set of reference curves was considered, gathering all expert curves but C i . A measure of fit F i (see equation ( 9)) was then defined as above for the curve C i , based on this reduced set. Testing of all expert curves in all DSA views provided an estimate of what the actual distribution of the fit measure was amongst experts.

Global acceptance rate for all test curves

Thanks to the above chi-square test, an accept/reject decision could be taken for each DSA view. Consideration of all the DSA views led to a percentage of acceptance of test curves. Unless this percentage was convincingly high, the question arose of the significance of the computed acceptance rate. An expected acceptance rate could be estimated through leave-one-out experiments on the expert curves as above. An accept/reject decision was taken for each expert curve C i , based on F i . Consideration of all expert curves provided an estimate of what the "natural" acceptance rate was amongst experts, setting a reference for the acceptance rate of new contours.

Application to 3DRA

The method formally described above was applied to a preliminary investigation: 3DRA was tested to delineate complex BAVM target shapes.

Clinical material

Eleven patients with BAVM were consecutively treated using embolization and RS at the University Hospital of Nancy (France) between June 2001 and January 2002. According to the complexity of the AVM shape, between 2 to 5 (2D) DSA views were considered for delineation (antero-posterior, lateral and oblique views). Overall, 37 DSA series were available for this study. In addition, one 3DRA volume was available for each patient. These images were acquired in stereotactic conditions, allowing for DSA/3DRA image registration.

Reference data

Three experts in interventional neuroradiology from the University Hospital of Nancy (France) and Karolinska Institute (Stockholm, Sweden) were involved in this study. They manually delineated twice at a one month interval the AVMs in the 37 DSA images. The delineation was done with the software used in clinical practice (Advantage Windows, GE Healthcare). The full DSA series was available so that hemodynamics could be leveraged to analyze the angioarchitecture (feeding arteries, nidus and draining veins) and better assess the actual nidus contours. Each time, the AVM nidus was reconstructed, based on the expert delineation in all the available views. This reference 3D shape was then reprojected on each available view. Thus, 6 contours of reference were available on each DSA image (see figure 3).

Test data

For each patient Maximum Intensity Projection (MIP) images of the 3DRA were generated with the same viewing parameters as used for the acquisition of DSA images (see figure 4). The same three experts were asked to delineate the AVMs in these MIP views, using the same delineation software as for DSA images. In the absence of hemodynamics information in 3DRA, the 3DRA volume was available to display using volume rendering software (Volume Viewer, GE Healthcare), in a view to analyze the angioarchitecture using full 3D information. Again, each AVM nidus was reconstructed using cones intersection and this new 3D shape reprojected onto the available views. These contours were then tested using the statistical model built from the six expert delineations.

AVM classification

Each expert was asked to label each AVM with A or B according to the difficulty encountered when delineating the nidus on DSA series: label A was given when no major difficulty was faced, whereas label B was affixed when the delineation required a very thorough analysis with difficult choices to make. Along with a B labelling, the experts were asked to list their reasons underlying their choice. Figure 4 shows an example of an Alabelled AVM, and a B-labelled AVM can be seen in figure 5.

Berger 7

Results

Test parameters

For each DSA view, a shape model was built based on the 6 reference contours. The number t of modes considered to compute F m was chosen to represent 95% of the variation in the reference contours.

In order to study the actual distribution for F , the experiment described in section 2.1.4 was performed. The training set enabled us to test 222 expert curves ( 3×2 expert curves on 37 DSA series). Each measure of fit was associated to a number of modes for the underlying model: 163 models had 3 modes, 53 had 2 modes, and the last 6 models had 4 modes. The distribution of the F i 's based on 3 modes was therefore studied with the hypothesis that either 3 or 4 degrees of freedom should be considered for the actual chisquare distribution (number of modes vs number of summed terms in equation ( 9)). As shown in figure 6, the latter hypothesis is the one that was eventually retained.

Therefore, t+ 1 degrees of freedom and a p-value of 0.05 were considered in the experiment.

AVM classification

Each AVM was classified within type A or B, based on the consensus among the 3 experts: 5 AVMs were labelled as A (17 DSA series) and 6 as B (20 DSA series). Most noted difficulties were the lack of clarity of the nidus due to diffuse compartments, neoangiogenesis or previous embolization and the lack of a conspicuous nidus contour because vessels are superimposed.

Global acceptance rate

Results about the acceptance rates are summarized in table 1.

DSA acceptance rates call for two comments. First, in the case where the fit measure F would faithfully follow a chi-square distribution, a p-value of 0.05 would let us expect a reference acceptance rate of about 95%. It is not the case here, as is also apparent in figure 6: the tail of the observed distribution shows outliers. Second, the better acceptance rate for label A AVMs than for label B AVMs is consistent with the estimation by the experts of an easier, and therefore less variable, delineation of label A AVMs.

3DRA acceptance rates show that global AVM delineation cannot be performed on 3DRA alone (table 1, column for labels A&B) but should be considered as a complement to DSA for rather easy cases (table 1, column for A-labelled AVMs).

Per patient acceptance rate

A patient by patient analysis (table 2) reveals a more interesting situation. DSA acceptance rates demonstrate that our chi-square test enabled us to correctly label the delineation difficulty for all patients but A3 and B4: A-labelled AVMs correspond to a DSA acceptance rate above 85%.

The situation for 3DRA is somewhat different. 3DRA acceptance rates are similar to the reference DSA rate for 5/11 patients (A2, A3, A4, B2 and B3). The 3DRA rate is even higher than the DSA rate for 3 patients (A4, B2 and B3). In contrast, this rate is drastically low for 3 patients (B1, B5 and B6). These cases will be discussed in the next section.

Discussion

DSA as a reference

DSA still remains the benchmark imaging modality to delineate BAVM for RS. Besides having unequalled spatial resolution, DSA is the only modality that provides sufficient temporal resolution to study the hemodynamics. This point is key in analyzing BAVM: the nidus must be delineated at an early phase of the arterial time, when the nidus is completely opacified, just at the beginning of the venous drainage. Capturing this instant in time requires an image acquisition rate of at least 2 images per second. If taken too early, the DSA image might not display the whole nidus. If taken too late, superimposed venous drainage would degrade the nidus visibility.

However, it has long been acknowledged that DSA has severe limitations. Many views are required to delineate a 3D target but it is often difficult to find enough views where the nidus is clearly visible. Superimposed images of adjacent vessels, uncertain contours due to neoangiogenesis or previous embolization, or inhomogeneous opacification are common. As a consequence, the epipolar consistency of the contours between views is a condition that is seldom fulfilled and is not easy to control [START_REF] Anxionnat | Intra-and inter-observer variability in the angiographic delineation of brain arterio-venous malformations (AVMs)[END_REF], in particular when oblique views are added (see figure 7).

Further, even when many views are available, the methods for reconstructing AVM from contours drawn in DSA are unable to deal with all types of concavities in the shape [START_REF] Zhang | Clinical significance of 3D reconstruction of arteriovenous malformation using digital subtraction angiography and its modification with CT information in stereotactic radiosurgery[END_REF] (see figure 1). More generally speaking, an infinite number of 3D shapes fit a given set of any number of projections. For example, in figure 1 (right), both the bold line convex shape and the concave shape have the exact same reprojection contours in all 5 views.

Interest of 3D

3D imaging thus seems essential for the accurate determination of BAVM target shape. CT and/or MR alone do not seem to be reliable enough to differentiate the AVM nidus from the other adjacent vascular structures even when injected volume acquisitions are used [START_REF] Zhang | Clinical significance of 3D reconstruction of arteriovenous malformation using digital subtraction angiography and its modification with CT information in stereotactic radiosurgery[END_REF][START_REF] Aoyama | Comparison of imaging modalities for the accurate delineation of arteriovenous malformation, with reference to stereotactic radiosurgery[END_REF]. Nevertheless these 3D imaging modalities are used by many teams to refine the initial approximate volume reconstructed from DSA [START_REF] Flickinger | A multi-institutional analysis of complication outcomes after arteriovenous malformation radiosurgery[END_REF][START_REF] Zhang | Clinical significance of 3D reconstruction of arteriovenous malformation using digital subtraction angiography and its modification with CT information in stereotactic radiosurgery[END_REF][START_REF] Bednarz | Combining stereotactic angiography and 3D time-of-flight magnetic resonance angiography in treatment planning for arteriovenous malformation radiosurgery[END_REF].

Three-dimensional rotational angiography (3DRA) is a novel imaging technique [START_REF] Wiesent | Enhanced 3-D reconstruction algorithm for C-arm systems suitable for interventional procedures[END_REF] offering a selective vascular 3D analysis with a high spatial resolution (0.25 mm). 3DRA provides very useful 3D information that has been clinically validated for endovascular treatment of intracranial aneurysms [START_REF] Anxionnat | Clinical value of three-dimensional digital subtraction angiography (3D DSA) in the therapeutic decision and endovascular treatment of intracranial aneurysms[END_REF]. It can also be used for a better understanding of the nidus shape as shown in a preliminary study [START_REF] Colombo | Three-dimensional angiography for radiosurgical treatment planning for arteriovenous malformations[END_REF] but this technique has not yet been validated for use on AVM.

All 3D modalities lack the temporal resolution that DSA provides. As a result, fast flow compartments or diffuse parts of the nidus might not be visible in the 3D image. On the other hand, vessels and nidus never become superimposed in 3D. In the framework of BAVM radiosurgery, validating a 3D modality comes down to validating the target shape that can be determined from the 3D image alone. The clinical outcome must definitely be investigated. However, such a study will take a long time since the incidence of BAVM is very low and the variability of cases is high [START_REF] Ogilvy | Recommendations for the management of arteriovenous malformations: a statement for healthcare professionals from a special writing group of the stroke council[END_REF][START_REF] Al-Shahi | A systematic review of the frequency and prognosis of arteriovenous malformations of the brain in adults[END_REF]. In addition, only underestimation of target volume can be proven by lack of treatment effect. Overestimation of target volume may actually show up as an apparent increase in treatment efficacy, since the likelihood of obliteration will increase with the higher dose of radiation to the periphery of the BAVM. The concomitant increased risk for undue radiation effects will not be appreciated since the target definition was a priori correct.

A first validation step can be taken by comparing the target delineated in the 3D modality to the target reconstructed from reference standard DSA. Various criteria were proposed to that aim. The volume (expressed in cc) [START_REF] Söderman | Volume measurement of cerebral arteriovenous malformations from angiography[END_REF][START_REF] Coste | 3D reconstruction of the encapsulating contour of arteriovenous malformations for radiosurgery using digital subtraction angiography[END_REF][START_REF] Aoyama | Comparison of imaging modalities for the accurate delineation of arteriovenous malformation, with reference to stereotactic radiosurgery[END_REF][START_REF] Buis | Stereotactic radiosurgery for brain AVMs: Role of interobserver variation in target definition on digital subtraction angiography[END_REF], the position of the center of mass [START_REF] Zhang | Clinical significance of 3D reconstruction of arteriovenous malformation using digital subtraction angiography and its modification with CT information in stereotactic radiosurgery[END_REF][START_REF] Aoyama | Comparison of imaging modalities for the accurate delineation of arteriovenous malformation, with reference to stereotactic radiosurgery[END_REF][START_REF] Buis | Stereotactic radiosurgery for brain AVMs: Role of interobserver variation in target definition on digital subtraction angiography[END_REF][START_REF] Tanaka | Initial experience with helical CT and 3D reconstruction in therapeutic planning of cerebral AVMs: Comparison with 3D timeof-flight MRA and digital subtraction angiography[END_REF][START_REF] Yamamoto | Differences in target outline delineation from CT scans of brain tumours using different methods and different observers[END_REF], or the largest or mean dimension [START_REF] Colombo | Three-dimensional angiography for radiosurgical treatment planning for arteriovenous malformations[END_REF][START_REF] Aoyama | Comparison of imaging modalities for the accurate delineation of arteriovenous malformation, with reference to stereotactic radiosurgery[END_REF] of the target will provide valuable information, but cannot ensure that it will be noticed when reconstructed targets have the wrong shape. More complex criteria were proposed to better take into account the target shape: ratio of the intersecting volume over the encompassing volume [START_REF] Buis | Stereotactic radiosurgery for brain AVMs: Role of interobserver variation in target definition on digital subtraction angiography[END_REF][START_REF] Yamamoto | Differences in target outline delineation from CT scans of brain tumours using different methods and different observers[END_REF], overlapping contours [START_REF] Zhang | Clinical significance of 3D reconstruction of arteriovenous malformation using digital subtraction angiography and its modification with CT information in stereotactic radiosurgery[END_REF], largest dimensions in several directions [START_REF] Tanaka | Initial experience with helical CT and 3D reconstruction in therapeutic planning of cerebral AVMs: Comparison with 3D timeof-flight MRA and digital subtraction angiography[END_REF]. More recently, surface to surface distance computation, associated to local analysis along the target surface was proposed for lung cancer delineation [START_REF] Steenbakkers | Reduction of observer variation using matched CT-PET for lung cancer delineation: A three-dimensional analysis[END_REF].

Such studies failed to provide a definite answer about the actual use of the investigated 3D modality for AVM radiosurgery. Practical and clinical arguments were called upon in the end for recommending it as a complement to DSA. Whereas the above criteria tend to a better analysis of the target shape, the characterization of the correct shapes they provide is not local enough to discriminate wrong shapes. Furthermore, the above studies faced a fundamental difficulty in taking the target reconstructed in 3D from DSA as a reference: though great care was taken in using complex and state-of-the-art reconstruction methods, such a reference should still be acknowledged as possibly false and subject to high intra-and inter-observer variability [START_REF] Al-Shahi | Observer agreement in the angiographic assessment of arteriovenous malformations of the brain[END_REF].

Paper contributions

A new validation methodology

Evaluating a segmentation with respect to expert delineations has long been recognized as a difficult problem [START_REF] Chalana | A methodology for evaluation of boundary detection algorithms on medical images[END_REF][START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): An algorithm for the validation of image segmentation[END_REF]. Various methods have been proposed for comparing given segmentation with those of a group of experts. The simplest approaches, in 2D, use metrics on curve distances [START_REF] Williams | Comparing the joint agreement of several raters with another rater[END_REF]. Other techniques compute the agreement of the curve under consideration with the joint agreement of the experts, provided that a convenient metric on curves is available [START_REF] Chalana | A methodology for evaluation of boundary detection algorithms on medical images[END_REF]. More sophisticated methods try to define a statistical model of the shape [START_REF] Cootes | Statistical grey-level models for object location and identification[END_REF] or a probabilistic framework allowing the definition of the true segmentation and a measure of performance to be associated to each segmentation [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): An algorithm for the validation of image segmentation[END_REF].

In our application, manual delineations suffer from observer bias and large inter-and intra-observer variability (figure 5, left): areas where the consensus is good alternate with areas where the expert delineations are highly variable, preventing us from defining a reliable metric between curves and from using the methods described in [START_REF] Chalana | A methodology for evaluation of boundary detection algorithms on medical images[END_REF]. The STAPLE methodology [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): An algorithm for the validation of image segmentation[END_REF] aims at estimating the true segmentation from the expert delineations. This methodology is well suited to region segmentation but needs to incorporate a model describing spatially correlated structures. Such assumptions are clearly not fulfilled in our case: due to large image inhomogeneities in the AVM appearance, AVM detection cannot be formulated as a region segmentation problem. Moreover, due to the large variability of the AVM shape, it is impossible to build a model that provides prior information on the expected shape.

The first contribution of this paper is a methodology for validating a 3D imaging modality for BAVM delineation. The strong points of this methodology are the following.

The argument that the comparison should be made in 2D whenever DSA is considered as a reference. Indeed, 3D targets reconstructed from delineations in DSA are known not to be trustworthy everywhere: only the projections of these 3D shapes back onto the DSA views (i.e. 2D curves) can actually be considered as reference data. As a consequence, making the comparison in 3D is not valid.

The contour curves to serve as reference should be consistent in 3D. It occurred that some medical experts chose to draw a larger contour for the nidus in areas where the image was not clear, thus making use of the fact that a different DSA view could bring a correction to the current view. As a consequence, the original DSA contours were not to be considered as reference data but rather the reprojection of the reconstructed 3D target. Figure 3 demonstrates how this reduced the variability in the delineation.

The design of a statistical acceptance test for a given test curve, based on a PCA of the expert curves. This statistical test relies on the computation of a measure of fit F as the sum of two terms (see equation ( 9)). The first term F m controls variations along the modes of the model, that is, along the parts of the nidus contour where the experts disagree, binding the distance from the mean curve to remain within the observed limits. The second term F r penalizes any variation in ways not seen in the training set, that is, where the experts agree [START_REF] Cootes | Statistical grey-level models for object location and identification[END_REF]. As a result, this test makes it possible to locally adapt the decision of acceptance/rejection.Theoretically, F should follow a chi-square law with t degrees of freedom. In practice, the global chi-square shape is usually a valid hypothesis but the number of degrees of freedom must be adapted to the training set [START_REF] Cootes | Statistical grey-level models for object location and identification[END_REF]. Leave-oneout experiments showed that t+ 1 degrees of freedom should be considered for the chisquare distribution modeling F in the case of BAVM contours.

In this paper, the clinical experiment focused on 3DRA but this methodology can be applied to any 3D modality.

Preliminary investigation of 3DRA

The second contribution is clinical. A database of AVM cases has been constructed and the delineation process has been carefully studied. Our experiments showed that the variability between experts can be very large. Physicians should be more aware of the importance of this variability and protocols for delineation should be standardized. Moreover, computed acceptance rates showed that BAVM delineation could not be performed on 3DRA alone. However, this modality is promising and should be considered as a complement to DSA.

The per-patient analysis (see table 2) gave example cases of this complementarity. Cases B2 and B3 were considered as difficult ones, but present with a better acceptance rate using 3DRA than DSA. These were AVMs with a rather large nidus, strongly remodeled by previous embolization but with limited venous drainage. Conversely, 3DRA was unable to faithfully depict the nidus of 3 B-labelled patients (B1, B5 and B6). The corresponding AVMs had in common to present with an important and fast venous drainage that superimposed onto the nidus in DSA views, and was reconstructed in 3DRA volumes. In other words, 3DRA is a better aid to understanding complex 3D shapes, such as remodeled nidii, but lacks the temporal resolution of DSA, which is necessary to distinguish the draining veins from the nidus.

However, this last point highlights one limitation of this study: the views chosen to delineate the nidus in 3DRA were the same as in DSA since our aim was to demonstrate the potential of our methodology. In order to fully demonstrate the capacity of 3DRA to display the nidus with all the required conspicuity, and in particular to distinguish the nidus from the draining veins, the experts should be allowed to take better advantage of 3DRA by choosing carefully adapted views, even, and especially, views that are not allowed by DSA, such as axial (or top) views. We believe that the third spatial dimension of 3DRA can compensate for its lack of temporal resolution. Moreover, the contrast medium injection protocol for 3DRA should be improved for cases with a large draining vein.

Future works

Several improvements and extensions can be made to our approach. The database must be enriched, especially with label A malformations. It would also be convenient if the construction of the statistical model of the shape could be facilitated. Short-term improvements consist in allowing curves instead of segments for cuts in our software interface used to match contours (see dark segments on figure 2). More robust and automated matching could be used but it either requires to have similar curves [START_REF] Tagare | Shape-based nonrigid correspondence with application to heart motion analysis[END_REF] or to define expected properties of the dense curve matching [START_REF] Davies | A minimum description length approach to statistical shape modeling[END_REF]. None of these conditions are met in our case.

Multimodality methods are under investigation in order to improve the current delineation process. One research direction is to use 3DRA to build virtual 2D views. Hence, the physician could potentially use an arbitrary number of views to delineate the AVM without further dose exposure to the patient. The accuracy of the target shape will be improved as the precision of the shape recovered with cone intersections is tightly related to the number of views that can be considered.

Outlining a 3D shape is classically performed by manually drawing contours in 2D cutplanes through the pathology. With AVMs, this is a very difficult task (see figure 8). In this paper, the 3DRA target volume was built by cone intersection from contours drawn in MIP views. Concavities in the target cannot be fully recovered with this method. Among other advantages, 3D modalities, and in particular 3DRA, should make it easier to segment the AVM nidus than when DSA is used. A further step to improve the target delineation is to use quasi-automated algorithms to reduce variability. We are currently working on evaluating active surfaces to that aim.

Figures and tables

Fig

  Fig. 1: AVM reconstruction using conventional angiography. Left: setup with 2 orthogonal views. Right: concavities cannot be recovered even with many views.

  Fig. 1: AVM reconstruction using conventional angiography. Left: setup with 2 orthogonal views. Right: concavities cannot be recovered even with many views.

Fig. 2 :

 2 Fig. 2: Semi-interactive interface to match curves: the segments in black were indicated by the user to define homologous sections of curves; the segments in gray show the point-topoint correspondence established between two curves within the set.

Fig. 3 :

 3 Fig. 3: (From left to right) 1) The set of 6 original expert delineations. 2) reconstructed surface in the same orientation as the DSA view 3) The new database using the reprojection of the 3D shape recovered from expert delineations for the same view. 4) Reprojection of the recovered 3D shape on a reference image in semi-transparent lightgray along with the expert delineation. Data for patient A2 in table 2.

Fig. 4 :

 4 Fig. 4: Example of A-labelled AVM. From left to right in AP orientation: DSA image; DSA image with all original expert delineations; 3DRA MIP image with all original expert delineations; 3DRA MIP. Little vessel superimposition and a limited venous drainage for a nidus with moderate flow speed. A high quality 3DRA reconstruction with no artefact displays an apparently even more conspicuous nidus. Data for patient A2 in table 2.

Fig. 5 :

 5 Fig.5: Example of B-labelled AVM. From left to right in AP orientation: DSA image; DSA image with all original expert delineations; 3DRA MIP image with all original expert delineations; 3DRA MIP. This large nidus appears unclear on most of its contour due to artery superimposition and large previously embolized areas (right and bottom). The venous drainage hampers the analysis of the upper part of the nidus in 3DRA. Experts excluded this part in 3DRA whereas it was always included in DSA. Moreover, the 3D reconstruction of the nidus presents with a poor quality especially in the previously embolized areas. Data for patient B3 in table 2.

Fig. 6 :

 6 Fig. 6: Distribution of the measures of fit for the expert curves in DSA views corresponding to models with 3 modes: the data histogram fits the chi-square distribution with 4 degrees of freedom.

Fig. 7 :

 7 Fig. 7: Contour consistency. Top row: DSA views with one original expert manual delineation (curve in white) and the reconstructed target (semi-transparent light grey) in antero-posterior, lateral and oblique orientations (from left to right). The reconstructed target in AP orientation is shown on the right. Second row: the DSA views with no graphics superimposed for reference. The reconstructed target in lateral orientation is shown on the right. Even though every contour seems correct independently on each DSA view, the large difference in the AP view between the initial contour and the reprojected target demonstrates their inconsistency in 3D. Data for patient A2 in table 2.

Fig. 8 :

 8 Fig.8: Difficulties encountered when reconstructing an AVM 3D shape from manual outlines in a stack of axial cut planes. (Left) A 3DRA axial plane through an AVM: the different vascular structures cannot be easily distinguished. (Middle) Two consecutive axial planes in T1 Gd-enhanced MRI together with the manually drawn AVM (Right) Reconstruction of an AVM by stacking up these manual outlines: the smoothness of the 3D shape is not ensured in the longitudinal direction.

Table 1

 1

		global acceptance rate	acceptance for	acceptance for
		(labels A & B)	label A AVMs	label B AVMs
	3DRA	69%	82%	58%
	DSA	84%	87%	81%

Table 1 :

 1 Acceptance rate of 3DRA contours (second line) compared to the acceptance rate of DSA (third line) computed as explained in section 2.1.5. DSA rates were established on 222 contours and 3DRA rates on 111 contours.

					Table 2
	Patient	A1 A2 A3 A4 A5 B1 B2 B3 B4 B5 B6
	Nb of DSA series	3	3	4	3	4	4	3	4	3	3	3
	DSA acceptance rate 89% 94% 79% 89% 88% 83% 72% 79% 89% 83% 83%
	3DRA acceptance rate 78% 89% 75% 100% 75% 17% 89% 100% 67% 33% 44%

Table 2 :

 2 Per patient acceptance rates: DSA rates are computed on 18 or 24 curves (6 curves per DSA series) and 3DRA rates on 9 or 12 curves (3 curves per 3DRA MIP view). The DSA data consisted of a total of 37 series; Patients A1 to A5 are A-labelled and B1 to B6 are B-labelled. Note the AP data for patient A2 are presented in figure4and AP data for patient B1 are presented in figure5.