Claude Le Bris
email: lebris@cermics.enpc.fr

Mazyar Mirrahimi
email: mazyar.mirrahimi@ensmp.fr

Herschel Rabitz
email: hrabitz@princeton.edu

Gabriel Turinici
email: gabriel.turinici@inria.fr

Hamiltonian identification for quantum systems : well posedness and numerical approaches

Keywords: 2000 Mathematics Subject Classification. 93B30, 65K10. December 18, 2005 Inverse problem, Quantum systems, Hamiltonian Identification, Optimal Identification * CLB, MM

This paper considers the inversion problem related to the manipulation of quantum systems using laser-matter interactions. The focus is on the identification of the field free Hamiltonian and/or the dipole moment of a quantum system. The evolution of the system is given by the Schrödinger equation. The available data are observations as a function of time corresponding to dynamics generated by electric fields. The well-posedness of the problem is proved, mainly focusing on the uniqueness of the solution. A numerical approach is also introduced with an illustration of its efficiency on a test problem.

Introduction

The ability of coherent light to manipulate molecular systems at the quantum scale has been demonstrated both theoretically and experimentally. Various methods have been shown to be relevant for different systems and settings, ranging from single atoms to large polyatomic molecules in the liquid phase [3-5, 7, 14, 15, 15-21, 23, 28, 29]. Many of these procedures are based on the possibility to perform a large number of experiments in a very small time frame. Thus, the output provided by these experiments can be used to correct the process and to identify more satisfactory control fields [10,[START_REF] Li | Optimal dynamic discrimination of similar molecules through quantum learning control[END_REF][START_REF] Phan | Learning control of quantum-mechanical systems by laboratory identification of effective input-output maps[END_REF].

The ability to rapidly generate a large amount of quantum dynamics data may also be used to extract more information about the possibly unknown parameters of the quantum system itself. For each test field (i.e., control), there is the possibility of performing many observations for deducing information about the system, and this process can often be carried out at a much faster rate than the associated numerical simulations of the dynamics. Moreover, the recent advances in laser technology provide the means for generating a very large class of test fields for such experiments. The inverse problem, called Hamiltonian identification has been formulated within an optimization framework in various settings [START_REF] Geremia | Optimal hamiltonian identification: The synthesis of quantum optimal control and quantum inversion[END_REF][START_REF] Kosut | Identification of quantum systems[END_REF][START_REF] Shenvi | Nonlinear kinetic parameter identification through map inversion[END_REF][START_REF] Tadi | Explicit method for parameter identification[END_REF]. However, few theoretical results are available to date that illuminate the relationship between the output measurements and the quality of the inversion [START_REF] Alis | On the inversion of quantum mechanical systems: Determining the amount and type of data for a unique solution[END_REF][START_REF] Zhu | Potential surfaces from the inversion of time dependent probability density data[END_REF]. To address these issues this paper introduces the theoretical grounding that formalizes the above problem. We then show the well-posedness of the inversion problem and also provide a numerical illustration for the case of a simple 4-level system. Consider a quantum system evolving according to the Schrödinger equation

i Ψ = (H 0 + V + ǫ(t) µ)Ψ, Ψ| t=0 = Ψ 0 , (1)
where H 0 is the kinetic energy operator, V the potential operator and µ the dipole moment operator coupling the system to a time-dependent external laser field ǫ(t). We assume this system is adequately described as finite dimensional, so that H 0 , V and µ are N × N Hermitian matrices with entries in C and Ψ(t) ∈ C N . Such a finite dimensional configuration typically represents an approximation of an infinite dimensional system. We denote by {φ i } N i=1 the eigenvectors of H 0 + V and {λ i } N i=1 the associated eigenvalues and throughout this paper H = H 0 + V is referred to as the Hamiltonian of the system in the absence of the laser field.

The goal is to determine information about H and µ from laboratory measurements on some observables of the system. Two different settings are considered for this identification problem depending on which parameters are to be identified and the nature of the information:

(S1): The Hamiltonian H is known and the goal is to identify the dipole moment µ. The so called populations along the eigenstates φ i , i.e. p i (t) = | φ i , Ψ(t) | 2 , i = 1, 2, ..., are measured for all instants t ≥ 0. This is performed with as many control amplitudes ǫ(t) as required. (S2): Neither the potential V nor the dipole moment µ are known and the goal is to identify them. Note that, by identifying H, we mean, throughout this article, identifying V , as H 0 is readily known. The eigenvalues of the Hamiltonian H = H 0 + V are also assumed to be known (this assumption is relevant in practice, see Remark 7). Here, we measure the populations p i along the states of a canonical basis

{e i } N i=1 : p i (t) = | e i , Ψ(t) | 2 , i = 1, 2, .
.. for all instants t > 0 and all control amplitudes ǫ(t). Two fundamental questions concerning the well-posedness of these problems arise: the existence and the uniqueness of the Hamiltonian, and/or the dipole moment, compatible with the given measurements. In addition, these questions may be formulated either in an infinite dimensional setting or in the finite dimensional setting, as above. It is to be noted that, in general, the internal Hamiltonian H is infinite dimensional, which gives by truncation the finite dimensional approximation that we analyze here. Besides, in practice high energy levels are rarely excited in experiments and thus the finite dimensional configuration seems to be relevant. Additional laboratory measurements could also considered beyond the populations. This set of measurements could contain the expectation value of a set of physical observables { Ôi } (here Ôi are Hermitian matrices with elements in C:

{ Ψ | Ôi Ψ | i = 1, 2, ...} where u | v = N i=1 u i v *
i is the Hermitian product). We only consider the uniqueness problem in the two settings mentioned above. This is the topic of Section 2. We do not address the existence question and we suppose that all measurements result from an existing physical system. In addition, the bilinear model ( 1) is assumed to be exact for this physical system (but V and µ may be unknown). Along with the questions of existence and uniqueness, the problem also entails finding an algorithm to extract as much information as possible from the measurements. Section 3 proposes a numerical algorithm similar to the approach in [START_REF] Geremia | Optimal hamiltonian identification: The synthesis of quantum optimal control and quantum inversion[END_REF] and tests its performance on a 4-level system.

To conclude the Introduction, we suppose throughout this article that the following assumptions hold true for the system under consideration:

(A1): Equation ( 1) is wavefunction controllable [START_REF] Ramakrishna | Controllability of molecular systems[END_REF]. (A2): The transitions of the Hamiltonian H are non-degenerate: λ i1 -λ j1 = λ i2 -λ j2 for (i 1 , j 1 ) = (i 2 , j 2 ) [START_REF] Turinici | Quantum wavefunction controllability[END_REF]. (A3): The diagonal part of the dipole moment µ, when written in the eigenbasis of the Hamiltonian H, is zero:

φ i | µ | φ i = 0, i = 1, .., N .
Some comments regarding the above assumptions are presented below.

The notion of wavefunction controllability of a quantum system, used in Assumption (A1), has been introduced in [START_REF] Ramakrishna | Controllability of molecular systems[END_REF], and a criterion of controllability involving the Lie algebra spanned by iH and iµ is given therein. Under an adequate assumption, this notion of controllability implies the existence of a time T 0 such that for any time T > T 0 and any pair of initial and target states on the unit sphere of the state space, there exists a piecewise constant control field ǫ(t) driving the system from the initial state to the target state in time T . For more details concerning the controllability of finite dimensional quantum systems, we refer to the Appendix A.

Assumption (A2) is relevant to local controllability [START_REF] Turinici | Quantum wavefunction controllability[END_REF], and it is a sufficient condition for global controllability for systems that are not decomposable into independent subsystems (i.e., the connectivity graph is connected in terms specified in [START_REF] Turinici | Quantum wavefunction controllability[END_REF]). Furthermore, related techniques may enable the use of weaker forms of (A2) (see the above cited works).

Assumption (A3) is necessary to avoid a systematic ambiguity in the identification process. Changing the dipole moment µ into µ + D, where D is a diagonal matrix (in the eigenbasis of H), does not change any of the observed populations. In this case, the state Ψ(t) is replaced by exp -iD t 0 ǫ(s)ds Ψ(t). Assumption (A3) can thus be made without loss of generality.

Remark 1. In taking laboratory data noise is a significant factor that has to be dealt with. In the formalization of the inverse problem above we do not explicitly consider a model for noise. Tolerance factors are introduced (see Section 3.4) which can be viewed as a numerical technique to help treat noise. Thus, this study should be viewed as a first step towards future works that fully take into account the noise.

Uniqueness

Setting (S1): extracting the dipole with the Hamiltonian known

To establish the well-posed nature of the problem it is necessary to prove that at most one dipole moment exists that is compatible with the measurements. Some degree of indeterminacy generally cannot be avoided if special symmetries arise in the system. Apart from this situation, uniqueness holds true, as stated by Theorem 1 below. The Hamiltonian H is known and becomes a diagonal matrix in terms of its eigenbasis with entries {λ i } N i=1 . The main result of this subsection is: Theorem 1. Suppose that there exist two dipole moments µ 1 and µ 2 , giving rise to two evolving states Ψ 1 and Ψ 2 respectively solving

i Ψ1 = (H + ǫ(t) µ 1 ) Ψ 1 , (2)
i Ψ2 = (H + ǫ(t) µ 2 ) Ψ 2 , (3)
that produce identical observations for all times t ∈ [0, ∞) and all fields ǫ(t):

| Ψ 1 (t) | φ i | 2 = | Ψ 2 (t) | φ i | 2 i = 1, .., N. (4)
Then under the assumptions (A1),(A2) and (A3), the two dipole moments are equal within some phase factors

{α i } N i=1 ⊂ R such that: ∀i, j = 1, .., N, (µ 1 ) ij = e i(αi-αj ) (µ 2 ) ij . (5)
The uncertainty (5) in the relative phases of the dipole matrix elements is due to the fact that any diagonal unitary transformation U : U = diag(e iα1 , e iα2 , .., e iα N ) leaves both the Hamiltonian H and the populations p i (t) = | Ψ(t) | φ i | 2 unchanged. We give an example below.

Example 1. Consider the following two systems given by the same Hamiltonian H:

H =   E 1 0 0 0 E 2 0 0 0 E 3   ,
but two different dipole moments µ -and µ + :

µ -=   0 -α 0 -α 0 β 0 β 0   µ + =   0 α 0 α 0 β 0 β 0   .
The two corresponding coupled systems of Schrödinger equations may be respectively written as follows

(Σ)    ı d dt Ψ 1 = E 1 Ψ 1 -α ǫ(t) Ψ 2 ı d dt Ψ 2 = E 2 Ψ 2 -α ǫ(t) Ψ 1 + β ǫ(t) Ψ 3 ı d dt Ψ 3 = E 3 Ψ 3 + β ǫ(t) Ψ 2
and

(Σ ′ )      ı d dt Ψ 1 = E 1 Ψ 1 + α ǫ(t) Ψ 2 ı d dt Ψ 2 = E 2 Ψ 2 + α ǫ(t) Ψ 1 + β ǫ(t) Ψ 3 ı d dt Ψ 3 = E 3 Ψ 3 + β ǫ(t) Ψ 2 Denote by Ψ = (Ψ 1 , Ψ 2 , Ψ 3 ) the solution of the system (Σ). Then Ψ = ( Ψ 1 , Ψ 2 , Ψ 3 ) = (Ψ 1 , -Ψ 2 , -Ψ 3 )
is the solution of the system (Σ ′ ). Thus, the wavefunctions Ψ and Ψ give rise to identical population observations for all control fields ǫ(t) at all instants t > 0 (4). So, the two dipole moments are solutions to the inversion problem.

Proof of Theorem 1:

Differentiating (4) with respect to time, we have (denoting by ℑ the imaginary part of a complex number):

ǫ(t) ℑ( µ 1 Ψ 1 (t) | φ i φ i | Ψ 1 (t) ) = ǫ(t) ℑ( µ 2 Ψ 2 (t) | φ i φ i | Ψ 2 (t) ). ( 6
)
From this equation we deduce

ℑ( µ 1 Ψ 1 (t) | φ i ) φ i | Ψ 1 (t) ) = ℑ( µ 2 Ψ 2 (t) | φ i φ i | Ψ 2 (t) ) ∀i = 1, 2, ... (7)
If for some particular time t, (7) cannot be deduced from (6) because ǫ(t) = 0, we can continuously extend the control field on (t, +∞) such that ǫ(t + s) = 0, s > 0 (without changing the solution of the Schrödinger equation on the interval [0, t]). As (6) holds true for any choice of ǫ(t) and as ǫ(t + s) = 0 for s > 0, equation ( 7) is satisfied for t + s with s > 0. Now by continuity and passing to the limit when s tends to zero, equation ( 7) is valid for t. An illustration of the procedure is presented in fig. 1. As ( 2) is assumed to be controllable, there exists a control ǫ(t) (defined on [0, T ]) that drives the system from Ψ 0 to

