Tom Douce

Michael Stern

Nicim Zagury

Patrice Bertet

Pérola Milman

Coupling a single nitrogen-vacancy center to a superconducting flux qubit in the far-off-resonance regime

Keywords: number(s): 03.67.Lx, 03.67.Bg, 42.50.Ex

We present a theoretical proposal to couple a single nitrogen-vacancy (NV) center to a superconducting flux qubit in the regime where both systems are off resonance. The coupling between both quantum devices is achieved through the strong driving of the flux qubit by a classical microwave field that creates dressed states with an experimentally controlled characteristic frequency. We discuss several applications such as controlling the NV center's state by manipulation of the flux qubit, performing the NV center full tomography and using the NV center as a quantum memory. The effect of decoherence and its consequences to the proposed applications are also analyzed. Our results provide a theoretical framework describing a promising hybrid system for quantum information processing, which combines the advantages of fast manipulation and long coherence times.

I. INTRODUCTION

Electronic spins in semiconductors such as nitrogenvacancy (NV) centers in diamond or impurities in silicon are characterized by their low decoherence rates s [1][2][3][4]. They would be excellent candidates for realizing an operating quantum processor if only there was a way to reliably entangle distant spins [5,6]. Recently, it has been proposed to solve this issue by using a superconducting circuit as a quantum bus [7,8]. Such hybrid systems combine the advantages of fast manipulation and long coherence times [9,[START_REF] Kurizki | Proc. Natl. Acad. Sci. USA[END_REF]. However, this approach would require reaching the strong coupling regime where the coupling strength g between the spin and the circuit is larger than their respective decoherence rates (g { s , circuit }). To achieve this goal, the coupling constant g should be increased by several orders of magnitude compared to the current state of the art where g is typically of the order of a few hertz [START_REF] Kubo | [END_REF]12].

The coupling constant g can be greatly increased by using flux qubits instead of a linear superconducting resonator. Indeed, flux qubits [13][14][15][16][17][18][19][20] are characterized by a macroscopic permanent current I P typically of the order of several hundreds of nanoampere. This current flows in the loop of the qubit and generates therefore a large magnetic dipole which allows an efficient coupling to spins. Bringing the spin at a distance of ∼20 nm from the flux qubit can be achieved by fabricating an ultranarrow superconducting wire and aligning it very precisely with the spin (see Fig. 1). This would allow increasing the coupling by several orders of magnitude and reaching a coupling constant g/2π ∼ 100 kHz, a value much larger than recently reported spin decoherence rates [1,2]. Flux qubits with decoherence rates smaller than 100 kHz have been demonstrated [21,22]. Capacitively-shunted flux qubits [23,24] with decoherence rates of order 10 kHz were also demonstrated recently [25], although their persistent current is in general lower. However, flux qubits suffer from a severe limitation. Their large magnetic dipole makes them very sensitive to flux noise which limits their coherence times. To cope with this problem, one should tune the flux threading the loop precisely at half a flux quantum ( = 0 /2), the so-called optimal point of the flux qubit [15]. At this point its energy does not depend on the flux at first order. The coupling of a spin-whose transition energy ω s is given by nature-with a flux qubit is optimal for ω s = . Unfortunately, reaching this target is a challenging task. The value of the so-called flux-qubit gap is an exponential function [17] of the parameters of the junctions which form the qubit and it is therefore poorly controlled. It is possible to design a tunable flux qubit by replacing one of the junctions by a superconducting quantum interference device (SQUID) [26]. However, coherence times required for single spin coupling have not been reported so far in this geometry. Applying an external magnetic field to tune the spin close to frequency also has many limitations. The critical field of aluminium is ∼100 G and therefore limits the tunability of the spin frequency to 100-200 MHz. As a consequence the two resonance frequencies ω s and may be highly detuned. Typically, one can expect that (ω s )/2π ∼ 300-500 MHz, prohibiting the use of the direct resonant interaction to entangle the two systems. In this work, we propose a theoretical strategy to cope with this issue and couple efficiently both systems even with such a detuning.

II. THE MODEL

A. Effective Hamiltonian

We now present our model, starting with NV centers, the flux qubit (FQ), and the coupling mechanism. NV centers are diamond color centers consisting of a nitrogen substitutional impurity next to a vacancy [27,28]. In their negatively charged state, they have an electronic ground state with a spin S = 1 which is described by the following Hamiltonian (neglecting the effect of strain [29]):

Ĥ s = DS 2 z + γ e S • B ext , ( 1
)
where S are the dimensionless spin-1 operators, D/(2π ) ≈ 2.88 GHz is the zero-field splitting ground state, and the rightmost term is the Zeeman interaction of the electronic spin of gyromagnetic ratio |γ e |/2π = 28 GHz/T. The quantization z axis is set in the nitrogen vacancy axis. B ext is an external magnetic field applied to lift the degeneracy of the {m S = ±1} levels. For an antiparallel field, the NV center is well described by a two-level system Hamiltonian in the {m S = 0,m S = +1} basis:

Ĥ s = ω s σ s z /2 , ( 2
)
where ω s = D + γ e B ext , σ s z = |1 1| -|0 0| with |0 and |1 corresponding, respectively, to m S = 0 (lower state) and m S = 1.

The flux qubit Hamiltonian can also be described by a two-level system [17]:

Ĥ f q = σ f q x /2 + σ f q z /2 , (3)
where = 2I p ( -0 /2). In order to reach the longest coherence times the flux qubit must be manipulated at the so-called optimal point with = 0 /2. Consequently we have = 0 and the flux qubit transition frequency is given by . The NV is in a position such that the field from the flux qubit is normal to the quantization axis of NV spin S. Therefore the Hamiltonian for the two interacting systems can be written as

Ĥ = ω s σ s z /2 + σ f q z /2 + gσ s x σ f q x , (4)
where the fourth term is the Zeeman interaction of the electronic spin due to the magnetic field generated by the flux qubit.

The coupling constant g is given by g γ e μ 0 I p /( √ 2 • 2πr) where μ 0 is the vacuum permeability and r the distance between the spin and the constriction shown in Fig. 1. By setting r = 15 nm and I p = 500 nA, we get g/2π 100 kHz.

The flux qubit is subjected to an additional microwave magnetic field of frequency ω in a direction normal to the plane of the persistent currents. The NV spin state has a much smaller magnetic dipole and therefore we will assume that it is not driven directly by the microwave field. Thus the Hamiltonian reads

Ĥ = ω s σ s z /2 + σ f q z /2 + g σ s x σ f q x + σ f q x cos(ωt), (5)
where = I P AB 0 is the Rabi frequency of the flux qubit of area A and persistent current I P driven by the classical microwave field B(t) = B 0 cos(ωt).

From (5) we can move to a frame rotating with the classical field's frequency ω, which is set to be the same as the flux qubit transition frequency-namely ω = . This is accomplished through the unitary operator transformation U s (t) = e -i t( σ s z +σ f q z )/2 U r (t), connecting the time evolution operators in the Schrödinger and rotating frame pictures. We define δ = ω sω = ω sand the raising and lowering operators in the σ z basis, σ± = ( σx ± i σy )/2. By setting |δ|, we identify two terms in the resulting Hamiltonian in the rotating frame:

Ĥr = Ĥ0 + Ĥint , ( 6
) with Ĥ0 = δ 2 σ s z + 2 σ f q x , Ĥint = g( σ s + e i t + σ s -e -i t )( σ f q + e i t + σ f q -e -i t
). ( 7) We can thus consider that the classical drive "dresses" the flux qubit [30,31]. From now on this representation will be referred to as the laboratory frame. In order to understand the physical meaning of this time-dependent interaction Hamiltonian, we plotted in Fig. 2 7). We can see the slow dynamics induced by the magnetic coupling while the inset displays the fast Rabi oscillations occurring with a much shorter time scale. The set of parameters is ω s = 28 800g, = 25 800g, and = -δ.

flux qubit and the spin. For the numerical simulation we have set ω s = 28 800g and = 25 800g so that δ = 3000g, with = -δ.

We now move to the interaction picture with respect to Ĥ0 in order to get an approximate expression for the interaction Hamiltonian. We have

ĤI = g 2 σ s + e i(2 +δ)t σ f q x -e i t σ f q z -i σ f q y 2 + e -i t σ f q z + i σ f q y 2 + σ s + e iδt σ f q x -e -i t σ f q z + i σ f q y 2 + e i t σ f q z -i σ f q y 2 + H.c . ( 8
) Let us define σ+,x = |+ -| = ( σz -i σy )/2 and σ-,x = |-+| = ( σz + i σy )/2, where |± = (|1 ± |0 )/ √ 2.
We recall that we focus on the case with |δ| in Eq. ( 8). Assuming also the conditions ( , ,δ g) we can perform the so-called rotating wave approximation and are led to the effective Hamiltonian:

Ĥeff (t) = ± g 2 σ s + σ f q ±,x e i( ±δ)t + σ s - σ f q ∓,x e -i( ±δ)t . ( 9
)
The effective Hamiltonian ( 9) describes a closed two-qubit system that can be diagonalized and studied analytically. It can lead to a number of applications that will be further detailed in the next sections. We stress that in the case of exact matching condition = ±δ the effective Hamiltonian becomes time independent. In the following we set = -δ so that Ĥeff = g 2 ( σ s

+ σ f q +,x + σ s - σ f q -,x
). Adapting the results to = δ is straightforward. Supposing that the spin-flux qubit bipartite system is initially in a separable pure state, it can be written as |ψ(0) = (cos θ |0 + e iϕ sin θ |1 )(cos α|+ + e iφ sin α|-). [START_REF] Kurizki | Proc. Natl. Acad. Sci. USA[END_REF] In [START_REF] Kurizki | Proc. Natl. Acad. Sci. USA[END_REF], we dropped subscripts and used a convention that will be kept the same throughout this manuscript: The spin's state is expressed in the basis of σ s z eigenstates |0 ,|1 , and the flux qubit state is expressed in the basis of σ f q x eigenstates |+ ,|-. This basis choice is well adapted to the coupling induced by (9). Moreover, it provides a clear distinction between each party's states without using auxiliary labels. The parameters θ and ϕ characterizing the spin's state are unknown and cannot be controlled or manipulated without the coupling (9). However, the parameters α and φ can be experimentally controlled and the flux qubit can be prepared in an arbitrary initial state by combining classical pulses with different phases, intensity, and duration [18]. The time evolution of the initial state [START_REF] Kurizki | Proc. Natl. Acad. Sci. USA[END_REF] subjected to Hamiltonian (9) is |ψ(t) = cos θ cos α|0,+ + e i(φ+ϕ) sin θ sin α|1,-

+ e iϕ sin θ cos α -i sin gt 2 |0,-+ cos gt 2 |1,+ + e iφ cos θ sin α cos gt 2 |0,--i sin gt 2 |1,+ . (11)
In order to test the validity of our model, we compare the time evolution induced by the effective Hamiltonian (9) to the evolution under the action of the exact Hamiltonian in the interaction picture (8). This can be done by comparing the numerical integration of the Schrödinger equation using (8) to the analytical diagonalization of ( 9). We can study as an example the case where θ = π/2 and α = 0. In this case, the initial state is |1 |+ using the basis convention. This state is an eigenstate of the free Hamiltonian Ĥ0 , so the interaction Hamiltonian is solely responsible for the dynamics. The time evolution of the population in |1 |+ under either Ĥeff or Ĥint is shown in Fig. 3, with = -δ (left) and = -δ + g (right). We can see that indeed, both curves display a very good overlap. Thus the effective Hamiltonian is a good approximation to the coupling between both systems, even if the resonance is not perfect, i.e., = ±δ, in this case leading to a reduction of the fringes visibility.

The tuning of the resonance can be done by continuously changing the intensity and phase of the classical dressing field until coupling between the NV center and the flux qubit is optimal, as in a "quantum radio." This point will be further discussed in Sec. III A.

We should also mention that similar ideas have been proposed to couple superconducting qubits with widely different frequencies [32][33][34].

B. Decoherence

The dynamics of the coupled spin-flux qubit system will be influenced by the presence of noise, that creates population losses and dephasing. In order to describe the proper dynamical equation for the density matrix of the system, we will start from a microscopic description of the flux qubit-environment interaction and then derive the corresponding Lindblad equation. Given the typical decay rates of a flux qubit compared to those of an NV center, we will neglect contributions to decoherence induced by the dephasing of the latter. Our microscopic model is thus made of two interactions: a dissipative process conveyed by the qubit σx operator and a dephasing one due to σz . More precisely we have

Ĥdiss = k γ k σ f q x ( bk + b † k ), (12)
and [35] Ĥdeph

= k λ k σ f q z ( ĉk + ĉ † k ), (13)
where the b( †) k 's and ĉ( †) k 's are bosonic annihilation (creation) operators of the modes of the environment. In our study we will assume the environment is in thermal equilibrium at zero temperature.

The flux qubit undergoes two independent nondissipative evolutions: The first is due to the microwave drive, at frequency , the second to the coupling with the NV center at frequency g/2. The latter is the one we are interested in while the former is included in the free Hamiltonian Ĥ0 . We derive the master equation based on these ingredients in the so-called Born-Markov approximation [START_REF] Carmichael | Statistical Methods in Quantum Optics[END_REF]. The Lindblad equation should model the effect of damping and dephasing in the basis of the free Hamiltonian eigenstates. Since the free Hamiltonian of the flux qubit is proportional to σx , the Lindblad operators will be composed of raising and lowering operators in the Pauli-σx basis. Specifically, we get the following dynamical equation for the system density matrix in the laboratory frame:

d ρ dt = -i[ Ĥ0 + Ĥeff , ρ] + x L σ f q x ρ + -L σ f q -,x ρ + + L σ f q +,x ρ, ( 14
)
where L is the Lindblad superoperator

L[ â] ρ = â ρ â † - 1 2 ( â † â ρ + ρ â † â).
Notice that in contrast to the standard approach the microwave drive rotates the decoherence basis:

The ground state is a thermal distribution in the {|-,|+ } basis. Even though the environment was assumed to be at zero temperature, the driving field leads to an effective heating of the bath. Moreover, the i 's decay rates are related to environmental spectral properties and thus to experimentally accessible quantities [START_REF] Ithier | [END_REF]. Let us define T 1 as the energy relaxation time and T ν the Rabi oscillations decoherence time which can be measured through independent experiments. Then we have

x = 1 4 T -1 1 , -= 1 4 T -1 1 , + = 1 4 T -1 1 + T -1 ν . ( 15
) One should notice the asymmetry between + and -that will have peculiar consequences later on. With the dynamical equation ( 14) we are able to numerically evaluate the time evolution of the system, including in our simulations realistic values for T 1 and T ν , experimentally determined, e.g., in [21]. There we have T 1 = 10-20 μs and T ν = 10-15 μs. In Fig. 4 we show the effect of decoherence on the reduced flux qubit state obtained after tracing over the spin state. We plot the probability of finding the flux qubit in the excited state, P 1 (t) = Tr spin ( 1 f q | ρ(t)|1 f q ), where Tr spin denotes the partial trace over the spin's degree of freedom.

We will now show how this formalism can be used in fundamental quantum information protocols.

III. APPLICATIONS IN QUANTUM INFORMATION PROTOCOLS