Ψ 1 (T ) = 1 √ N N j=1 φ j . (8)
Note that since the system is controllable any final state is accessible, and the form in ( 8) is chosen for its convenience in the following analysis. Using (4) for t = T , we have

Ψ 2 (T ) = 1 √ N N j=1 e iαj φ j , for some N -tuple (α 1 , ..., α N ) ∈ R N t t 2 1 ε Control field intensity 0 0 0 continuation at t continuation at t' Figure 1.
Illustration of the derivation of eqn ( 7) from ( 6). The field ǫ of ( 6) is assumed to be the null field ǫ ≡ 0 between t 1 and t 2 . In order to obtain Eqn. [START_REF] Chen | Competitive tracking of molecular objectives described by quantum mechanics[END_REF] we construct for each time (e.g., t and t ′ ) a non-trivial continuation. Evaluating [START_REF] Chen | Competitive tracking of molecular objectives described by quantum mechanics[END_REF] at time t = T and using (µ 1 ) ii = (µ 2 ) ii = 0 yields:

j =i ℑ( µ 1 φ j |φ i ) = j =i ℑ(e i(αj -αi) µ 2 φ j |φ i ).
Once the system in equation (2) (resp. equation ( 3)) is in the state1

√ N N j=1 φ j (resp. Ψ 2 (T ) = 1 √ N N j=1
e iαj φ j ), we choose the control term ǫ(T + s) to be zero for s > 0.