We will see in the following how a proper choice of the parameters α and φ can lead to a number of applications relying on the coupling between the two-level systems, such as manipulation and measurement of the spin's state.

A. Detecting spin-qubit coupling

The first step in view of exploring the flux qubit-single spin coupling is to detect the coupling itself. Here we study a protocol enabling this in the presence of decoherence, which introduces noise and loss of information about the flux qubit's state.

We first provide the exact calculation neglecting dissipation and restricting to zero detuning, i.e., = -δ, focusing on the interaction picture. We define P i 1 (t) = Tr spin ( 1 f q | ρi (t)|1 f q ), where ρ i (t) is the density matrix of the whole system computed in the frame rotating with Ĥ0 . Based on Eq. ( 11) we have

P i 1 (t) = 1 2 1 + sin 2α cos φ cos gt 2 . ( 16
)
We recall that the main goal of the present protocol is to detect the existence of the coupling. In this regard, the key point of the dynamics of the population P i 1 (t) is that it does not depend on the NV center state, which is a priori unknown, but rather on the flux qubit state which is controllable. This feature will not be preserved when including decoherence to the dynamics. Hence we can set α = π/4 (i.e., starting with the flux qubit in the excited state |1 ) to get the highest visibility. In that case P i 1 (t) reads

P i 1 (t) = 1 2 1 + cos gt 2 . ( 17
)
This equation is valid provided the flux qubit and the spin are in perfect resonance = -δ. Furthermore, the resonance condition ensures that P i 1 (t) corresponds to the envelope of the fast oscillations when moving back to the laboratory frame (see Fig 5,top left). Equation (17) corresponds to an oscillation with unit visibility. The fast oscillations at frequency come from the classical driving field which induces the free Hamiltonian of Eq. ( 7). They are therefore of no interest regarding the coupling with the spin. In practice, exact resonance must be determined experimentally. In this respect, driving the flux qubit off-resonance reduces the visibility of the coupling oscillations, as can be seen in Fig. 5. So one practical way to find out the resonance is by modifying and searching for high visibility oscillations. It is possible to gradually vary until the slowly oscillating terms reach a population inversion dynamicspractically cross the P 1 = 1/2 horizontal line, which uniquely characterizes the resonance condition. Figure 5 shows the time evolution of P 1 (t), defined in Sec. II B, for different values of the detuning. It is clearly possible to identify when the resonance condition is matched. Moreover, the amplitude of the oscillations characterizing the coupling are very sensitive to the detuning-typically for values as low as g/2-which guarantees the precision of this protocol.

However, decoherence may significantly damp the oscillations and destroy the expected signal related to the population measurement of state |1 . Furthermore, as will be shown later, the dynamics gets to depend on the spin's initial state. In order to characterize the effects of decoherence we studied the behavior of P 1 (t) in the presence of decoherence for different values of the system coherence times and detuning (see Fig. 6). We set the initial state of the spin to the ground state |0 and based the numerical analysis on the Lindblad equation derived in Eq. ( 14). It is clear from Fig. 6 that even when including decoherence, the coupling of the flux qubit with the spin can be identified through the measurement of oscillations in the excitation probability of the flux qubit. The sensitivity to the detuning is sustained, which enables searching for the resonance condition by gradually varying .

One may ask how low the coherence times can be in order to see the slow oscillations due to the coupling with the spin. The coupling should be strong enough with respect to decoherence to allow for the experimental observation of at least one oscillation. In Fig. 7, we show the oscillations at resonance condition for decreasing values of T 1 , while T ν is set to 15 μs. The initial state is |0 s |1 f q . Surprisingly enough, the figures reveal that the coupling can be identified for rather low values of T 1 such as 2 μs, that is, gT 1 0.2.

An interesting feature of the coupling also appears when dissipation is considered: Depending on the spin initial state, it can either protect the flux qubit from the effects of decoherence or enhance them. This effect is illustrated in Fig. 8 where two different initial states of the spin are considered, leading to different damping behaviors. This is a consequence of the specific form of the master equation Eq. ( 14), which exhibits an asymmetry between the dissipation terms corresponding to σ +,x and σ -,x . We could indeed verify that the dependence on the spin's initial state vanishes when artificially imposing + = -.

B. Spin state initialization

A first application of the coupling described in (9), that is also the first of the DiVincenzo [38] criteria to define a qubit, is the spin state initialization.

By setting α = 0 in [START_REF] Kurizki | Proc. Natl. Acad. Sci. USA[END_REF], we have that the initial state of the hybrid system is given by

|ψ(0) = (cos θ |0 + e iϕ sin θ |1 )|+ . ( 18
)
From the time evolution induced by the coupling between the spin and the flux qubit, after t = π/g, this state will be set to state |0 . From this point on, the NV center can either be manipulated and set to an arbitrary state, or can be used as a quantum memory, as developed in the next subsection.

We now analyze the impact of decoherence on the initialization protocol. The results are shown in the left column of Table I. Even though we consider the NV center itself to be perfectly isolated from the environment, the purity is reduced by the interaction with the flux qubit and leads to the preparation of the nonpure state instead of the target one |0 . In Table I we call fidelity Tr flux ( 0 s | ρ(π/g)|0 s ), where Tr flux denotes the partial trace over the flux qubit's degree of freedom. In the presence of high decoherence rates, one way to improve the state initialization is by repeating the protocol. After a spin-flux qubit interaction time of t = π/g, one measures the flux qubit state disregarding the results of the measurement. Mathematically, this measurement projects the flux qubit and destroys entanglement with the spin, which corresponds to TABLE I. Average fidelities for the initialization protocol. Details on how they were computed can be found in the main text.

Typical times

Average fidelity Average fidelity (T 1 ,T ν ) (in μs) of the protocol after five iterations (10,10) 0.92 0.95 (20,15) 0.96 0.97 tracing over the measurement eigenbasis. Then one may again prepare the flux qubit in |+ . We then let the spin and the flux qubit interact for the same duration. The repetition of the protocol improves the state initialization procedure, until it reaches an asymptotic value independent of the initial spin state-right column of Table I.

The values of the fidelity between the final reduced density matrix of the spin and the projector on its ground state-see Table I-were numerically calculated based on the dynamical evolution Eq. (14). In order to compute those values, we repeated the protocol described in the former paragraph and then averaged the fidelity we obtained over 500 Haar-random initial spin states.

C. Quantum memory and spin state manipulation

We have shown in Sec. III B how to initially prepare the spin in state |0 . From this initialized state, we show now how the spin can be used as a quantum memory, encoding in a long lifetime quantum system the state of the flux qubit. We will also see how to adapt this strategy to manipulate the spin's state, realizing arbitrary single-qubit rotations.

We start with the quantum memory protocol. If the initial spin state is |0 , we have that θ = 0 in [START_REF] Kurizki | Proc. Natl. Acad. Sci. USA[END_REF] and the flux qubit is prepared in an arbitrary state that should be perfectly transferred to the spin.

|ψ(0) = |0 (cos α|+ + e iφ sin α|-). (19)
Using [START_REF] Kubo | [END_REF], we see that after a time t = π/g with the initial condition (19), the SC quantum state will be completely transferred to the spin state. Thus, the latter can play the role of a quantum memory, since it has a longer coherence time than the flux qubit. Figure 9 shows the fidelity of the final spin state after this protocol with respect to the initial flux qubit state when decoherence is also considered. The computed fidelity does not depend on the phase φ, because the decoherence process itself is φ independent. The maximum fidelity is obtained for α = 0, that is, for an initial state |0 |+ which is an eigenstate of the coupling Hamiltonian of eigenvalue 0. The minimum is for α = π/2 which implies a complete excitation transfer and is therefore the most likely to be impacted by decoherence.

The strategy above enables the spin state manipulation as well. Suppose that one wants to apply a given rotation to an arbitrary spin state, that can be expressed as in (18). Any rotation can be associated with the following transformation of the basis states of the NV center:

|0 = | 0 → cos β|0 + e iχ sin β|1 , ( 20
) |1 = | 1 → cos β|1 -e -iχ sin β|0 .
This transformation can be achieved by the following protocol: In a first step, the NV center state is reinitialized, as described in Sec. III B. Thus, the state of the final coupled system composed by the spin and flux qubit is given by

|ψ i = |0 (cos θ |+ + e iϕ sin θ |-), (21)
i.e., the state of the spin is completely transferred to the flux qubit. After state transfer, the coupling between the spin and the flux qubit can be stopped by turning off the intense dressing classical pulse. The flux qubit can then be manipulated by Finally, the dressing microwave field can be turned on, coupling again the spin to the flux qubit. After t = π/g, we obtain the final state:

|ψ r = cos θ | 0 + e iϕ sin θ | 1 , ( 23
)
which corresponds exactly to the realization of an arbitrary rotation to the initial state (18).

D. NV center state tomography

The previously discussed strategies can also be used for realizing the full spin state tomography. Since we have shown that it is possible to transfer the-unknown-state of the NV center to the flux qubit, one can, after this operation, simply realize the full tomography of the flux qubit.

Suppose now the initial state is such that α = 0 and φ = 0, so the flux qubit is in |+ . Then for a pulse duration t = π/g, the state becomes, according to Eq. [START_REF] Kubo | [END_REF], |ψ(t) = |0 (cos θ |+ + ie iϕ sin θ |-).

(

) 24
The unknown state of the NV center has been transferred to the accessible flux qubit. Full tomography of the latter yields perfect knowledge about the initial state of the NV-center.

We have also studied the fidelity of this protocol in the presence of decoherence, finding that it displays the same behavior as the one of the initialization protocol.

IV. CONCLUSIONS

We have shown that a hybrid system, composed of a directly inaccessible spin and a superconducting flux qubit, can be effectively coupled by driving the flux qubit with an intense classical microwave field even in the limit where both two-level systems have far-off-resonant characteristic frequencies. The coupling is created by the dressing of the flux qubit by the classical field. Such dressing leads to the possibility of tuning the dressed eigenvalues to the frequency difference between the two quantum devices, a process that can be described by an effective Hamiltonian. The possibility of coupling such devices, that present complementary advantages with respect to quantum information processing, leads to a number of applications, discussed in the present paper. We have developed protocols to manipulate the spin state, use it as a quantum memory, and realize its full tomography. In all protocols, a detailed study of the effects of decoherence in the dressed system was included, establishing limits on the expected fidelities according to decoherence rates compatible to the state of the art. Our results serve as a road map to promising experiments using hybrid quantum devices. An interesting perspective is studying how the flux qubit can serve as a data bus, intermediating the coupling between the spin and the quantum field of a resonator, both in the strong and the ultrastrong limits [39,40].

FIG. 1 .

 1 FIG. 1. (Color online) Schematic view of a flux qubit coupled to a single spin. The flux qubit consists of a superconducting loop with a constriction and intersected by four Josephson junctions (inred). Three of the junctions are identical while the fourth junction is smaller than others. When the flux threading the loop is close to half a flux quantum, the system behaves as a two-level system The nitrogen vacancy center represented as a golden arrow is situated in the close vicinity (r = 15 nm) of the constriction to maximize the coupling between the two systems. The application of an external magnetic field B ext in the direction of the NV axis lifts the degeneracy between the m S = ±1 states.

1 FIG. 2 .

 12 FIG. 2. (Color online)Probability of finding the composite system in state |1 |1 as a function of time (in units of 2π/g). We numerically solve the Schrödinger equation with Hamiltonians Ĥ0 and Ĥint from Eq.(7). We can see the slow dynamics induced by the magnetic coupling while the inset displays the fast Rabi oscillations occurring with a much shorter time scale. The set of parameters is ω s = 28 800g, = 25 800g, and = -δ.

  FIG.3. (Color online) Probability of finding the composite system in state |1 |+ as a function of time (in units of 2π/g). In dashed blue we show the exact time evolution of population while in continuous red the one obtained using the effective Hamiltonian (9). (Left) Case of perfect resonance, = -δ. (Right) Nonresonant coupling δ + = g. We can see that in both cases the effective evolution matches the exact one with very good accuracy.

1 FIG. 4 .

 14 FIG. 4. (Color online) Probability of finding the flux qubit in the excited state as a function of time (in units of 2π/g), for a spin initial state |1 . Time evolution is computed in the laboratory frame. The plot is at resonance condition; decoherence is T 1 = 20 μs, T ν = 15 μs. (Inset) Zoom on the fast Rabi oscillations for very short times.

1 FIG. 5 .

 15 FIG. 5. (Color online) Probability of finding the reduced flux qubit density matrix in the excited state |1 as a function of time (in units of 2π/g). Decoherence here is not taken into account. From top left and increasing clockwise, values of the detuning are + δ = 0,g/2,g, and 2g.

1 FIG. 6 .

 16 FIG. 6. (Color online) Probability of finding the reduced flux qubit density matrix in the excited state |1 as a function of time (in units of 2π/g). Decoherence is modeled through T 1 = 20 μs, T ν = 15 μs. Black crosses show the envelope of the fast oscillations in the absence of coupling g = 0. From top left, and increasing clockwise, values of the detuning are + δ = 0,g/2,g, and 2g.

1 FIG. 7 .

 17 FIG. 7. (Color online) Probability of finding the reduced flux qubit density matrix in the excited state |1 as a function of time (in units of 2π/g). T ν = 15 μs while T 1 is 20, 10, 5, and 2 μs (from top left, decreasing clockwise). Black crosses show the envelope of the fast oscillations in the absence of coupling g = 0.

1 FIG. 8 .

 18 FIG. 8. (Color online) Probability of finding the reduced flux qubit density matrix in the excited state |1 as a function of time (in units of 2π/g). Decoherence is modeled through T 1 = 20 μs, T ν = 15 μs. Black crosses show the envelope of the fast oscillations in the absence of coupling g = 0. (Left) Spin initial state is |0 . (Right) Spin initial state is |1 .

FidelityFIG. 9 .

 9 FIG. 9. (Color online) Fidelity of the NV center reduced density matrix after the protocol with the initial flux qubit state (a) as a function of α and φ and (b) as a function of α only. (a) Shows that the fidelity of the protocol is actually φ independent. Values of decoherence for the two plots in (b) are the same as in TableI.

ACKNOWLEDGMENTS P.B. acknowledges support from the European Community's Seventh Framework Programme (FP7/2007-2013) through European Research Council Grant No. 615767 (CIRQUSS). N.Z. acknowledges support from the Brazilian agency CNPq and the National Institute of Science and Technology for Quantum Information.


Armelle Bauer

Ali-Hamadi Dicko

Franc ¸ois Faure

Olivier Palombi

Jocelyne Troccaz

Ladaf

Ljk

Anatomical Mirroring: Real-time User-specific Anatomy in Motion Using a Commodity Depth Camera

Keywords: User-specific anatomy, Augmented Human, Real-Time, Motion Capture, Augmented Reality, Markerless Device. Concepts:, Computing methodologies → Motion capture; Mixed / augmented reality;

Figure 1: Using a single Kinect sensor output (color map (a), motion capture skeleton (b), point cloud (c)) and a generic 3D anatomical model (d), a user-specific anatomy is generated and animated in real-time (e).

Introduction