Using the same arguments as above, equation ( 7) still holds:

j =i ℑ(e iωij s µ 1 φ j |φ i ) = j =i ℑ(e iωij s e i(αj -αi) µ 2 φ j |φ i ).
where ω ij = λ j -λ i . Let us take θ 1 ij and θ 2 ij to be, respectively, the phases of µ 1 φ j | φ i and µ 2 φ j | φ i . 1 Thus, for all s > 0,

j =i sin(ω ij s + θ 1 ij )| µ 1 φ j | φ i | = j =i sin(ω ij s + α j -α i + θ 2 ij )| µ 2 φ j | φ i |. (9)
The non-degeneracy assumption (A2) then implies

| µ 1 φ j | φ i | = | µ 2 φ j | φ i |,
which means that the entries ij of the dipole moments µ 1 and µ 2 are equal up to a phase factor ∀i, j, there exists

θ ij such that (µ 1 ) ij = e iθij (µ 2 ) ij where θ ij = θ 2 ij -θ 1 ij . In addition, equation (9) yields sin(ω ij s + θ 1 ij ) = sin(ω ij s + θ 2 ij + α j -α i ) ∀i, j = 1, ..., N and ∀s > 0. ( 10
)
Equation (10) leads to

e iθij = e i(θ 2 ij -θ 1 ij ) = e i(αi-αj ) which concludes the proof.
Remark 2. As implied by the above proof, the dynamical time needed for the inversion is "heuristically" of the same order of magnitude as the shortest control time T 0 for the system. Note however that the proof will still work choosing a different Ψ 1 (t) in [START_REF] Cohen-Tannoudji | Mécanique Quantique, volume I& II[END_REF] provided it has nonzero populations in all levels. It was proved [27, Lemma 2], that under the hypothesis above the time T necessary to reach at least one such point can be made arbitrarily small. If no practical considerations are taken into account this would imply that the time required for inversion could be arbitrarily small.

Setting (S2): both the Hamiltonian and the dipole moment are unknown

This subsection proves the following main result: Theorem 2. Suppose that there exist two Hamiltonians H 1 and H 2 and two dipole moments µ 1 and µ 2 , giving rise to two evolving states Ψ 1 and Ψ 2 respectively solving

i Ψ1 = (H 1 + ǫ(t) µ 1 ) Ψ 1 , (11)
i Ψ2 = (H 2 + ǫ(t) µ 2 ) Ψ 2 , (12)
that produce identical observations for all times t ∈ [0, ∞) and all fields ǫ(t)

| Ψ 1 (t) | e i | 2 = | Ψ 2 (t) | e i | 2 i = 1, .., N, (13)
where {e i } N i=1 is the canonical basis of the state space C N . Assume (A1), (A2), (A3) for both Hamiltonians H 1 and H 2 (see however Remarks 3 and 4 below) and assume the following additional condition: (A4): the Hamiltonians H 1 and H 2 have the same eigenvalues. Assume also (see Remark 8 below), that (A5): there does not exist a subspace of dimension one or two spanned by the vectors {e i }, which remains invariant during the free evolution (ǫ ≡ 0) of the first system (H 1 and µ 1 ).

Then, there exists a set of phases {α i } N i=1 such that, for all 1 ≤ i, j ≤ N ,

(µ 1 ) ij = e i(αi-αj ) (µ 2 ) ij , (14)
(H 1 ) ij = e i(αi-αj ) (H 2 ) ij . (15)
A number of remarks are in order.

Remark 3. Assumptions (A1) and (A2) for, say, H 1 , along with Assumption (A4) and ( 13) for both, imply that (A1) also holds for H 2 .

Remark 4. Assumption (A2) for one of the two Hamiltonians, along with Assumption (A4), establishes (A2) for the other Hamiltonian.

Remark 5. As in Subsection 2.1, the uncertainty in the relative phases of the matrix elements is due to a diagonal unitary transformation U : U = diag(e iα1 , e iα2 , .., e iαN ) which leaves the measurement of populations unchanged. Remark 6. Note that, in contrast with the setting of Theorem 1, the observations here refer to populations along a canonical basis [START_REF] Kullback | On information and sufficiency[END_REF]. As the Hamiltonian is not known, the measurement of the populations in the eigenstates of H is not available a priori. However, even if the latter populations were available, their knowledge would still not allow for the identification of H and µ. Such a setting would then be similar to Section 2.1, and H and µ will only be identified up to a unitary transformation of the eigenbasis of H.