The emergence of commodity depth cameras such as Kinect sensors motivates new educational, medical and healthcare applica-tions. However, previous studies show that raw Kinect data cannot be easily employed in human motion tracking [START_REF] Pfister | Comparative abilities of microsoft kinect and vicon 3d motion capture for gait analysis[END_REF][START_REF] Malinowski | On using the microsoft kinect sensors to determine the lengths of the arm and leg bones of a human subject in motion[END_REF]. In this paper, a new calibration and motion capture sufficiently accurate for AR applications are introduced and demonstrated by superimposing internal anatomy on the user's color map in real-time. At calibration time, the length and width of body segments are estimated based on specific body poses and silhouettes. A novel anatomically sound deformer is applied to fit a high-quality generic 3D biomechanical model in order to generate a user-specific anatomical model. At run-time, our model tracks bone motions based on the Kinect body joints output, while enforcing anatomical plausibility rules such as constant lengths and joint limits. Our user study preliminary shows that the precision is sufficient to superimpose the user-specific anatomical model onto the color image, using linear blend skinning.

The paper is organized as follows: Section 1 briefly survey related work. Section 2 introduces a body size measurement procedure based on multiple poses and silhouette points and an anatomically sound deformer well adapted to Kinect outputs. Section 3 describes how it is animated based on robust motion capture using anatomical constraints. Section 4 goes through results and the validation process. Section 5 finally concludes by presenting possible applications of this work and under development features.

Related Work

Nowadays, human body modeling and tracking are widely studied for a variety of applications such as motion capture or morphometric studies.

Skin Registration is the most accessible approach to generate a wide range of human bodies. Most studies are based on skin statistical models generated by a shape and pose database [START_REF] Helten | Personalization and evaluation of a real-time depth-based full body tracker[END_REF]]. [START_REF] Gilles | Frame-based elastic models[END_REF]] use frame-based skinning methods to deform a generic skin to fit at best the user data. Other approaches using point cloud [START_REF] Li | 3d self-portraits[END_REF] or multiposition silhouettes [START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF]] may also be used to reconstruct the body skin. Most often, raw data come from acquisition of people wearing clothes and this may lead to non-realistic bodies (part proportions, etc). [Bȃlan and Black 2008] and [START_REF] Zeng | Estimation of human body shape and cloth field in front of a kinect[END_REF] intend to find ways to pass through these limitations. Since they rely only on skin models and do not include internal anatomy, those methods may result in unrealistic skin twisting.

Anatomy Registration

The most accurate subject-specific anatomy registration methods come from the medical imaging fields [START_REF] Sotiras | Deformable Medical Image Registration: A Survey[END_REF]]. However, 3D medical images are not easilly used in a non medical context and are not adapted to realtime capture. Several other methods have been proposed. [START_REF] Quah | 3d modeling of humans with skeletons from uncalibrated wide baseline views[END_REF]] present a pose-dependent method to register a 3D anatomical model onto 2D images. Based on key points, they register skin and skeleton (no soft tissues). However this method gives static results. Using Kinect point cloud, [START_REF] Zhu | Adaptable anatomical models for realistic bone motion reconstruction[END_REF]] register user-specific skin and skeleton during motion. [START_REF] Dicko | Anatomy transfer[END_REF] as [START_REF] Gilles | Creating and animating subject-specific anatomical models[END_REF], present a pose-dependent method to transfer 3D anatomical model to a target skin. This method is static and time consuming. The method introduced by [START_REF] Saito | Computational bodybuilding: Anatomically-based modeling of human bodies[END_REF]] achieves nearinteractive run time skin, skeleton and soft tissue 3D model editing.

Fig. 2 summerizes this comparison between state of the art anatomy registration methods.

Non pose dependent

Automatic method

Soft Tissue

Skin and Skeleton

Figure 2: Comparison between state of the art anatomy registration methods. Legend : green means that the characteristic is totally handled by the method, red that it is not, and orange that it is partly handled.

User Tracking

In [START_REF] Pfister | Comparative abilities of microsoft kinect and vicon 3d motion capture for gait analysis[END_REF]] the authors assess that the basic Kinect body tracking results are enough for basic motion measurements such as stride timing, joint angles during motion, but are far beyond Vicon cameras in terms of software and hardware.

The tracking algorithm used in this paper is based on the Kinect SDK animated skeleton which is really noisy. Whereas we add constraints to upgrade the tracking, [START_REF] Meng | Kinect for interactive ar anatomy learning[END_REF]] ask the user to pinpoint anatomical key points to help positioning the data. [START_REF] Shen | Exemplar-based human action pose correction and tagging[END_REF]] use an example-based method to learn how to correct initially tracked poses.

Because body tracking is a critical step, other methods like [START_REF] Zhou | Posture reconstruction using kinect with a probabilistic model[END_REF] or [START_REF] Wei | Accurate realtime fullbody motion capture using a single depth camera[END_REF]] use the Kinect depth map and implement their own posture registration process using probabilities or pose estimations. [START_REF] Zhu | Adaptable anatomical models for realistic bone motion reconstruction[END_REF]] use multi-Kinect depth maps and anatomical knowledge to enhance realistic limb motions.

Nowadays, in the game industry, sports and fitness training applications using depth map tracking devices are commonly used (eg. Nike Kinect+, Get Fit With Mel B, Your Shape, etc...). To our knowledge, the best tracking games are based on the Kinect technology. All of these games shows the user depth map or silhouette only.

With AR, a precision and a realism constraint are added compared to this field state of the art by presenting the anatomy superimposed onto the user's color map. Fig. 3 highlights Kinect tracking problems such as disconnected bones head and overlaps between bones. AR Systems In the last few years, the number of AR applications increased in the medical education field [START_REF] Kamphuis | Augmented reality in medical education?[END_REF]].

The Magic Mirror [START_REF] Blum | mirracle: An augmented reality magic mirror system for anatomy education[END_REF] superimposes statically CT scans of the abdomen onto the user's image. The Digital Mirror [Maitre 2014] shows full body CT scans but does not superimpose them on the user image. In these two cases, data follow the user's motion but are not deformed with respect to these motions.

The Anatomical Mirror [START_REF] Borner | Anatomie spiegel[END_REF] allows fullbody motion by using the Kinect tracking, but it displays animated generic 3D models while we show a user-specific one.

Thanks to the use of anatomical knowledge, we significantly improve AR realism and anatomy motion plausibility with respect to [START_REF] Bauer | Interactive visualization of muscle activity during limb movements: Towards enhanced anatomy learning[END_REF]] and [START_REF] Bauer | Living book of anatomy project : See your insides in motion! In SA'15 Emerging Technologies[END_REF] in the Living Book of Anatomy project. Data Validation Validation of anatomical data requires in-vivo measurements, the simplest way is to use as ground truth body measurements [START_REF] Dao | Fast and automatic body circular measurement based on a single kinect[END_REF] or/and anatomical landmarks [START_REF] Espitia-Contreras | Development of a kinect-based anthropometric measurement application[END_REF]] taken directly onto the user's body. The study made by [START_REF] Malinowski | On using the microsoft kinect sensors to determine the lengths of the arm and leg bones of a human subject in motion[END_REF] gives limb bones length during motion and compare them with ground truth body measurements.

Using user body anatomical landmarks introduces measurement errors due to body position and skin curvature. We decided to use MRI data as ground truth: in addition to externally visible specific anatomical points to be able to obtain internal specific points (eg. femoral head of bone).

User-Specific Anatomy

We present a novel approach using Kinect SDK outputs (color map, body tracking skeleton and point cloud) and a 3D reference model including skin surface and internal anatomy (skeleton, muscles, organs, etc) to generate user-specific anatomical data.

The method consists of four steps. First, the user-specific body segment lengths and widths are computed using the Kinect SDK outputs (see Section 2.1) to define a list of 3D key points. In the second step the generic skin is deformed based on key points and the partial user's point cloud (Section 2.2). The third step consists in transfering the reference skeleton inside the user-specific skin (Section 2.3). Finally the soft tissue between the bones and the skin is determinded using Laplacian interpolation in a way similar to [START_REF] Dicko | Anatomy transfer[END_REF]]. These different steps are summarized in Fig. 5. Step 2: skin registration skin (j)

Step 3: skeleton registration

Step To ease the understanding of the rest of this section, descriptions of each type of deformation skeleton used are provided below:

• Kinect body tracking skeleton: composed of 25 joints, this animation skeleton is given by the Kinect SDK.

• skin registration skeleton: composed of 22 control frames (Fig. 8, red dots) defined on the generic 3d skin and corresponding to some of the Kinect body tracking joints, and 18 control points defined on the generic 3d skin contour (Fig. 8, green dots) that corresponds to the silhouette key points computed in Section 2.2. This system is used for skin registration.

• Internal anatomy registration skeleton: composed of 96 joint contraints between bones and 373 control frames position (see Fig. 11). This system is used to keep anatomical consistency during internal skeleton registration and is defined based on anatomical rules presented in Section 2.3.

Body size

The Kinect SDK provides a simple body tracking skeleton, without temporal coherence: links in-between segments may have different lengths at each frame. At calibration time: starting from a T-pose (Fig. 6.a), the user flexes his or her elbows (Fig. 6.b) and knees (Fig. 6.c). This allows us to estimate the lengths of upper limb and lower limb segments. The user silhouette and the body tracking skeleton given by Kinect are needed to compute body measurements (see Fig. 7.b) and define the 18 key points used for skin registration, as presented in Section 2.2.

The Kinect body tracking skeleton is mapped from camera space to image space using Kinect SDK tools.

A key point corresponds to the intersection between the user's silhouette edge pixel and a perpendicular line computed using a Bresenham algorithm. For robustness, we have designed a silhouette detection criterion: an edge pixel is defined by a black pixel followed by three white pixels to avoid silhouette holes.

For each key point, the Bresenham algorithm is initialized using the middle of in-between link segments as starting point and the perpendicular vector as the direction to follow. For instance using the point in-between shoulder and elbow link gives us upper arm width.

The 2d key points found are mapped from image space to camera space using Kinect SDK tools, Fig. 7.c shows the key points we use (eg. body tracking joints points, silhouette head points, silhouette waist points, etc. ...). Due to clothing and occlusion, some dimensions might be unreliable, especially thigh widths. Firstly, by assuming the human body symmetric along the sagittal plane, small errors in limb lengths are avoided. For each limb the average length value is used as real length in both sides. Other key point positions are infered based on the user silhouette and basic anatomical knowledge. Based on an average human body, we defined ratios between body parts. For instance, knowing that the thigh measurement should be half of the hip measurement, the thigh width can be infered. Some validations are shown in Section 4.

Skin registration

The skin registration method is based on the silhouette key points computed in Section 2. rigged using frame-based elastic deformers [START_REF] Gilles | Frame-based elastic models[END_REF]] corresponding to the Kinect body tracking skeleton joints (red dots in Fig. 8). Each skin vertex is controled by several frames, using linear blend skinning. The skinning weights are computed using Voronoi shape functions as in [START_REF] Faure | Sparse Meshless Models of Complex Deformable Solids[END_REF]]. The silhouette key points (green dots in Fig. 8) are mapped onto the skin to optimize the final result.

Instead of using global affine transformations (12DOFs) as in [START_REF] Dicko | Anatomy transfer[END_REF]; we use 9DOFs scalable rigids as frames, each bone matrix combines 3 translation, 3 rotation and 3 scale. The advantage over affine control frames is obtaining a better nonuniform local scaling to avoid shearing artefacts.

The skin model is registered to the target by minimizing a weighted sum of three energies [START_REF] Gilles | Frame-based elastic models[END_REF][START_REF] Gilles | Frame-based interactive simulation of complex deformable objects[END_REF]] using an implicit solver.

The predominant energy E skeleton attracts the control frames of the template to the bones of the user-specific model (red points in Fig. 8). Then the energy E keypoint attracts the silhouette points (green points in Fig. 8). Minimizing these two first energies scales the limbs, the torso, the neck and the head of the generic model according to the target body measurements as illustrated in Fig. 9.a and b.

The energy E cloudpoint attracts the skin to the target point cloud using an ICP apporach to define the correspondance. The forces are propagated from the skin vertices to the degrees of freedom: skeleton control frames (see Fig. 8 and Fig. 9.c).

Thanks to the fact that a small set of control frames are used, awkward configurations are avoided and no smoothness or kinematic constraint terms are needed.

Internal Anatomy Registration

User-specific anatomy reconstruction is divided in two sub-parts: anatomical skeleton registration and soft tissue registration. Soft tissue are deformed as described in [START_REF] Dicko | Anatomy transfer[END_REF]]; here the only focus is on internal skeleton registration. Inputs are the 3D reference of skin and skeleton model and the estimate of the user skin registered obtained in Section 2.2.

First, our method uses a volumetric interpolation to estimate the user anatomical skeleton. As in [START_REF] Dicko | Anatomy transfer[END_REF]], the use of Laplacian interpolation (Fig. 12.a) with as boundary condition the transformations between the two skins ensures that all the internal anatomy is bounded inside the user's skin after transfer.

A major limitation of the Anatomy Transfer [START_REF] Dicko | Anatomy transfer[END_REF] is the fact that the joint structure of the generic model is not maintained. Nothing prevents a bone from passing through another one (Fig. 12.b) or from being disconnected from a bone to which it should be connected (for instance ribs and thoracic vertebra, or ulna and humerus around the elbow joint, see Fig. 12.c). To keep correct joint structures and avoid these issues, joint constraints between the elements of our elastic bone model are added. The joint location, kinematics and limits are set according to [START_REF] Nordin | Basic biomechanics of the musculoskeletal system[END_REF]] (see Fig. 11). Thus, the internal anatomy registration skeleton is defined using frame based elastic deformations [START_REF] Gilles | Creating and animating subject-specific anatomical models[END_REF]] with weights computed using a Voronoi shape function as in [START_REF] Faure | Sparse Meshless Models of Complex Deformable Solids[END_REF] to smoothly propagate all long the bone each control frame transformation. 9DOFs scalable rigids for the control are used to keep head bone consistency as it is in the generic model. This guarantees that the bone heads can only translate, rotate and scale, and thus they keep a similar type of shape as in the generic bone model. The list of anatomical rules used to define the internal anatomy registration skeleton follows:

• R01: Keep long bones straightness (no bending or twisting)

• R02: Keep 3D model consistency: the complete set of entities is transferred to avoid holes To avoid bending bones (Fig. 12.d), an alignment constraint is added between the two bone heads. This constraint restrains the possible displacements between the control frames in only one direction defined by the line between them (see yellow lines in Fig. 11).

Thereby, the control frames can translate in one direction, but can still scale in all three directions. This alignment constraint is applied to long bones only.

It has been shown in [START_REF] Zhu | Adaptable anatomical models for realistic bone motion reconstruction[END_REF]] and similar approaches has been explored in [START_REF] Saito | Computational bodybuilding: Anatomically-based modeling of human bodies[END_REF]] that non-uniform scaling can be used to get more plausible bone deformations. This is why we introduced more control frames per anatomical bone. The number of frames varies according to bone type, the goal being to give enough deformability to each (for the registration process) while keeping good computation times (see blue dots in Fig. 11). For the short bones such as carpal bones, one frame per bone is used. For the long bones such as the femur two frames per bone are needed: one at the center of each bone head. For the flat bones such as the ribs three frames per bone are defined to keep ribs close to the skin in terms of curvature: two on bone heads (eg. close to the joints ribvertebra and rib-sternum), and one between the two others (middle of the rib). For bones with more complex shape such as vertebrae three frames per bone allows enough deformability to register the model while avoiding overlaps (e.g overlaps between facet joints, and spinous process of two different vertebrae). The complexity of the skull deserves a special treatment: use of five control frames for the whole skull deformation.