Remark 7. The Assumption (A4) is physically relevant, as spectroscopy allows for identifying the eigenvalues of the Hamiltonian, and thus discriminating between two systems that do not share the same ones. In fact spectroscopy only gives eigenvalue differences, not the absolute values. This overall unknown additive factor is not seen by the measurements and has no impact on the identification result.

Remark 8. Assumption (A5) is made here for simplicity. Theorem 3 relaxes this assumption and shows that the only necessary assumption, replacing (A5), is that the system be of dimension strictly greater than 2.

Remark 9. As will be clear from the proof below, one may consider any arbitrary fixed basis, independent of the system, instead of the canonical basis {e i } N i=1 .

Proof of Theorem 2:

Differentiating equation ( 13) with respect to time produces

ℑ( H 1 Ψ 1 | e i e i | Ψ 1 )+ǫ(t)ℑ( µ 1 Ψ 1 | e i e i | Ψ 1 ) = ℑ( H 2 Ψ 2 | e i e i | Ψ 2 ) + ǫ(t)ℑ( µ 2 Ψ 2 | e i e i | Ψ 2 ) (16)
Consider an eigenbasis (φ 1 i ) N i=1 of H 1 and (φ 2 i ) N i=1 an eigenbasis of H 2 . Each vector e j of the canonical basis is decomposed into these two eigenbases

e j = c 1 j,i φ 1 i = c 2 j,i φ 2 i ∀j ∈ 1, .., N .
Both matrices c 1 and c 2 are unitary.

The proof is divided into three steps. In the first step, we prove the equality |c 1 i,j | = |c 2 i,j | for i, j = 1, .., N , thus showing that the entries are equal within a phase factor. In the second step, we show that such phases satisfy a special relation (see [START_REF] Ramakrishna | Controllability of molecular systems[END_REF] below). This yields [START_REF] Levis | Selective bond dissociation and rearrangement with optimally tailored, strong-field laser pulses[END_REF]. The third step concludes with an analysis of the dipole moment and shows [START_REF] Bris | Towards efficient numerical approaches for quantum control[END_REF].

Step 1: As the physical system is assumed to be controllable, we may choose ǫ(t) such that Ψ 1 (T ) = e 1 and ǫ(T ) = 0. 2 Using (13), we have Ψ 2 (T ) = e iθ e 1 for some phase θ. Then let the system freely evolve after the time T with the control ǫ(T + s) = 0 for s > 0. Thus,

Ψ 1 (T + s) = c 1 1,i e -isλi φ 1 i , Ψ 2 (T + s) = c 2 1,i e -isλi+iθ φ 2 i . (17)
Equation ( 16) and ǫ(T + s) = 0 for s ≥ 0 leads to

ℑ( H 1 Ψ 1 (T + s) | e 1 e 1 | Ψ 1 (T + s) ) = ℑ( H 2 Ψ 2 (T + s) | e 1 e 1 | Ψ 2 (T + s) ) s ≥ 0. ( 18
)
2 When a time reversible system is controllable, it is always possible, for any two states Ψ 1 and Ψ 2 , to find a control ǫ 1,2 that steers Ψ 1 to Ψ 2 and is null at the final time: ǫ 1,2 (T ) = 0. For example, by propagating Ψ 2 backwards in time to some Ψ 3 at time Tη, η > 0, by free evolution, we may then steer Ψ 1 to Ψ 3 by an appropriate control field, and then let Ψ 3 freely evolve to Ψ 2 .

Additionally,

H 1 Ψ 1 (T + s) | e 1 e 1 | Ψ 1 (T + s) = λ i c 1 1,i e -iλis φ 1 i | c 1 1,i φ 1 i c 1 1,i φ 1 i | c 1 1,i e -iλis φ 1 i = i λ i |c 1 1,i | 2 e -iλis   j |c 1 1,j | 2 e iλj s   = i,j λ i |c 1 1,i | 2 |c 1 1,j | 2 e iωij s , (recall that ω ij = λ j -λ i ) = N -1 i=1 N j=i+1 (λ i e iωij s + λ j e -iωij s )|c 1 1,i | 2 |c 1 1,j | 2 + N i=1 λ i |c 1 1,i | 4 .
As the same computations hold true for H 2 and Ψ 2 , equation ( 18) leads to

N -1 i=1 N j=i+1 |c 1 1,i | 2 |c 1 1,j | 2 (ω ij sin(ω ij s)) = N -1 i=1 N j=i+1 |c 2 1,i | 2 |c 2 1,j | 2 (ω ij sin(ω ij s)).
Assumption (A2) then implies

|c 1 1,i | |c 1 1,j | = |c 2 1,i | |c 2 1,j | ∀i = j, i, j = 1, .., N. (19)
Consider now the first line of the matrix c 1 . There are three cases to be considered:

(1) there exists more that two non-zero entries;

(2) there are exactly two non-zero entries;

(3) there is exactly one non-zero coordinate. In case (1), we denote c 1 1,i , c 1 1,j and c 1 1,k as three non-zero entries. Using equation [START_REF] Mirrahimi | Reference trajectory tracking for locally designed coherent quantum controls[END_REF] for these terms produces

|c 1 1,i | |c 1 1,j | = |c 2 1,i | |c 2 1,j |, |c 1 1,i | |c 1 1,k | = |c 2 1,i | |c 2 1,k |, |c 1 1,j | |c 1 1,k | = |c 2 1,j | |c 2 1,k |, implying that |c 1 1,i | = |c 2 1,i |. (20)
One can argue an analogous relation for all other non-zero entries of the first line of the matrices.

In case (2), c 1 1,i and c 1 1,j are the two non-zero entries. In the i-th column of the matrix c 1 there is necessarily a non-zero coordinate c 1 k,i other than c 1 1,i . We have |c 1 1,i | < 1 and the norm of column i is one, since c 1 is a unitary matrix.

Consider now the k-th line. Again because the matrix c 1 is unitary, we have

c 1 1,i c 1 * k,i + c 1 1,j c 1 * k,j = 0.
Since the terms c 1 1,i , c 1 1,j and c 1 k,i are non-zero, the entry c 1 k,j must be also non-zero. So in the k-th line, we have at least two non-zero entries c 1 k,i and c 1 k,j . Using (A5), we know there are at least three non-zero entries in this line. Otherwise, the subspace of dimension 2 spanned by {e 1 , e k } would coincide with the eigenspace spanned by {φ 1 i , φ 1 j }. Arguing on line k, as we did above on line 1 when studying case (1), leads to

|c 1 k,l | = |c 2 k,l | ∀l = 1, .., N.
Thus, there exists some phases θ i and θ j such that:

c 2 k,i = e iθi c 1 k,i , c 2 k,j = e iθj c 1 k,j .
Using that both c 1 and c 2 are unitary matrices yields

c 1 1,i c 1 * k,i + c 1 1,j c 1 * k,j = 0 , c 2 1,i c 1 * k,i e -iθi + c 2 1,j c 1 * k,j e -iθj = 0, and |c 1 1,i | = |c 2 1,i | , |c 1 1,j | = |c 2 1,j |.
We finally show that case (3) never occurs. Denote by c 1 1,i the only non-zero element of the first line. With the matrix c 1 being unitary, we have |c 1 1,i | = 1, thus e 1 = e iθi φ 1 i for some phase θ i . But, Assumption (A5) of the theorem forbids this situation: the subspace of dimension one spanned by e 1 would coincide with the eigenspace spanned by φ 1 i . In both cases (1) and ( 2), we have obtained that all entries of the first line of c 1 and c 2 coincide in modulus

|c 1 1,i | = |c 2 1,i |,
for all indexes i. Arguing on each line of the matrices as we did on the first line, we obtain the existence of a set {θ ij |i, j = 1, 2, .., N }, such that

c 2 i,j = e iθij c 1 i,j . (21)
Step 2: We now show that there exists two N -tuples {α i } and {β j } such that, for all 1 ≤ i, j ≤ N ,

θ ij = α i + β j . (22)
To prove [START_REF] Ramakrishna | Controllability of molecular systems[END_REF], consider again equation ( 16) and assume that Ψ 1 (T ) = e 1 , which implies that Ψ 2 (T ) = e iθ e 1 .

Leaving the system evolving freely (with ǫ(T + s) = 0) after the time T , the relations in [START_REF] Maday | New formulations of monotonically convergent quantum control algorithms[END_REF] will still be valid. Equation ( 16) yields

ℑ( H 1 Ψ 1 (T + s) | e k e k | Ψ 1 (T + s) ) = ℑ( H 2 Ψ 2 (T + s) | e k e k | Ψ 2 (T + s) ) s ≥ 0. ( 23
)
The same computations as above yield

c 1 1,i c 1 * k,i c 1 * 1,j c 1 k,j = c 2 1,i c 2 * k,i c 2 * 1,j c 2 k,j ∀i, j = 1, .., N. (24)
This relation together with [START_REF] Rabitz | Perspective. shaped laser pulses as reagents[END_REF] implies that

θ 1i -θ ki + θ kj -θ 1j = 0 ∀i, j = 1, .., N.
More generally, using any e k instead of e 1 , we obtain

θ li -θ ki + θ kj -θ lj = 0 ∀i, j, k, l = 1, .., N, (25)
and [START_REF] Ramakrishna | Controllability of molecular systems[END_REF] follows upon setting

α i = θ i1 , β j = θ 1j -θ 11 .
The uniqueness of the Hamiltonian up to a diagonal unitary transformation has thus been proved. At most requiring the application of this unitary transformation to [START_REF] Kullback | Information Theory and Statistics[END_REF], we may assume H 1 = H 2 . We now prove µ 1 = µ 2 , which will establish [START_REF] Bris | Towards efficient numerical approaches for quantum control[END_REF] in general.

Step 3: We denote by σ 1 i (γ), σ 2 i (γ) respectively the i-th eigenvalue of H + γµ 1 and H + γµ 2 , with γ being a real constant, and furthermore

ω 1 ij (γ) = σ 1 j (γ) -σ 1 i (γ) and ω 2 ij (γ) = σ 2 j (γ) -σ 2 i (γ).
Using (A3) (for γ = 0) and the continuity of the spectrum with respect to γ around γ = 0, for all positive sufficiently small γ, we know that

ω 1 ij (γ) = ω 2 lk (γ) =⇒ (i, j) = (l, k).
We now choose a control field to be of the form γ + ǫ(t). In this way the Hamiltonian H for the two systems changes to H + γµ 1 and H + γµ 2 , while the dipole moments remain unchanged. Next, the same operations may be performed as in Step 1. Assumption (A5) is still valid for H + γµ 1 when γ is sufficiently small. The only modification in the argument of Step 1, lies in the fact that the transitions ω 1 ij (γ) and ω 2 ij (γ) are not necessarily the same. It suffices then to remark that in this new situation, each term of both sides of ( 19) may be again identified. This is due to the fact that, even if ω 1 ij = ω 2 ij , there does not exist any other transition in resonance with ω 1 ij and ω 2 ij . We thus obtain, instead of ( 19),

|c 1 1,i (γ)| |c 1 1,j (γ)| = 0 = |c 2 1,i (γ)| |c 2 1,j (γ)| ∀i = j.
Likewise, equation ( 24) is changed to

c 1 1,i (γ)c 1 * k,i (γ)c 1 * 1,j (γ)c 1 k,j (γ) = 0 = c 2 1,i (γ)c 2 * k,i (γ)c 2 * 1,j (γ)c 2 k,j (γ) ∀i, j = 1, .., N.
We are able to conclude that, for γ sufficiently small, there exists {α n (γ)}, {β m (γ)} :

(φ 1 m (γ)) n = e iαn(γ) e iβm(γ) (φ 2 m (γ)) n
where φ 1 m (γ) and φ 2 m (γ) are respectively the n-th eigenstate of H + γµ 1 and H + γµ 2 . Next we use non-degenerate perturbation theory (see e.g. [8, Section XI.B]) to write

(φ 1 m (γ)) n = (φ m ) n + k =m γ (µ 1 ) k,m λ m -λ k (φ k ) n + O(γ 2 ) (φ 2 m (γ)) n = (φ m ) n + k =m γ (µ 2 ) k,m λ m -λ k (φ k ) n + O(γ 2 )
In view of the order zero term, we have e iαn(γ) e iβm(γ) → 1 as γ → 0.

Using the terms up to the first order yields