User Tracking

A single Kinect (markerless depth sensor) is used to perform body tracking. To reduce tracking noise, we record Kinect data in daylight, Kinect gives better results with background and ground matte materials. We observed that if the user's ground reflection is too visible, the Kinect includes it as part of the user silhouette which leads to lower limb length errors. The Kinect position is 60cm off ground for good lower-limb tracking results as determined in [START_REF] Pfister | Comparative abilities of microsoft kinect and vicon 3d motion capture for gait analysis[END_REF].

Because Kinect segments the depth map to compute body tracking joints at each frame, the in-between link distances change from frame to frame. This leads to disconnected articulations anatomical skeleton (on the limbs) or elongated meshes (on the torso zone).

We present the pipeline of our enhanced body tracking system in Fig. 13. Firstly we define a hierarchical body tracking system by constraining the limb lengths and by recomputing joint orientations (see Section 3.1 for more details). To smooth out small tracking noise, we apply a Kalman filter onto the joint positions. Joint orientations are recomputed from the filtered joint positions. Then we anatomically constrain the joint orientations: more details are given in Section 3.2.

Hierarchical body tracking system

Our hierarchical body tracking system is composed of 25 joints according to the Kinect SDK body tracking system. To define each joint f , the position and the orientation of its par- ent p is required. To overcome this, we begin by computing the joints from the root (spine base joint) to the leaves (eg. hand tips, foot joints and head joint). The root joint is defined by keeping the filtered Kinect position and orientation. The advantage of using a hierarchical skeleton is to obtain the body pose at each time t using only the joint rotations. We use the current Kinect body tracking skeleton to retrieve these rotations. Most of the time, orientations given by Kinect are incorrect so we decided to recompute them using in-between link directions by finding the smallest rotation R between initial direction (f c(t0)) and current direction (f c(t)), see Fig. 14.b. The 3x3 rotation matrix R is defined by the rotation of angle α around axis. For more details, see equation 1 and 2.

(a) (b) (c)

p(t 0 ) p(t) p(t) f(t 0 ) f(t) f(t) c(t 0 ) c(t) c(t) R fc(t0) fc(t)
α = asin -→ f c(t0) ∧ -→ f c(t) (1) axis = -→ f c(t0) ∧ -→ f c(t) -→ f c(t0) ∧ -→ f c(t)
(2) Fig. 14.c shows our hierarchical body skeleton system at step t.

Anatomically contrained joint orientations

To correct non-anatomically plausible behaviors due to tracking errors, each Kinect hierarchical body tracking joint orientation is constrained by limiting the number of possible rotations based on anatomical motion knowledge (eg. knee joint can be aproximated as a 1DOF joint, whereas the hip joint is a 3DOFs joint). This is done by constraining a given quaternion using Euler-angle constraints to find the closest rotation matrix defined only with valid axis within the joint limits. Computation is made using the Geometric Tools library [START_REF] Eberly | Constrained quaternions using euler angles[END_REF]]. Fig. 15.a illustrates in red a raw Kinect tracking and in grey the result after applying this constraint. To add even more anatomical plausibility to the result, joint limits are added to each rotation axis.

Results and Validation

To our knowledge, dealing with realistic anatomy visualization and motion is one of the most complex AR system ever because superimposing 3D anatomical data onto the user's color map reveals all the user measurement and tracking errors.

Our calibration method is a little time consuming (1-2sec for skin registration,15-30sec for skeleton registration and 30-60sec for soft tissue registration) but allows us to obtain a 3D model with accurate user measurements; moreover the motion capture pipeline, even with the introduction of delay during quick motions, leads to realistic and stable user tracking.

Thanks to these two features, the presented method allows a realistic experience for understanding anatomy. The described method is implemented in C++ and runs on a commodity laptop (Intel CoreI7 processor at 3 GHz, Nvidia Quadro K2100M and 8GB of RAM).

The real-time AR visualization runs between 35 to 62 fps depending on the 3D feedback: full-body musculoskeletal system (49211 vertices, 95189 faces) will run at 35 fps whereas internal organs (20144 vertices, 39491 faces) will run at 62 fps. The computational bottleneck of our system is the quality of the 3D model (number of faces and vertices) alongside the quality of the user color map (Kinect gives a high definition color map, which is reloaded at each frame).

The visual feedback can be performed on a commodity laptop screen or can be projected onto a 1.50m/2.0m screen for a demo display (see right side of Fig. 16).

Fig. 16 presents snapshots of the provided visualization. In a first set of experiments, the motion sequences were acquired for 4 men with an average height of 1.70m, and 3 women with an average height of 1.60m. To get uniform results we work with Kinect sequences made in similar environment conditions (daylight, background material reflections, Kinect position, etc...).

Fig. 17 presents two tracking data of the same user wearing different clothing and with different hair styles. It can be seen on the right side that the registered skeleton for these two datasets are almost identical; the red one is a little bigger (1.2% in the limbs lengths and 2.5% in torso widths) than the other one (green). This comparison allows the validation of our skin registration process (see Section 2.2).

Validation with MRI

The major contribution of our work, and also the most critical point is the closeness between the user-specific anatomy generated and the user's own. As explained in Section 1, using MRI data as ground truth allows us to obtain external as well as internal specific anatomical points for validation purpose.

Fig. 18 present MRI data of two users in front and lateral views side to side with the corresponding 3D user-specific registered anatomies. The internal anatomy registration skeleton introduced in Section 2.3 is used to set the 3D model in a similar pose as user's in MRI data. After comparing body height, we found an average error of 1.5% between 3D and MRI data which is quite accurate. 3D data being always smaller than MRI data, this error is due to limited skull deformations. In Kinect body tracking skeleton data, we observe a lot of change in limbs lengths. Thus, we pinpoint anatomical specific points (long bones protuberences) onto the MRI and onto the user-specific associated 3D model. With these specific points, we compare ulna and tibia lengths in real and 3D data. We suffer an average 5.2% error in limb lengths, most of the time the user-specific 3D model lacks a few centimeters. This percentage seems quite acceptable taking into account Kinect raw data noisiness.

To evaluate the torso body part realism, we propose to compare the user-specific 3D model and full-body MRI data by comparing the distance between left and right humerus bone heads. The average error between 3D and MRI data is rather small: 1.5%.

Fig. 10 shows that the point cloud and skin are fairly close; our generic skin being registered for a woman or a man, but what about internal anatomy? We know that women hips are in average larger than men to allow birth. The 3.1% error between MRI and 3D data (the 3D data being always bigger than the MRI data) demon-strates that the distances between left and right femoral bones head difference between women and men is well transcribed in internal anatomy.

Using the lateral view, we pinpoint specific points to find rib cage depth. The user-specific rib cage is always bigger than the MRI data (around 20% bigger). It may be due to the difference in posture during acquisition (for MRI data, the user is lying whereas for Kinect data the user is standing). It may also be due to the use of a partial point cloud instead of a complete one. Due to front view capture, we observe depth errors in the skull as well: the skull is about 12% bigger in depth in 3D than in MRI data.

Validation with User Study

In a second set of experiments, are involved 20 different subjects with no motor problems and working everyday on tools involving medicine or medical imaging. For each subject, we captured a range of full body motions involving upper and lower limb motion as well as torso motion.