φ m + k =m γ (µ 1 ) k,m λ m -λ k φ k = e iβm(γ) U (γ)φ m + k =m γ (µ 2 ) k,m λ m -λ k φ k ( 26
)
where U (γ) is the diagonal unitary transformation U (γ) = diag(e iα1(γ) , e iα2(γ) , ..., e iαN (γ) ) written in the basis of {e i } N i=1 . Taking the Hermitian product of ( 26) with φ m produces

e iβm(γ) U (γ)φ m , φ m = 1.
Since U (γ) is a unitary transformation, implying that U (γ)φ m is of norm one, it follows that U (γ)φ m = e iα(γ) φ m . Finally, taking the Hermitian product of (26) with φ k with k = m we obtain

(µ 1 ) k,n λ n -λ k = (µ 2 ) k,n λ n -λ k ∀k, n = 1, .., N
which shows µ 1 = µ 2 and concludes the proof.

Remark 10. The phases {β j } which multiply the columns of the matrix c 1 do not affect the matrix H 1 (as they are the phases which multiply the different eigenvectors of the Hamiltonian). The phases {α i } correspond to the only non-identifiable part of the system, such that (H 2 ) i,j = e i(αi-αj ) (H 1 ) i,j .

Assumption (A5) forbids the situation where a subsystem of dimension two is decoupled from the rest of the system under the free evolution (ǫ(t) ≡ 0). Theorem 3 relaxes this assumption and gives a more general result under the only additional condition that the system is of dimension greater than two.

Theorem 3. The situation is the same as that of Theorem 2, also assuming (A1), (A2), (A3), (A4). But, instead of Assumption (A5), we suppose that the system dimension N is larger than or equal to 3. Then, the same conclusion holds.

Remark 11. It is necessary to consider the dimension of the system to be greater than two, as a system of dimension two admits some additional symmetry. One can easily check that the two systems

i d dt ψ ψ 2 = α z z * β + ǫ(t) γ ξ ξ * δ ψ 1 ψ 2 ,
and i d dt

ψ1 ψ2 = (λ 1 + λ 2 ) 1 0 0 1 - α z * z β -ǫ(t) γ ξ * ξ δ ψ1 ψ2 ,
where λ 1 and λ 2 are the eigenvalues of the matrix

H 1 = α z z * β
give rise to the same populations for all times t > 0. Indeed, if (ψ 1 , ψ 2 ) is a solution of the first system, ( ψ1 , ψ2 ) = (e -i(λ1+λ2)t ψ * 1 , e -i(λ1+λ2)t ψ * 2 ) is a solution of the second system. In addition, the matrix

H 2 = (λ 1 + λ 2 ) 1 0 0 1 - α z * z β
has λ 1 and λ 2 as eigenvalues like H 1 .

Proof of Theorem 3

The only situations which need to be treated are:

(1) there is a subspace of dimension 2 spanned by two vectors in {e i } N i=1 invariant under the free evolution of the two systems (for simplicity we take these two elements to be e 1 and e 2 );

(2) there is an element of the canonical basis {e i } N i=1 (which we take for simplicity to be e 1 ) coinciding with an eigenstate of the free Hamiltonian H 1 .

In the first case the matrices H 1 and H 2 can be written

H 1 = (M 1 ) 2×2 0 0 (N 1 ) (N -2)×(N -2) , H 2 = (M 2 ) 2×2 0 0 (N 2 ) (N -2)×(N -2)
The assumption of controllability of the first system implies that the dipole moment µ 1 must couple the subsystem of dimension two spanned by e 1 and e 2 and the rest of the system. We can thus assume that, for example, the coefficient (µ 1 ) 1k is non-zero (for some 2 < k ≤ N ). Next, as in the proof of Theorem 2, we consider the matrices H 1 + γµ 1 and H 2 + γµ 2 for small positive parameters γ. The matrix H 1 + γµ 1 satisfies the assumption (A5) for a sufficiently small parameter γ = 0. Using non-degenerate perturbation theory yields

(φ 1 m (γ)) n = (φ 1 m ) n + k =m γ (µ 1 ) k,m λ m -λ k (φ 1 k ) n + O(γ 2 ) ( 27
) (φ 2 m (γ)) n = (φ 2 m ) n + k =m γ (µ 2 ) k,m λ m -λ k (φ 2 k ) n + O(γ 2 ), (28)
and, arguing as above, ∀γ > 0, there exists {α n (γ)} and {β m (γ)} :

(φ 1 m (γ)) n = e iαn(γ) e iβm(γ) (φ 2 m (γ)) n
The phase parameters α n (γ) and β m (γ) take their value in the compact interval [0, 2π]. Letting γ go to zero permits extracting a subsequence (γ i ) ∞ i=1 ց 0 such that

∀n, m = 1..N, α n (γ i ) → α n and β m (γ i ) → β m as i → ∞
for some phase parameters α n and β m in [0, 2π]. Passing to the limit in equations ( 27) and ( 28) implies

(φ 1 m ) n = e iαn e iβm (φ 2 m ) n .
This concludes the proof of [START_REF] Levis | Selective bond dissociation and rearrangement with optimally tailored, strong-field laser pulses[END_REF]. Similarly, considering the terms of first order k =m γ

(µ1) k,m λm-λ k (φ 1 k ) n and k =m γ (µ2) k,m λm-λ k (φ 2 k
) n , we may show that ( 14) holds. The second case corresponds to a situation where e 1 coincides with an eigenstate φ 1 k of the Hamiltonian H 1 . Once again, consider the two systems (H 1 + γµ 1 , µ 1 ) and (H 2 + γµ 2 , µ 2 ) for a small parameter γ. The first system is still controllable. Furthermore, it either satisfies the assumption (A5) or can be included in the first situation discussed above. As the system (H 1 + γµ 1 , µ 1 ) is identifiable passing to the limit when γ tends to zero, then the initial system (H 1 , µ 1 ) is also identifiable, which concludes the proof.

Numerical applications

An optimal identification machine

This section presents a numerical approach for the identification in the same spirit as [START_REF] Geremia | Optimal hamiltonian identification: The synthesis of quantum optimal control and quantum inversion[END_REF]. We consider the inverse problem posed in Section 1 in setting (S2), with data consisting of the populations p i in terms of the vectors of the canonical basis {e i } N i=1 and identify both the Hamiltonian H and the dipole moment µ. For simplicity, this section uses the terminology identification of parameters, where the parameters are the entries of the Hamiltonian and the dipole moment. Such parameters reside in R N ×N × R N ×N .

The numerical approach consists in iteratively correcting the control field ǫ(t) in order to restrict as much as possible the set of parameters giving rise to measurements compatible with the experimental data. As indicated in [START_REF] Geremia | Optimal hamiltonian identification: The synthesis of quantum optimal control and quantum inversion[END_REF], the approach therefore consists in a two-step iteration procedure

• inversion step: produce the field, submit the data for the inversion process and obtain a family of Hamiltonians and dipole moments compatible with the measurements; • control selection step: select new trial fields and return to the first step if the results are not satisfactory.

A global exploration of the entire set of Hamiltonians and dipole moments consistent with the data is not feasible in practice. The set A of admissible parameters (i.e., admissible matrix elements) is assumed to lie in a cube centered at the true (unknown) Hamiltonian and dipole moment, in the space

R N ×N × R N ×N , A = [H -, H + ] × [µ -, µ + ], i.e. (H, µ) ∈ A if and only if H ij ∈ [(H -) ij , (H + ) ij ] and µ ij ∈ [(µ -) ij , (µ + ) ij ]
for i, j = 1, ..., N . In real situations, the parameters will be known up to some uncertainties, which is consistent with this picture.

A single control field ǫ(t) generally gives rise to a family of different admissible parameters that all reproduce the data to a predefined tolerance threshold. Hopefully, through iteration, the solution family can be reduced to a narrow distribution of Hamiltonians and dipole moments. The identification algorithm consists in minimizing, over all possible control fields ǫ(t), the domain of the admissible parameters. To compute the diameter of the domain, one approach is to enumerate, as comprehensively as possible, the full family of Hamiltonians and dipole moments consistent with the observed data [START_REF] Geremia | Optimal hamiltonian identification: The synthesis of quantum optimal control and quantum inversion[END_REF].

With a view to decreasing the computational cost, an second approach is proposed, based on an alternative definition of the diameter. We only require the characterization of two different pairs (Hamiltonian, dipole moment): one near (H -, µ -) and the other one near (H + , µ + ). The distance between those is used as an estimate of the domain diameter.

The inversion step

Identifying the Hamiltonians that reproduce the laboratory data set, p

(lab) k = {p (lab) k,1 , ..., p (lab)
k,N } for the k-th trial field, and which are as close as possible to the extremities of A requires minimizing both the difference between the laboratory data and the calculated observables, p k,i [H, µ, ǫ k (t)] and the distance of the parameters to the limits of the set A. For this purpose, we introduce the cost functional

J(H, µ; ǫ k , H, µ) = N i=1 |p k,i [H, µ, ǫ k (t)] -p (lab) k,i | 2 |p (lab) k,i | 2 + α dist(H, H) + α dist(µ, µ) (29)
where α ∈ R is a penalization factor. In the present numerical tests α is set to 0.05. The " relative distance" is defined by

dist(H, K) = N i,j=1 |H ij -K ij | 2 |K ij | 2 . ( 30
)
Note that this is not a metric (in the usual topological sense) because it is not symmetric. Minimizing the cost functionals J(., .; ǫ k , H -, µ -) and J(., .; ǫ k , H + , µ + ) on the set A, provides two different solutions to the inversion process near the two extremities of the set A

(H inf (ǫ k , H -, µ -), µ inf (ǫ k , H -, µ -)) = argmin (H,µ) J(H, µ; ǫ k , H -, µ -) (H sup (ǫ k , H + , µ + ), µ sup (ǫ k , H + , µ + )) = argmin (H,µ) J(H, µ; ǫ k , H + , µ + ). (31)
Various approaches may be used for these minimization problems. Here we use a modified Nelder-Mead simplex [START_REF] Barton | Nelder-mead simplex modifications for simulation optimization[END_REF] for both. In order to benefit from the previous iterates, and to improve the quality of the inversion, a third external loop is added to the procedure allowing for changing the limits of the set A when some "good" candidates are found: H -and H + are replaced by H inf and H sup whenever some appropriate conditions are satisfied.

Remark 12. Different functionals may be considered to improve the definition in Eqn. [START_REF] Zhu | Potential surfaces from the inversion of time dependent probability density data[END_REF] such as the Kullback-Leibler relative entropy [START_REF] Kullback | Information Theory and Statistics[END_REF][START_REF] Kullback | On information and sufficiency[END_REF] which conveniently interprets the entries |K ij | and |H ij | in a probabilistic framework. These considerations, which are expected to be especially relevant in practice, will not be developed in detail here as we only focus on the feasibility of the data inversion.

Control selection step

The control selection step searches for a new trial field ǫ k+1 that will decrease the current value of the cost-functional

J (ǫ k ; H -, H + , µ -, µ + ) = dist (H inf (ǫ k , H -, µ -), H sup (ǫ k , H + , µ + )) + dist (µ inf (ǫ k , H -, µ -), µ sup (ǫ k , H + , µ + )) , (32)
to produce more discriminatory data. This optimization also might be done in alternative ways. More details on the approach are given in the next section.

Numerical illustration

This section considers the 4-level system where

P =     0 1 -1 1 -1 0 1 1 1 -1 0 -1 -1 -1 1 0    
is an anti-Hermitian matrix. As mentioned in Section 2, spectroscopy allows us to assume the diagonal matrix D known, again up to an overall additive constant. Thus the identification objective is to find the dipole moment µ and the rotation matrix P (or its logarithm P).

We will suppose that µ and H are initially known within 10% relative error, which corresponds to a 10% relative error on µ and about a 5% relative error on P. This initial domain size is reasonable given the quality of ab initio quantum mechanical techniques to provide such estimates.

The laser fields are chosen of the form:

ǫ(t) = β σ(t) i =j sin(ω ij t + θ ij )
where β is a fixed constant, σ(t) is a Gaussian envelope and ω ij = λ j -λ i is the transition frequency between the eigenstates i and j of H. The phases θ ij are the control parameters. More general fields can be employed with amplitudes β ij also introduced as control parameters.

the diameter of the parameter domain corresponding to the final laser field being 6.258 × 10 -3 . We then have J(P, µ; ǫ c , P + , µ + ) = J(P, µ; ǫ, P -, µ -) = 0

The solutions µ and P are determined to be which corresponds to a neighborhood of the physical matrices with 2% relative error.