The group is composed of 13 men between 24 and 54 years old (average height: 181cm, average weight: 82.6kg), and 7 women between 22 and 44 years old (average height: 164cm, average weight: 61.7kg This user study was designed to evaluate the believability of our system. We defined eight criteria which evaluation is given in table 2 to evaluate the quality of our mirror-like AR system.

Body position range (criterion C01) corresponds to motions while standing, crouching or sitting. In most cases, the results are well received. For other cases, limitations are directly connected to Kinect occlusion limitations.

Body orientation range (criterion C02) corresponds to body orientation from Kinect point of view: eg. facing, profile, 3/4, back.

When Kinect raw data are occluded or self-occluded, our system returns false motion poses: the more occlusion in Kinect raw data, the more errors we will have. A major topic of future work is to be able to handle important occlusion zones.

Motion range (criterion C03) defines simple motions like Flexion/extension of the knee, as well as complex motions in the extremities like finger motion or supination/pronation of the arm. We obtain high motion quality for simple motions, for complex motions we are limited by Kinect: this criterion is still in need of improvements. The Kinect SDK outputs a small number of joints which limits the body motion possibilities (eg. spine bending). Head Tracking could be improved by using Kinect facial tracking.

For Motion fluidity and delay (criterion C04) and Motion consistency (criterion C05), the goal is reached. Motion consistency refers to the absence of outliers during motion. We should state the fact that part of the visual latency that might occur comes from the low frame rate of the color map display.

Motion plausibility (criterion C06) corresponds essentially to joint DOFs and angular limits. For this criterion we obtain different results depending on the body segment studied. For instance, it is more easy to implement constraint for 1DOF joints than for 3DOFs joints such as spine or shoulders joints due to motion range.

Anatomy realism (criterion C07) gives a feedback on the registration method by focusing on limb length and torso width. For this criterion, people with professional knowledge in anatomy were the only ones to rate the user-specific anatomy as average.

For almost everyone, the Augmented reality (criterion C08) results were of good level. The overall quality can even be increased with mesh texturing. Table 2: User study complied results according to the quality criteria for a mirror-like augmented reality system. For each criterion: number of user having a bad/average/good evaluation of the criterion.

Conclusion

We present the first live system of personalized anatomy in motion. Superimposing the anatomy onto the user's image allows us to create a real-time augmented reality experience. The attached video (see https://youtu.be/Ip17-Vaqqos) illustrates the application pipeline and shows AR results of our system. We believe that the basic Kinect body tracking enhanced with our method is sufficiently accurate for our needs.

The system could be extended in different ways. Most users claimed in the user study that the overall quality of AR is of good level. However, some artefacts are still visible during motion and future work will be done to ensure that the 3D user-specific data always lie within the user's silhouette; we could apply a silhouette retargeting as in [START_REF] Zhou | Parametric reshaping of human bodies in images[END_REF]] to correct our hierarchical body tracking system. The addition of biomechanical simulations could allow to get more realistic deformations of soft tissues and organs but this could be at the cost of interactivity.

To show full body muscular activity for every possible body motion, inverse dynamics [START_REF] Murai | Musculoskeletal-see-through mirror: Computational modeling and algorithm for whole-body muscle activity visualization in real time[END_REF] will also be developed.

An improvement in the skin registration can be done by reducing the 9DOFs controllers to 6DOFs (3 rotations and 3 scales). This can be done by exploiting appropriately the hierarchical structure and would allow more robustness and skin consistency around body joints.

Our work is designed to be used as a tool for anatomy learning for medical and sports students. This is why in the future it is planned to display anatomical educational content (text, images, videos, etc) in addition to the AR visualization. This system could also be used as a way to communicate between medical practionners and their patients, about surgery, rehabilitation or any other health issue. Novel artistic content might also be produced using our technology, as well as interactive advertising.

Our system has been featured as a live demo during two conferences and at the Consumer Electronic Show. More than 400 people have been able to test it out. Most of them enjoyed the experience and a lot of them recommended it and came back with others.

Figure 3 :

 3 Figure 3: Kinect body tracking system and user-specific anatomy.

Fig. 4

 4 Fig.4 summerizes comparisons between state-of-the-art demos and our work.

Figure 4 :

 4 Figure 4: Comparison between anatomical mirror-like applications. Legend : green for good, orange for average, red for bad.

Figure 5 :

 5 Figure 5: Pipeline of the user-specific anatomy generation.

Figure 6 :

 6 Figure 6: Calibration process: (a) to have global proportions (head and torso), (b) for real upper limb parts lengths, (c) for real lower limb parts lengths.

Figure 7 :

 7 Figure 7: (a): skeleton key points. (b): body measurements key points. (c): 3D key points used in skin registration.

Figure 8 :

 8 Figure 8: Skin registration. Red dots: origins of control frames; green dots: silhouette key points; blue dots: Kinect point cloud.

Figure 9 :

 9 Figure 9: Skin registration result at the end of each step of the optimization process. (a): minimizing E skeleton . (b): minimizing E skeleton and E keypoint . (c): minimizing the three energies.

Fig

  Fig.10presents the skin results after registration with the corresponding Kinect point cloud. By using E cloudpoint , the torso skin is slightly deformed to refine the model: in the same way, the user being a woman or a man.

Figure 10 :

 10 Figure 10: Kinect point cloud and corresponding registered skin. Top: 1.55m female. Bottom: 1.85m male.

Figure 11 :

 11 Figure 11: Right arm internal anatomy registration skeleton. Blue dots for control frame positions; yellow lines and middle bone frames for alignment constraints; other frames for joint contraints.

Figure 12 :

 12 Figure 12: (a): Laplacian Interpolation, (b): registration without joint contraint (overlaps between bones), (c): registration without joint contraint (bone heads disconnected), (d): registration without alignment constraint (bent and twisted bones), (e): our method.

Figure 13 :

 13 Figure 13: Enhanced body tracking pipeline.

Figure 14 :

 14 Figure 14: (a): our hierarchical body tracking skeleton at (t0). (b): Kinect body tracking skeleton at (t). (c): our result.

Figure 15 :

 15 Figure 15: Kinect data in red and corrected in grey. (a): off angular limits rotation. (b): rotation axis error (DOFs).

Figure 16 :

 16 Figure 16: Left: system set-up. Right: snapshots of results.

Figure 17 :

 17 Figure 17: For the same user with different clothing and hair style (Left), we obtain almost identical results (Right).

Figure 18 :

 18 Figure18: Morphometric measurements (green lines) used to compare our results with ground truth MRI data. Top: 1.55m female. Bottom: 1.85m male.

Table 1 :

 1 ). Table1gives global informations about each subject. Global informations about each user study subject.

	Gender Age Body Mass Index Anatomy Knowledge
	male	44	32.4	poor knowledge
	male	47	29.3	professional
	male	31	27.4	professional
	male	41	26.3	poor knowledge
	female	44	26	average knowledge
	male	41	25.5	average knowledge
	male	24	25.1	poor knowledge
	female	28	24.4	professional
	male	30	24.3	average knowledge
	male	30	24.3	professional
	male	27	23.7	average knowledge
	male	54	23.6	average knowledge
	female	22	23.5	poor knowledge
	female	34	22.2	average knowledge
	female	32	22	professional
	female	24	21.8	poor knowledge
	male	29	21.6	average knowledge
	male	28	21.6	average knowledge
	female	31	21.3	poor knowledge
	male	31	20.5	average knowledge

Acknowledgments

Many thanks to Estelle Charleroy our Graphist Artist, and to Olivier Carre for engineering advices. We also would like to thanks all people involved in the user study, and all people helped testing the application. MRI data used for validation are from the A2H study promoted by Grenoble University Hospital. This work has been partially supported by the LabEx PERSYVAL-LAB (ANR-11-LABX-0025).