P =     0
Remark 14. The penalty constant α and the tolerance constants C 1 , C 2 , C 3 need to be appropriately set. The convergence properties and the results are sensitive to their values. When these parameters are set to small values the algorithm will have difficulty recovering good solutions and may be easily trapped in a local minimum.

In contrast, values that are too large slow down the convergence, or even prevent it. At this time there is no systematic means to determine these constants, but adaptations of the algorithm to render it less sensitive to the values of parameters seem related to a more convenient definition of J (see Remark 12).

We conclude by emphasizing there is certainly room for improvement in the above numerical strategy. The present strategy only aims at demonstrating the feasibility of the approach. Definite conclusions on how this strategy performs for large size systems remains for future assessment.

i

  Ψ = (H + ǫ(t) µ)Ψ

The phase of a complex number re iφ for r > 0 is defined to be φ ∈ [0,

2π).

Acknowledgments:

The authors thank Robert L. Kosut and Pierre Rouchon for interesting discussions.

We now adapt the identification algorithm presented above for this problem. In practice we modify the cost functional J(P, µ; ǫ k , P -, µ -) defined in equation [START_REF] Zhu | A rapid monotonically convergent iteration algorithm for quantum optimal control over the expectation value of a positive definite operator[END_REF], as follows:

Here C 1 is a fixed small tolerance constant (taken initially as C 1 = .05 and decreased during the iterations as described below). The cost functional J (ǫ; P -, P + , µ -, µ + ) is defined by the iterations:

(1) Use a simplex algorithm to minimize the cost functional J(P, µ; ǫ, P -, µ -) with respect to P and µ, until J(P, µ; ǫ, P -, µ -) < C 2 (typically C 2 = 5 × 10 -3 ). This yields some nearly optimal parameters P inf (ǫ; P -, µ -) and µ inf (ǫ; P -, µ -). ( 2) Use a simplex algorithm to minimize the cost functional J(P, µ; ǫ, P + , µ + ), until J(P, µ; ǫ,

This yields some P sup (ǫ; P + , µ + ) and µ sup (ǫ; P -, µ -). (3) Set the cost functional

The following algorithm is used to identify the Hamiltonian and the dipole moment:

(1) Set D 1 = dist(P + , P -) and D 2 = dist(µ + , µ -).

(2) Initialize the laser field ǫ(t) by randomly choosing the phases:

(3) Use a simplex algorithm in order to minimize the cost function J (ǫ; µ -, µ + , P -, P + ) with respect to ǫ(t) and to find a ǫ(µ -, µ + , P -, P + ) and its corresponding Hamiltonians and dipole moments: P inf (ǫ), P sup (ǫ), µ inf (ǫ) and µ sup (ǫ). At each iteration of the simplex compute dist(P sup (ǫ), P inf (ǫ)) and verify if it is less than D 1 or not. If it is the case replace P -by P inf (ǫ) and P + by P sup (ǫ), and then update D 1 . At the same time we update C 1 by dividing it by √ 2 (such an update leads to a more restrictive stopping criterion for the internal simplex). The same operations are done for the dipole moments: i.e. compute dist(µ sup (ǫ), µ inf (ǫ)) and verify if it is less than D 2 or not. If it is the case we replace µ -by µ inf (ǫ) and µ + by µ sup (ǫ), update D 2 and C 1 and continue the simplex. The simplex stops when the two following criteria are satisfied:

• J (ǫ) < C 3 • J(P inf , µ inf ; ǫ, P -, µ -) = 0 or J(P sup , µ sup ; ǫ, P + , µ + ) = 0 where C 3 is a small tolerance constant (here C 3 = .01). This algorithm has been tested on system (33).

Remark 13. The tolerance constants C 1 and C 2 are relevant for the numerical resolution of the optimization problem but can also mimick the impacts of noise. The parameter C 1 is a measure of the data error and it decreases from a large value (convenient from the numerical point of view) to a final non-null value reflecting the nature of the data errors being finite.

Note also that C 1 is relevant for ensuring that the inverted values remain inside some given range where the search takes place. For a perfect numerical optimization method no bounds would be required since the main result of this paper guarantees that the solution is unique. But, once the numerical effort is taken into account the bounds not only work to accelerate convergence, but also help to take into account the non-uniqueness induced by the presence of the noise. In practical situations, with noise included, the bounds on H 0 and µ should be set off as separate from the data error bound. That is, C 1 would be split into two constants, one to deal with data errors and the other with the parameter bounds.

The simulation time is T = 20000 which corresponds to nearly 100 periods of the transition with the smallest frequency. The constant amplitude β of the laser fields is fixed to be 0.1 and the Gaussian envelope is σ(t) = exp(0.7 × 10 -7 (t -T 2 ) 2 ). After 7 iterations of the laser field ǫ the external simplex algorithm converges,

Appendix A. Controllability

This short appendix has for subject to clarify the controllability assumption A1 needed for the uniqueness result of the Section 2. At a first step, we will give a definition for this controllability criteria and then we will provide a brief review over the existing results on this subject. We restrict ourselves to the finite dimensional configuration being the case of interest in this paper. For the infinite dimensional setting however, much less results are available up to date.

A finite dimensional quantum system (1) is said to be wavefunction controllable if every pair of initial and final state Ψ 0 and Ψ 1 in S N -1 (the unit sphere of C N ) there exists a control function ǫ and a time T > 0 such that the solution of (1) at time T , with initial state Ψ 0 , is Ψ(T ) = Ψ 1 .

In [START_REF] Ramakrishna | Controllability of molecular systems[END_REF] a very natural sufficient condition for the wavefunction controllability of a finite dimensional quantum system is given:

• If the Lie algebra spanned by the skew-Hermitian matrices (H 0 + V )/i and µ/i contains the Lie algebra of all N × N skew-Hermitian matrices, then the system (1) is wavefunction controllable. In a more recent work Turinici and Rabitz [START_REF] Turinici | Quantum wavefunction controllability[END_REF] provided an easier to check sufficient condition utilizing the connectivity graph of the dipole moment µ.

Finally Albertini and D'Alessandro [START_REF] Albertini | Notions of controllability for bilinear multilevel quantum systems[END_REF] gave a necessary and sufficient condition also based on the Lie algebra spanned by (H 0 + V )/i and µ/i characterizing the wave function controllability of the system (1).