Beata Taudul

M Bowen

M Alouani

Impact of single and double oxygen vacancies on electronic transport in Fe/MgO/Fe magnetic tunnel junctions

Keywords: numbers: 71.15.Mb, 73.20.-r, 73.40.-c, 73.40.Qv

The combination of a low tunnelling barrier height and a large tunnelling magnetoresistance (TMR) ratio in MgO-class magnetic tunnel junctions (MTJs) has enabled next-generation information storage and bio-inspired computing solutions thanks to the spin transfer torque effect. Recent literature has proposed that this synergistic combination arises from the electronic properties of oxygen vacancies. To explicitly understand their impact on spin-polarized transport, we have computed the electronic and transport properties of single (F centers) and paired (M centers) oxygen vacancies using density functional theory and the projector augmented wave (PAW) method. These point defects can generate energy level positions of 0.4 eV with respect to the Fermi level for FeCo electrodes irrespective of the defect's spatial position within the MgO barrier, and of the orientation of the M center. These defects promote a strong decrease in the conductance of the spin up channel in the MTJ's parallel (P) magnetic state that mainly accounts for an order-of-magnitude drop in TMR, from ≈10000% in the ideal case toward values more in line with experiment. When placed in the middle layer of the MgO barrier, the F center introduces additional P ↑ transmission away from the Γ point. This scattering lowers TMR to 145%. In contrast, the M center merely broadens this transmission around Γ, thereby boosting TMR to 315%. Rotating a M center so as to partly point along the transmission direction sharpens transmission around Γ, further increasing TMR to 1423%. When these defects are placed at the MTJ interface, the transmission and ensuing TMR, which reaches ≈4000%, suggest that such junctions behave as would an ideal MTJ, only with a much lower TMR. Our results thus theoretically reconcile the concurrent observations of high TMR and low barrier heights, in line with experimental preparation techniques such as post-deposition oxidation of metallic Mg, which can generate oxygen vacancies at the lower MTJ interface, and annealing which can promote M centers over F centers. Our theory is also in line with an origin of perpendicular magnetic anisotropy in terms of oxygen vacancies at MTJ interfaces. The effective size of these vacancies sets a limit for both the barrier thickness, in line with experiment, as well as for the MTJ's lateral dimension. Our work provides a much-needed theoretical basis to move beyond the mostly unsuspected, fortuitous defect engineering of spintronic performance that has thus far propelled MgO-based spintronics and its applications.

I. INTRODUCTION

Spintronic research exploits both charge and spin degrees of freedom in solid-state systems [1][2][3][START_REF] Peng | Wiley Encyclopedia of Electrical and Electronics Engineering[END_REF]. A widely studied spintronic device is the magnetic tunnel junction (MTJ), composed of two ferromagnetic metallic electrodes separated by an ultrathin dielectric.

The electrical resistance of the MTJ depends on the relative orientation of the electrode magnetizations, which can be controlled by an external magnetic field or a spin-polarized current [START_REF] Locatelli | [END_REF]. This change in resistance is called tunnelling magnetoresistance (TMR) and is defined as

T M R = R AP -R P R P , (1)
where R P and R AP are the resistances for parallel (P) and antiparallel (AP) configurations of the two magnetizations of the electrodes.

Initial MTJs with an Al 2 O 3 tunnel barrier exhibited a maximum TMR of 70% at room temperature [6,7]. Theoretical studies of magnetotransport across Fe/ZnSe/Fe [8] and Fe/MgO/Fe [8][9][10][11][12] MTJs revealed that the TMR can be greatly increased if the amorphous barrier is replaced by a crystalline one, such that certain orbitals with a high spin polarization in the electrodes preferentially tunnel across the barrier. Nowadays, textured FeCoB/MgO/FeCoB MTJs with TMR values above 600% at room temperature [13] offer promising prospects for data read-out, storage and processing, magnetic sensors [14][15][16].

Despite the importance of these MTJ technologies, an understanding of exactly how the device operates remains a work in progress. Indeed, the TMR effect is a complex phenomenon that depends strongly on the electronic structure of the electrodes, the properties of the insulating barrier and on the chemical bonding at the MTJ interface. As an illustration, consider how the success of achieving high TMR concurrently with low barrier heights required [START_REF] Halisdemir | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF] to implement spin transfer torque toward these MTJ technologies [START_REF] Locatelli | [END_REF][14][15][16] imply that structural defects, which may lower TMR from the 10000% theoretical prediction, may actually play a beneficial spintronic role here.

Several causes for an effective deviation from the MTJ's ideal structure have been considered. Experiments often reveal the presence of interface oxidation, which alters the nature of chemical bonding at the interface between the ferromagnetic electrodes and the MgO spacer and can degrade TMR [START_REF] Meyerheim | [END_REF]. The combination of theoretical [19] and experimental [20] studies showed that, even if one includes the interface disorder or the oxidation of interfacial Fe layer, the drastic drop of TMR cannot be fully explained by this mechanism alone.

As another cause, atomic diffusion may occur during the sample preparation and annealing [21]. In particular, boron diffusion into the MgO barrier (forming boron oxides), or its segregation at the CoFe/MgO interface, has been studied [22][23][24][25], but this is not always the case [26][27][28]. Rather, at a proper annealing temperature, boron does not diffuse into MgO but rather goes further away from the interfaces. Even if boron diffuses into MgO, it was shown theoretically that this should not create additional states within the MgO band gap [29].

Finally, another source can be imperfections in the MgO spacer itself, such as grain boundaries and point defects. The impact of the grain boundaries on the electronic structure and on the transport is difficult to address. Nonetheless, it was shown by Mizuguchi et al. [30] that the tunnelling current flows uniformly despite the existence of the grain boundaries and hence the device performance is not affected considerably by this kinds of defects. Moreover, the combined experimental and theoretical investigations of Bean et al. [31] showed that grain boundaries can cause a decrease of the effective barrier of MgO but this band gap decrease can not explain the observed low barrier heights [32].

Point defects, on the other hand, can promote localized states within the band gap of MgO, giving rise to a variety of interesting optical, catalytic and transport properties that are absent in the ideal crystalline material [33]. The most plausible imperfections are oxygen and magnesium vacancies, denoted F and V centers respectively. They can appear in a neutral, singly charged or doubly charged state that is denoted as F + , F 2+ , V -V 2-, respectively. Moreover, two point defects can form a paired vacancy: two F centers form a F 2 pair of oxygen vacancies, which is a M center when they are nearest-neighbor on the oxygen sublattice. A F center can also combine with a V center to form a MgO vacancy.

As discussed by Gibson et al. [34], oxygen vacancies exhibit the lowest formation energy, which implies that this species of defects is more likely to occur in MgO. This defect species promotes localized states in the band gap of MgO and can affect the optical and the electrical properties of the dielectric [32,[34][35][36]. As a result, the barrier heights encountered by the propagating electrons are locally reduced [37,38]. The electrons can then tunnel through the barrier using scattering states at the Fermi level, modified by the defect states, which are different from the scattering states for an ideal barrier as it can be seen later in Fig. 9.

This would explain experimental reports of a barrier height in MgO MTJs that is far below the nominal value of 3.9 eV (see Tab. I).

MTJ TMR (%)

Barrier height (eV)

Fe/MgO/Fe 130 (190 1K ) 0.38/0.82 [39] Fe/MgO/Fe 180 (247 20K ) 0.39 [40] FeCo/MgO/FeCo 120-220 1.1-1.7 [41] FeCoB/MgO/FeCoB 100 0.62/0.5 [32] Fe/MgO/FeCo 23 4.2K / 20 70K 0.9 [42] Although oxygen vacancies within MgO appear to play an important role not only toward MTJ performance, but also spin transfer torque [START_REF] Halisdemir | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF], their clear identification [43,44] and impact on the tunnelling current has remained a work in progress. According to theory [35,36,45], single oxygen vacancies should create barrier heights of about 1.1 eV for the tunnelling electrons and decrease the resulting TMR. Even if we consider more Ftype vacancies within MgO [45], the general conclusion is that these vacancies should degrade the TMR ratio [39]. On the other hand, in the presence of a 0.4 eV barrier height, coherent transport seem to be spintronically favorable [32,40,46]. Indeed, it was suggested by Schleicher et al. [32] that this barrier arise from M centers, which preserve coherent transport according to McKenna and Blumberg [46]. Only in a recent work, we have shown an unified experimental/theoretical picture of the potential landscape due to F and M centers [47,48].

By studying the complex band structure of MgO with F and M centers we showed that the M centers reduce the barrier height in agreement with our experimental findings [47]. In a different study, we have shown that the energy positions of the F and M centers, with respect to the Fermi level, can be tuned by the choice of the type of the electrodes. In particular, this energy shift can be as large as 0.4 eV when the iron electrode is substituted by cobalt and was explained in terms of the change in the electrode work function [48].

In this work, we extend our previous studies [47,48] by analyzing the electronic properties of F and M centers in MgO, and their impact on spin and symmetry polarized transport in Fe/MgO/Fe junctions employing density functional theory. The transport is calculated within the Landauer-Büttiker formalism as implemented in PWcond [49]. We show that the position of defect levels with respect to the Fermi level of the MTJs is robust against the type of exchange and correlation functional, and do not depend on the defect's spatial position in the barrier. Our results indicate that M centers can account for the experimental barrier height of 0.4 eV, and promote improved TMR relative to F centers. We also find that the defect's position within the barrier thickness strongly conditions spintronic performance.

Relative to an ideal MTJ, this performance is mostly unaffected when the F or M center is close to an interface, but is reduced when the vacancy is moved onto the barrier's middle monolayer. In that case, rotating the M center restores stronger TMR. When judicious, we discuss how our theoretical framework of spintronic tunnelling across MgO in the presence of F and M centers adheres to experiment. The spatial extent of the F and M center suggest that a MTJ with a lateral size ≈ 2 nm can still exhibit high TMR [47].

The paper is organized as follows: In section II we present the details of our calculations and methodology. In sec. III we compute and discuss the electronic ground state properties of F and M centers either in bulk MgO or incorporated into MTJs. We compare the defect level positions obtained theoretically with the experimental data and discuss its change with respect to the types of electrodes. In sec. IV we show the results of transmission calculations for Fe/MgO/Fe MTJ with oxygen vacancies generated within MgO spacer. We explain the importance of the geometrical position of the vacancy with respect to the interface and the orientation of the defect plane for the M center. In the last section, we conclude the paper with a general discussion.

II. METHODOLOGY

A. Ground state calculations

To calculate the electronic structure of MgO with F and M centers, denoted as F(M)-MgO, we created the F/M centers by removing one/two neutral oxygen atoms from a simple cubic supercell containing 64 atoms. We used the experimental MgO lattice constant of a MgO = 4.21 Å. These calculations were performed using VASP package [50,51] based on the projector augmented wave (PAW) method [52] and the Perdew, Burke, Enzerhof (PBE) generalized gradient approximation [53] for the exchange-correlation potential. The kinetic energy cutoff value of 500 eV for the plane wave basis set and the convergence criterion for the total energy of 1µeV is used. The structures with defects were relaxed by requiring that the forces acting on atoms be less than 0.001 eV/Å. Due to the large size of the supercell, we found that a k-point mesh of 4×4×4 using the Methfessel-Paxton method with a smearing of τ = 0.2 eV yields a satisfactory convergence of localized states resulting from these defects. Since it is well known that the generalized gradient approximation (GGA) underestimates the size of the band gap, we also used the hybrid Heyd, Scuseria, Ernzerhof (HSE) functional [54,55] to correct the band gap and verify whether the defect level positions with respect to the Fermi level depends on the type of functional employed.

The HSE hybrid functional mixes a portion of exact Fock exchange with that of DFT using an adjustable parameter µ. We found that, by increasing the Fock exchange to 43% in HSE06, we can reproduce the experimental band gap of MgO. This parameter is then used to calculate the defect levels in the MgO MTJs employing HSE06.

To determine the positions of the defect levels with respect to the Fermi energy and compare them with experiment, we used more realistic Fe/MgO/Fe MTJs. At the metal/insulator interface, the metallic electrode and the electronic transfer between the two materials and the metal induced gap states (MIGS) in the band gap of MgO will set the Fermi level position for the junction and establish the energy position of the defect levels accordingly. The geometry of the MgO/Fe supercell was based on the experimental results: the Fe conventional unit cell was rotated by a 45 • with respect to that of MgO to match the lattice constants of both materials and avoid strains in the structure: a MgO = √ 2a Fe . In addition, oxygen atoms were placed on top of Fe atoms and the Fe-O distance at the interface was fixed to 2.17 Å following previous theoretical predictions [35,56]. It is important to notice that the measured Fe-O distance is in the range of 2 Å[57, 58] to 2.2 Å [START_REF] Meyerheim | [END_REF]. We fixed the lattice constant of MgO to its experimental value (a MgO = 4.21 Å) and adjusted the lateral lattice parameter of the electrodes to it. This choice reflects the experimental evidence [59] that the annealing of FeCoB/MgO-based MTJs led to a recrystallization of the electrode/barrier interfaces so as to adopt the MgO lattice constant. The lattice parameter along the z axis was rescaled accordingly. The structure of the junction and relevant parameters are indicated in Fig. 1.

Note because of the high CPU cost of HSE functional, the relaxation of the atoms around the F and M vacancies were conducted using only the PBE functional and no significant atomic motion was found in agreement with previous calculation of oxygen vacancies in MgO [68]. We have also found that the the lattice relaxation of the MgO supercell did not change significantly the DOS with or without oxygen vacancies (not shown), however this lattice relaxation does not correspond to the real situation where the MgO slab is inserted between two semi-infinite metallic electrodes. Our calculation uses periodically repeated Fe/MgO/Fe slabs. To relax the geometry along the z direction one needs to study very large supercells as the inter-distance of many Fe and MgO layers near the interface will vary as it is shown in our previous study [60]. The full relaxation of the MTJs is CPU intensive a very cumbersome and will just complicate the study of the effect of F and M centers on the transport properties. As the aim of the present work is to compare the transport properties of F and M centers and show that the M centers are important for achieving higher TMR.

We believe that this physics will not be significantly affected by the structural relaxation especially because it is found that the atomic relaxations around the F and M centers are negligible. The experimental FeCoB electrodes are initially amorphous alloys whose interface with MgO adopts the latter's (001) texture upon annealing [59]. However, the exact arrangement of Fe and Co atoms is not really known. Moreover, the role of boron or its influence on the electronic structure of FeCo electrode is not clear. Therefore in our work we decided to consider a simplified scenario in which the FeCo electrode has the same structure as that of Fe. For the FeCo/MgO calculations, we have alternated the Fe and Co layers and chosen

Co at the interface with MgO. Nonetheless, even the simplified structure can give us some insights into the impact of the type of electrode on the defect levels.

For the calculations involving full junctions, the kinetic energy cutoff value and the convergence criterion were the same as for the bulk MgO, but the value of τ was decreased to 0.1 eV and the spin polarized version of PBE-GGA was used. To compute the electronic structure of oxygen vacancies, the lateral directions of MTJ were doubled and periodic boundary conditions applied in all directions. The atoms around the oxygen vacancies were allowed to relax as described above.

B. Transport calculations

The ballistic conductance is calculated using Landauer-Büttiker [61,62] formula

G(E F ) = e 2 h n,k,σ T n,k,σ (E F ), (2)
where T n,k,σ (E F ) is the total transmission at the Fermi energy and the summation is over all bands n crossing the E F , for each k point and spin σ. Note that since the spin-orbit coupling is very small in the Fe/MgO/Fe MTJs and we are concerned only on the conductance at a very small bias, we ignored the spin mixing processes. The electron transmission was evaluated using the scattering based approach with a plane wave basis set and ultrasoft pseudopotential (USPP) scheme as implemented in the PWcond [49] module of the Quantum Espresso (QE) [63] package. Note since the ballistic conductance is not yet implemented in VASP, we had to use QE to compute it. However, we had to make sure that both code generated the same electronic structure for MTJ electronic structure.

The conductance was evaluated between two semi-infinite electrodes connected by a scattering region that contains an insulating MgO spacer and a part of the leads on each side of the spacer. To compute the electron transmission at a given energy E, we first calculated the total energy of the ground state properties with the PWscf code from QE package [63] and determined the effective potential. We then constructed the generalized Bloch states, including propagating and evanescent states, as a solution of Kohn-Sham equations at energy E for the infinite periodic leads, and the results were used to construct the scattering states and compute the transmission across the entire system. Moreover, in the spin density functional picture, electrons of different spin move independently in their different selfconsistent potentials [49]. Therefore in this approach spin flip events are not included and the total transmission is the sum of the two spin channels such that

T (E) = T ↑ (E) + T ↓ (E).
To find how many electrode layers should be contained in the scattering region, we studied the changes of the electrostatic potential in the scattering region. To ensure that the electron wave function changes smoothly at the interface between the bulk of the electrode and the scattering region, the part of the leads in the scattering region has to be big enough so that the changes induced in the electrostatic potential due to the interaction with MgO are contained entirely within the scattering region. If not, an artificial potential that scatters the incoming electrons could be present and might affect the results. By comparing the total electrostatic potential for the scattering region to the total potential of the bulk electrode, we found that from the 2nd-3rd monolayer (ML) of Fe, the bulk electrostatic potential is restored. In order to guarantee a proper geometrical matching between the scattering region and the electrodes, we used 4 ML of Fe on the left side of MgO and 5 ML of Fe on the right side. For the defect calculations, we doubled the lateral size of the junctions. When the antiparallel alignment of the electrodes was considered, the size of the junction along the z direction was doubled, such that the composition of the supercell was Fe(P)/MgO/Fe(AP)/MgO/F(P).

In the ground state calculations with PWscf code, the cutoff energy values for the plane wave basis set and the electron density were set to 40 Ry and 400 Ry, respectively. The electronic occupations were broadened using a Gaussian smearing technique with a smearing parameter τ = 0.02 Ry. The total energy convergence threshold was set to 10 -8 Ry and the electron density mixing parameter to 0.1. Since we needed to use the same form of pseudopotential and the corresponding exchange-correlation functional for all atoms in the junction, we chose the Perdew and Wang (PW91) generalized gradient functional [64] in a spin-polarized form already generated and available in the QE library. For the ground state calculations of the ferromagnetic alignment of the electrodes, we used a k-point mesh of 5×5×1, while for the antiferromagnetic alignment the same k-point grid was slightly shifted out of the Γ point in order to speed up the convergence. Indeed to calculate the antiferromagnetic state we doubled the unit cell. This results in almost folded band structure producing numerically degenerate states along the high symmetry directions of the Brillouin zone. To avoid working with the eigenvectors of degenerate bands we shifted the k-point grid by half a step in each direction and required that the grid preserves the crystal symmetry.

An important factor in the transmission calculations is the convergence of the 2D basis set used in the PWcond. Here, two parameters control the basis set: (i) ewind defines the energy window for reducing the 2D plane wave basis set in the transverse xy plane, and

(ii) epsproj is a threshold for the 2D basis set reduction. The default values for the two are ewind = 1 Ry and epsproj = 0.001, but in general, the larger the ewind, the smaller the epsproj, the higher the accuracy of the transmission. However, the increase in the transmission accuracy increases the computational cost and a suitable compromise should be found. These parameters were tested by examining the complex band structure (CBS) of bulk Fe and MgO [65,66]. We found that ewind = 3, epsroj = 10 -6 were sufficient to converge the CBS and hence they were used also to compute the transmission. In addition, the transmission was evaluated as a function of the number of k points in the 2D BZ.

We tested meshes of 20×20, 30×30, 50×50 and 80×80 k-points and kept the 50×50 which showed a well-converged transmission.

III. ELECTRONIC STRUCTURE A. MgO bulk

Here, we briefly describe the ground state electronic properties of M center within MgO bulk and MgO incorporated in the Fe(FeCo)/MgO MTJs. The removal of two neighboring oxygen atoms from an MgO supercell results in the creation of two occupied energy levels below E F . The electrons that were transferred from Mg and remain after the oxygen removal are trapped and are mostly localized on the vacancy sites. Atomic relaxation showed that the atoms around the defect are only slightly distorted and the electronic structure of the defect is essentially not affected. This is due to the fact that the electron distribution on the vacancies resembles that of an oxygen O 2-ion. The distortions around the defect are therefore neglected [67]. In addition, the generation of double vacancy creates states near the minimum of conduction bands, which we refer as the excited states of M center (see Ref. [47] and Ref. [67] for details). These defect states could produce important features in optical spectra due to their localization within the defect region.

To understand the nature of the M center levels, we plot in Fig. 2 We plot as insets to Fig. 2 the electronic spatial distribution for each of the M centers ground states. It is clear that the electron distribution for the M 1 state resembles a bonding-like state, while it is anti-bonding for the M 2 state. As in the case of bond formation between atoms, the coupling between two F centers creates a bonding state with a lower energy, and an anti-bonding state with higher energy, with respect to the original F state energy level.

Indeed, the F center peak is always positioned in between two M center ground states. As a consequence, the energy levels created by the F center is always higher than that associated with the M 1 state.

We emphasize that the nature of the MgO growth matters when comparing experiments with theoretical studies on nearly ideal MgO [68]. Annealing-induced control over the M/F ratio of centers is achieved only for dc-sputtered MgO, i.e., steps of Mg sputtering in Ar atmosphere plus oxidation to form the MgO; see Ref. [47] and not rf-sputtered MgO, i.e., one step of MgO sputtering in Ar atmosphere; see Ref. [69]. This is because dc-sputtered MgO has many more defects (including oxygen vacancies) than rf-sputtered MgO, in line with the resistance area (RA) products obtained (3 orders of magnitude lower for dc-MgO than rf-sputtered MgO). These experimental observations of a combined a) lack of control over the M/F ratio of centers and b) much higher RA products in rf-sputtered MgO MTJs are in line with the high diffusion barrier reported for the near pristine MgO in Ref. [68]. decay with the number of MgO layers. As a result the band gap of MgO disappears at the interface because of states coming from Fe. From the third layer the band gap of a bulk MgO is restored. Moreover, the ferromagnetic electrode induces spin polarization in the neighboring MgO layers and the difference in the DOS of spin up and down electrons can be clearly seen.

For the Fe/MgO heterostructures, the M 1 and M 2 states are located respectively at -1.7 eV and -0.7 eV below the Fermi level, while the F center level is at -1.2 eV. As expected the energy level associated with the F center is in between the M 2 and M 1 levels. In the case of the FeCo electrodes the defect levels are shifted towards the Fermi level by about 0.5 eV.

It is worth noticing that not only are defect levels shifted, but also the whole band structure of MgO is rigidly shifted towards higher energies. This shift can be understood considering the 0.5 eV difference in the work function between Fe and Co. As before, the vacancy affects also the closest MgO layers up to 3 ML of MgO along the direction perpendicular to the M-center plane. We also studied the effect of shifting the M center within the MgO (7 ML) spacer and found that the defect energy level remains practically unchanged when the vacancy approaches the interface (data not shown) [67]. However, if the vacancy is placed on the interfacial MgO layer the DOS associated with the F/M center is washed out due to the strong interaction with the ferromagnetic electrode. Clearly, when the defect is closer to the interface, the hybridization between the two types of materials is strongly affected, which in turn influences the position of the Fermi level. However, it is surprising that this modification led only to a small differences in the Fermi level positions of the order of 0.09 eV.

Thus far, we have considered only the situation where the M center is in the plane parallel to the interfaces. In that case, we found that the effective size of the M center in the direction perpendicular to the interface reaches up to 3 ML of MgO due to charge transfer onto adjacent MLs. We can also rotate the M center such that it can be partially aligned along the z direction and shared between two neighboring MgO layers. In that case, the effective M center size reaches 4 ML. Regardless of the orientation of the defect plane, the level position of the M center remains practically unchanged and similar results for the DOS are also obtained (not shown) [67]. Even though the changes in the DOS upon shifting or rotating the M center within the MgO spacer are not significant, we will show later that these changes have a huge impact on the transmission as hinted by complex band structure calculations [47]. Thus, based on our calculations, we can state that the computed energy range for energy levels associated with an M center can be associated with those measured experimentally at 0.4 eV, due to the M 2 state located ≈0.4 eV below the Fermi level for a FeCo/MgO/FeCo MTJ.

In the case of FeCo/MgO/FeCo MTJ, the interface termination will not affect the results substantially. What would change is the shape of DOS at the interfacial layer due to hybridization with Fe rather than Co. We have showed in Ref. [48] that the defect level positions can shift slightly for different electrodes; the cases of pure Fe and pure Co electrodes are the two extreme limits. In the case of CoFe/MgO termination, the defect levels should be slightly shifted towards lower energies with respect to FeCo/MgO termination.

The transmission for the spin up channel in parallel electrode configuration for all types of electrodes, (Fe, Co and FeCo/Co) has in general the same shape and similar intensity, bigger differences are found in the spin down channel and AP configuration. The last two are affected by interfacial resonant states and are more sensitive to the changes at the interface.

In this study we did not consider a random distribution of oxygen vacancies because of the prohibitive computational cost for computing the transport properties of large supercells. Here our aim is to compare the F and M center as a function of its position in the barrier region. We considered how, by controlling the position of a precise oxygen vacancy species, one may alter magnetotransport. This reasoning is valid in the limit of a single nanoconduction channel that dominates transport, and in the limit of laterally nanoscaled junctions (down to 2 nm). It will be useful in future work to consider different vacancy configurations and concentration and their effect on TMR using a tight-binding model as performed Tsymbal et al. [70] and Klyukin et al. [71].

We first calculated the transmission at the Fermi level for junctions with an ideal, 5 MLthick MgO structure in order to examine the impact of introducing oxygen vacancies. The results found for the 7 ML spacer are similar to these for 5 ML and will be only briefly discussed.

Fig. 6 presents the transmission in the two dimensional Brillouin Zone (2D BZ) for the parallel electrode magnetization for the spin up and the spin down electron channels (left and middle panels) and the corresponding transmission for the antiparallel configuration (right panel). In agreement with previous theoretical predictions [56,72], we found that the majority electron transmission is centered around the Γ point and dominated by the ∆ 1 symmetry. The transmission for the minority channel occurs basically at the edges of the 2D BZ and is much smaller than for the majority channel. The transmission in the AP configuration is a mixture of features seen in both spin channels. By summing the transmission over the BZ for each channel, we obtain the conductance (see eq. 2) and the resulting TMR. Tab. II summarizes the transmission results for junctions with 5 and 7 ML of MgO. As expected, the transmission decays exponentially with the thickness of the MgO spacer and hence drops by at least one order of magnitude when passing from 5 to 7 ML of MgO and the TMR increases with the number of MgO layers.

This reflects the favorable symmetry filtering across MgO(001) of ∆ 1 electrons with a high spin polarization at the bcc(001) Fe electrodes's Fermi level. For this reason, the ∆ 1 channel is blocked in the MTJ's AP magnetic state, such that transmission is ensuring by ∆ 5 and ∆ 2 symmetry channels, which appear for both spin populations at E F in Fe. When the MgO thickness is increased, the contributions to the conductance from the strongly attenuated ∆ 5

and ∆ 2 channels become smaller. This leads to a bigger overall difference in the transmission between the P and the AP configurations and causes the increase in TMR. According to literature [45], the TMR value should continue to grow at least up to 13 ML of MgO due to the dominant ∆ 1 contribution.

B. F/M center in the middle layer of MgO

In the next step, we introduced single and double oxygen vacancies in the middle layer of the MgO spacer. Fig. 7(c) and 8(c) show the corresponding 2D BZ transmission for F and M centers, respectively. The BZ transmission distribution for P↓ is only slightly affected by the presence of the vacancies, resulting in a small increase in the total amplitude with respect to the ideal case (see Table II). However, in the spin up channel, a clear distinction in the P↑ transmission between the F and the M centers can be made. It appears that the F center scatters the propagating electrons to states with higher k-vectors. As a result, the transmission has a minimum at the Γ point and occurs mostly along k x and k y directions with maxima at the edges of the 2D BZ. The electrons are scattered symmetrically in each direction due to spherical symmetry of a single oxygen vacancy. On the other hand, the P↑ transmission in the presence of the M center becomes broadened in the 2D BZ but maintains a symmetric maximum at the vicinity of the Γ point. This clearly suggests that, while the transmission across a F center is reduced by an order of magnitude due to transport across k = 0 states, coherent transport that preserves spin and symmetry of the electron wave function is still possible when M center is present in the MgO spacer.

We found that, while F and M centers promote a reduction in total P↑ transmission (see Tab. II), both centers promote an increase in both the P↓ and AP total transmissions.

Furthermore, the transmission distribution in the AP configuration changes significantly from that of an ideal MgO junction (see Fig. 6). While introducing defects reduces the TMR, the TMR is higher for M centers compared to F centers. We found similar trends for the 7 ML , (see Tab. II). Again, if we increase the number of MgO layers (here from 5 to 7 ML), the TMR also increases regardless of the defect type.

In the theoretical model structures, the increase of MgO thickness results in a decrease of the defect concentration. The concentrations for F and M centers are respectively 2.5% and 5% for MTJs with 5 ML MgO and 1.78% and 3.57% for the 7 ML MgO. Nonetheless, as we have seen in the cases of 64 and 216 atom supercells with M(F) center concentration of 6.25%(3.13%) and 1.85%(0.93%) respectively, the reduced size of the supercell did not affect significantly the defect electronic properties. In the case of MTJs with increasing MgO thickness, the resulting TMR in presence of defects increases (Tab. II). This indicates the manifestation of symmetry filtering effects and the decreased impact of the defect states.

However, the TMR for ideal MTJ is one order of magnitude larger than that of MTJ with oxygen vacancies and there is a clear distinction between F and M center even for smaller defect concentrations. Since the M center promotes higher TMR than the F centers in MgO MTJs with FeCo electrodes, these transmission results can account for the simultaneous experimental occurrence of high TMR alongside 0.4 eV barrier heights. They also confirm the initial assumption that coherent transport can be preserved when a M center is present.

Note that the defect level positions discussed previously were evaluated using the VASP code with the PAW basis set. To verify the robustness of these results, we switched to a plane wave basis set in conjunction with an ultrasoft pseudopotential approach. While the shape of the layer-projected DOS is practically the same, we noticed a small shift of about 0.15 eV of the F and M 1 states towards lower energies. We then examined how this shift can influence the conductance by examining the transmission in the energy window E F ± 0.1 eV. In the case of the spin up transmission, the 2D BZ distribution and the amplitude of the transmission for all structures remains practically the same. However, some changes were observed in the spin down transmission. The likely cause is the presence of minority interfacial resonant states (IRS) [10]. This discrepancy should not influence the generality of the results presented since the contributions from the spin down channel to the G P are much smaller than those of the spin up channel.

C. Effect of shifting the vacancy on the transmission

We now examine the impact of the position of the F and M centers within the barrier on spin-polarized transmission. Although we did not observe any significant change in the layer-projected DOS upon moving the vacancy to the interface, the transmission was nevertheless profoundly altered. Fig. 7(a) and 8(a) respectively show the 2D BZ spindependent transmissions with F and M centers generated in the interfacial MgO layer. Here, the M center always remains in the plane parallel to the interfaces. Interestingly, we found that the transmission distribution, for both defect types, is almost the same as for the ideal junction with the peaks amplitude very close to the ideal case (compare with Fig. 6). The calculated TMR reaches about 4261% and 3911% for the F and M center, respectively, i.e.

are of the same order of magnitude as for the ideal junction (see Tab. III).

When vacancies are on the second layer from the interface, the transmission decreases and we observe additional sharp spikes in the P↑ channel (Figs. To understand the changes in the transmission when varying the F/M centers position, we analyzed a real space distribution of scattering states at the Γ point. We discuss here the MTJ's P magnetic state, focusing on the spin up channel since its transmission strongly drives the ensuing spintronic performance. In the spin up channel we focus on the ∆ 1 symmetry since is has the smallest attenuation rate within MgO barrier and the biggest impact on the resulting transmission. Fig. 9 shows the density of a ∆ 1 scattering state, summed over the xy plane, as a function of the position z along the transport direction for various defect configurations. Clearly, in the presence of vacancies, the amplitude of the ∆ 1 channel is decreased with respect to the ideal case. Interestingly, M centers systematically yield a higher transmission amplitude than F centers at all defect positions within the barrier. MgO layer. In that case, the amplitude of the scattering stated is just slightly lowered with respect to the ideal situation, and the ∆ 1 channel is still transmitted from the left to the right electrode. Moreover, these graphs indicate that the further from the interface a vacancy is, the bigger the difference between F and M centers. The difference in distribution of the ∆ 1 channel for F and M center in the third layer of MgO (Fig. 10) can explain the resulting values of TMR, 145% and 304% respectively, which underscores the synergistic spintronic role of M centers compared to F centers. All the data are normalized and the same logarithmic scale is used.

The overall picture is that the transmission of the spintronically crucial ∆ 1 spin up channel in the MTJ's P magnetic state in the presence of vacancies 1) is higher in presence of an M center rather than an F center for a defect positioned at the center of the barrier, and 2) is close to that of an ideal junction when either vacancy type is positioned near the interface. This second point sheds precious light into how a MgO-class MTJ can experimentally exhibit both high TMR and a low barrier height. Indeed, the MgO barrier is often formed atop the FeCoB metallic surface by sputtering metallic Mg, followed by an oxidation step [73,74].

Avoiding the oxidation of the lower FeCoB interface can naturally lead to the presence of oxygen vacancies within the first ML of MgO. Interfacial oxygen vacancies also play a role in promoting perpendicular magnetic anisotropy for ultrathin ferromagnetic films in MTJs [75,76]. Our results show that, counterintuitively, such states can maintain nearideal levels of TMR and promote the low barrier height needed [START_REF] Halisdemir | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF] for spin transfer torque.

Finally, our theoretical calculations indicate that a the formation energy of a M(F) center created at the interfacial MgO layer is lower by 2 eV (1 eV) than that of one one located on the center ML of a 5ML-thick MgO layer. This indicates that, indeed, interfacial defects are most probable to occur.

D. Rotating the M center

We now examine the impact on transmission of rotating the M center plane so that it is shared between two adjacent MgO layers. Prior complex band structure calculations [47] indicate that if the M center is located on two neighboring MgO layers parallel to the interface, the attenuation coefficient for the ∆ 1 channel can be slightly smaller or comparable to the ideal case. To verify whether this attenuation is reflected in the transmission, we considered a symmetric junction with 6 ML of MgO, such that the M center is shared between the two middle layers. We also increased the number of electrode layers included in the scattering region to ensure proper geometrical matching at the interfaces.

Fig. 12 presents the transmission in the MTJ's P magnetic state for both spin channels, and in the AP magnetic state. Comparing with the ideal case (Fig. 6), we find that both the P↓ and AP transmissions are practically unaffected by the defect. The spin up transmission is even more concentrated around the Γ point than before. The TMR value reaches 1423%, which is as high as when the F/M center is positioned next to the interfacial layer.

If we compare the spin up transmission in Fig. 12 and 8(c), we infer that the shape of the transmission reflects the symmetry/orientation of the M center. Indeed, when the M center is generated in a MgO plane parallel to the interfaces, the two oxygen vacancies lie along the diagonal, and a propagating electron simultaneously encounters both oxygen vacancies. This explains the elongation of the transmission peak along the diagonal of the plane in Fig. 8(c). On the other hand, when the M center is partially along the transport direction, i.e., in the yz plane, such that the two oxygen vacancies are in adjacent xy planes, the propagating electron reaches the first oxygen vacancy and then the second. As a result, the transmission is now along the k y direction in the BZ (Fig. 12). The transmission peaks are of same intensity because, owing to the symmetrical MgO spacer, electrons propagating from the left and the right electrodes see the same potential landscape.

This preservation of high TMR thanks to a P↑ transmission channel that is concentrated Incorporating either a F or M center within a Fe/MgO/Fe magnetic tunnel junction can decrease the transmission of the P↑ channel, while increasing somewhat that of the P↓ and AP channels. As a result the theoretical TMR can drop by up to two orders of magnitude, from 10000% to 100%. Overall, M centers tend to maintain a transmission maximum at the Γ point for the P↑ channel, with only a small broadening, while F centers introduce scattering to higher k-vectors, thereby decreasing the channel's conductance. F and M centers induce only small increase in the P↓ and AP conductances. Consequently, the TMR is generally higher for transport across M centers than F centers. Since the calculated formation energy of an M center of 14.53 eV is lower than that of two F centers (14.65 eV) by 0.12 eV [67], annealing can induce the preferential presence of M centers over F centers, which in turn promotes higher spintronic performance [47]. This annealing-induced fusion of two F centers into a M center is possible only for a large nominal vacancy concentration as found in MgO grown by dc-sputtering Mg and post-deposition oxidation [47], rather than rf-sputtering MgO [69], in line with the high diffusion barrier for an oxygen vacancy in ideal MgO [68].

Our study indicates that the position of F and M centers crucially impacts magnetotransport. Compared to the case of an ideal junction, defects located on the interfacial MgO layer induce practically no change in either the shape of the transmission distribution or its amplitude. The resulting TMR reaches around 4000%, and the system amounts to an ideal MTJ with a barrier of reduced height and effective thickness. Moving the defect away from the interface reduces the P↑ transmission, and thus TMR, especially for the F center. This theoretical insight is compatible with the likely presence in experiments of oxygen vacancies at the lower MTJ interface when the MgO barrier is grown by oxidizing thin layers of metallic Mg deposited atop the lower ferromagnetic metallic electrode while avoiding the latter's oxidation [73,74]. It is also in line with the role of interfacial oxygen vacancies in promoting perpendicular magnetic anisotropy in the adjacent ultrathin ferromagnetic films [75,76] of MgO-class MTJs with perpendicular magnetization. Finally, we find that orienting a M center at the barrier center so as to partly point along the transmission direction yields TMR 1000%.

Our study thus identifies conditions on the nature and positioning within the MgO of single/double oxygen vacancies so as to obtain predicted TMR values in excess of 1000% in MTJs with low barrier heights, in line with TMR amplitudes reported experimentally [13].

Our theoretical results thus reconcile the simultaneous presence of high TMR and low barrier heights in MgO-class MTJs by ascribing them to the presence of oxygen vacancies. Looking ahead, the respectively 3 ML and 3-4 ML effective physical size of the F and M centers condition not only the minimum barrier thickness for sizeable TMR (around 3 ML [77]), but also the MTJ's lateral size. Technological progress has enabled the demonstration of working MTJs with a lateral size down to 4.3 nm [78]. Experiments are thus approaching the 7-8 ML (i.e ≈2 nm) limit estimated for a M center to retain its electronic properties [47]. Our work provides a much-needed theoretical basis to move beyond the mostly unsuspected, fortuitous defect engineering of spintronic performance that has thus far propelled MgObased spintronics and its applications.

FIG. 1 .

 1 FIG. 1. Schematic representation of Fe/MgO junction with lattice parameters as indicated (Red atoms are Fe, blue O and grey Mg).

FIG. 2 .

 2 FIG. 2. Band structure and orbital projected DOS for M-MgO showing the symmetries of the bands. Insets show the electron distribution (linear scale) for each of the ground state energy levels of the M center. Red color indicates the region with the highest electron density 0.1 electron/Å 3 while blue denotes zero density.

Figure 3 FIG. 3 . 2 FIG. 4 .

 3324 Figure3presents the spatial distribution of the electron density for both the ground (panels a/c) and excited (panels b/d) states of M-MgO. In panels a/b (c/d), a 214-atom (62-atom) supercell was used. We observe how neighboring oxygen vacancies hybridize to create an M center. As expected from the band structure plots, the electrons remaining after oxygen removal are localized on the vacancy sites and the electrons are distributed among the vacancies. Since the M center's excited state lies within the conduction band states, a nonzero electron density is present on atoms far from the defect. The electron density plots also reveal a hybridization between the M center states and the nearest oxygen ions, thereby showing that the defect level's excited state is indeed mostly of oxygen p character.The spatial electron density of the M center is fully isolated from the periodic images in the 214-atom supercell. Thus, the lateral extent of the M center spans 1 ML on either side of the oxygen vacancy sites that define the M center. On the other hand, spatial overlap develops between the electron density of M centers in the 62-atom supercell calculation. This means that M centers separated by 2 ML of MgO will experience electronic interactions.The electron distribution also indicates that the ground states of the M center are mostly s-like, while the excited states are of p-like character.

FIG. 5 .

 5 FIG. 5. Spin-polarized layer-projected DOS for Fe/F-MgO/Fe, Fe/M-MgO/Fe, FeCo/F-MgO/FeCo and FeCo/M-MgO/FeCo system with the F and M centers placed on the 4th layer.

FIG. 6 .

 6 FIG. 6. Transmission in the two dimensional Brillouin zone (2D BZ) across an ideal Fe/MgO(5ML)/Fe junction in its P magnetic state for (a) the spin up channel, (b) the spin down channel, and (c) in its AP magnetic state.

  7(b) and 8(b)). The P↓ and AP transmission distributions are only slightly affected. Note that the layer alternation also causes a rotation of the M center within the xy plane when we go from one layer to the next one. This explains the observed rotation in the transmission amplitude in the 2D BZ (compare for example panel (c) and (d) in Fig.8).

4 FIG. 7 .

 47 FIG. 7. Parallel alignment spin up (left) spin down (middle) and antiparallel (right) transmissions in the 2D BZ for Fe/F-MgO/Fe junction for 5 ML of MgO with F center in a) first ML (TMR=4261%), b) second ML (TMR=1239%) and c) third ML (TMR=145%).

4 FIG. 8 .

 48 FIG. 8. Parallel alignment spin up (left) spin down (middle) and antiparallel (right) transmissions in the 2D BZ for Fe/M-MgO/Fe junction for 5 ML of MgO with M center in a) first ML (TMR=3911%), b) second ML (TMR=1135%) and c) third ML (TMR=315%).

Fig. 10

 10 Fig.10shows the ∆ 1 scattering states at the Fermi level across a Fe/MgO(5ML)/Fe junction for the ideal case and the various positions of the F and M centers. All the data are normalized and the same logarithmic scale is used for comparison purposes.

FIG. 9 .

 9 FIG. 9. Spin up ∆ 1 scattering state distribution along the transport direction z for various vacancy configurations of Fe/F-MgO(5ML)/Fe and Fe/M-MgO(5ML)/Fe in the MTJ's P magnetic state.

  FIG. 10. 2D representation of the spin up ∆ 1 scattering states in the x, z plane for various vacancy configurations of . All the data are normalized and the same logarithmic scale is used (see Fig.9.

FIG. 11 .

 11 FIG. 11.Transmission in the two dimensional Brillouin zone (2D BZ) across an ideal Fe/MgO(5ML)/Fe junction in its P magnetic state for (a) the spin up channel, (b) the spin down channel, and (c) in its AP magnetic state.

TABLE I .

 I Experimental TMR and barrier heights for MTJs based on MgO.

TABLE II .

 II Total spin polarized transmissions ×10 4 and TMR for Fe/MgO/Fe, Fe/F-MgO/Fe and Fe/M-MgO/Fe junctions, each with 5 and 7 ML of MgO. The F/M center is located in the middle layer.

	P-UP P-DOWN AP TMR [%]
	Ideal 5ML 79.0	0.46	1.0	7850
	F (5ML) 7.21	0.63	3.2	145
	M (5ML) 17.1	1.47	4.5	315
	Ideal 7ML 5.3	0.003 0.03 15770
	F (7ML) 0.12	0.006 0.03	304
	M (7ML) 0.62	0.007 0.04 1624

TABLE III .

 III Total spin polarized transmission×104 and TMR for Fe/F-MgO(5ML)/Fe and Fe/M-MgO(5ML)/Fe junctions with vacancies shifted within MgO layers.

	P-UP P-DOWN AP TMR [%]
	F in 1st ML 54.5	0.4	1.26 4261
	F in 2nd ML 20.3	0.4	1.55 1239
	F in 3nd ML 7.21	0.63	3.2	145
	M in 1st ML 67.6	0.6	1.70 3911
	M in 2nd ML 30.3	0.5	2.49 1135
	M in 3nd ML 17.1	1.47	4.5	315

ACKNOWLEDGMENTS

This work was performed using HPC resources from the Strasbourg Mesocenter and from the GENCI-CINES Grant gem1100. MB acknowledges funding from the Agence Nationale de la Recherche (ANR-14-CE26-0009-01) and the Labex NIE "Symmix" project (ANR-11-LABX-0058 NIE) and thanks F. Schleicher, D. Lacour, M. Hehn, F. Montaigne, S. Boukari and W. Weber for useful discussions.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

[1] I. Žutić, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004).


J M Aubry
email: jmaubry@math.cnrs.fr

C Butucea

K Meziani
email: meziani@math.jussieu.fr

State estimation in quantum homodyne tomography with noisy data

Keywords: density matrix, Gaussian noise, L 2 -risk, nonparametric estimation, pattern functions, quantum homodyne tomography, quantum state, Radon transform, Wigner function AMS 2000 subject classifications : 62G05, 62G20, 81V80

In the framework of noisy quantum homodyne tomography with efficiency parameter 0 < η ≤ 1, we propose two estimators of a quantum state whose density matrix elements ρ m,n decrease like e -B(m+n) r/2 , for fixed known B > 0 and 0 < r ≤ 2. The first procedure estimates the matrix coefficients by a projection method on the pattern functions (that we introduce here for 0 < η ≤ 1/2), the second procedure is a kernel estimator of the associated Wigner function. We compute the convergence rates of these estimators, in L 2 risk.

Introduction

Experiments in quantum optics consist in creating, manipulating and measuring quantum states of light. The technique called quantum homodyne tomography allows to retrieve partial, noisy information from which the state is to be recovered : this is the subject of the present chapter.

Quantum states

Mathematically, the main concepts of quantum mechanics are formulated in the language of selfadjoint operators acting on Hilbert spaces. To every quantum system one can associate a complex Hilbert space H whose vectors represent the wave functions of the system. These vectors are identified to projection operators, or pure states. In general, a state is a mixture of pure states described by a compact operator ρ on H having the following properties :

1. Selfadjoint : ρ = ρ * , where ρ * is the adjoint of ρ.

2. Positive : ρ ≥ 0, or equivalently ψ, ρψ ≥ 0 for all ψ ∈ H.

Trace one : tr

(ρ) = 1.
When H is separable, endowed with a countable orthonormal basis, the operator ρ is identified to a density matrix [ρ m,n ] m,n∈N .

The positivity property implies that all the eigenvalues of ρ are nonegative and by the trace property, they sum up to one. In the case of the finite dimensional Hilbert space C d , the density matrix is simply a positive semi-definite d × d matrix of trace one. Our setup from now on will be H = L 2 (R), in which case we employ the orthonormal Fock basis made of the Hermite functions

h m (x) := (2 m m! √ π) -1 2 H m (x)e -x 2 2 (1)
where H m (x) := (-1) m e x 2 d m dx m e -x 2 is the m-th Hermite polynomial. Generalizations to higher dimensions are straightforward.

To each state ρ corresponds a Wigner distribution W ρ , which is defined via its Fourier transform in the way indicated by equation (2) : W ρ (u, v) := e -i(uq+vp) W ρ (q, p)dqdp := Tr ρ exp(-iuQ -ivP)

where Q and P are canonically conjugate observables (e.g. electric and magnetic fields) satisfying the commutation relation [Q, P] = i (we assume a choice of units such that = 1). It is easily checked that W ρ is real-valued, has integral R 2 W ρ (q, p)dqdp = 1 and uniform bound |W ρ (q, p)| ≤ 1 π . For any φ ∈ R, the Wigner distribution allows one to easily recover the probability density x → p ρ (x, φ) of Q cos φ + P sin φ by

p ρ (x, φ) = R[W ρ ](x, φ), ( 3
)
where R is the Radon transform defined in equation ( 4)

R[W ρ ](x, φ) = ∞ -∞
W ρ (x cos φ -t sin φ, x sin φ + t cos φ)dt.

Moreover, the correspondence between ρ and W ρ is one to one and isometric with respect to the L 2 norms as in equation ( 5) :

W ρ 2 2 := |W ρ (q, p)| 2 dqdp = 1 2π ρ 2 2 := 1 2π ∞ j,k=0 |ρ jk | 2 . (5)
From now on we denote by •, • and • the usual Euclidian scalar product and norm, while C(•) will denote positive constants depending on parameters given in the parentheses.

We suppose that the unknown state belongs to the class R(B, r) for B > 0 and 0 < r ≤ 2 defined by R(B, r) := {ρ quantum state : |ρ m,n | ≤ exp(-B(m + n) r/2 )}. [START_REF] Cavalier | Efficient estimation of a density in a problem of tomography[END_REF] For simplicity, we have chosen to express the results relative to a class which is the intersection of the (positive) ball of radius 1 in some Banach space with the hyperplane tr(ρ) = 1.

Another radius for the class would only change the constant C in front of the asymptotic rates of convergence that we will find.

As it will be made precise in Propositions 1 and 2, quantum states in the class given in [START_REF] Cavalier | Efficient estimation of a density in a problem of tomography[END_REF] have fast decreasing and very smooth Wigner functions. From the physical point of view, the choice of such a class of Wigner functions seems to be quite reasonable considering that typical states ρ prepared in the laboratory do satisfy this type of condition.

Statistical model

Let us describe the statistical model. Consider (X 1 , Φ 1 ), . . . , (X n , Φ n ) independent identically distributed random variables with values in R × [0, π] and distribution P ρ having density p ρ (x, φ) (given by (3) with respect to 1 π λ, λ being the Lebesgue measure on R×[0, π]. The aim is to recover the density matrix ρ and the Wigner function W ρ from the observations.

However, there is a slight complication. What we observe are not the variables (X ℓ , Φ ℓ ) but the noisy ones (Y ℓ , Φ ℓ ), where

Y ℓ := √ ηX ℓ + (1 -η)/2 ξ ℓ , (7)
with ξ ℓ a sequence of independent identically distributed standard Gaussians which are independent of all (X j , Φ j ). The detection efficiency parameter 0 < η ≤ 1 is known from the calibration of the apparatus and we denote by N η the centered Gaussian density of variance (1 -η)/2, and N η its Fourier transform. Then the density p η ρ of (Y ℓ , Φ ℓ ) is given by the convolution of the density

p ρ (•/ √ η, φ)/ √ η with N η p η ρ (y, φ) = ∞ -∞ 1 √ η p ρ y -x √ η , φ N η (x)dx =: 1 √ η p ρ • √ η , φ * N η (y).
In the Fourier domain this relation becomes

F 1 [p η ρ (•, φ)](t) = F 1 [p ρ (•, φ)](t √ η) N η (t), (8)
where F 1 denotes the Fourier transform with respect to the first variable.

The theoretical foundation of quantum homodyne tomography was outlined in [START_REF] Vogel | Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase[END_REF] and has inspired the first experiments determining the quantum state of a light field, initially with optical pulses in [START_REF] Smithey | Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomogra-26 phy : Application to squeezed states and the vacuum[END_REF][START_REF] Smithey | Experimental determination of number-phase uncertainty relations[END_REF][START_REF] Leonhardt | Measuring the Quantum State of Light[END_REF]. The reconstruction of the density from averages of data has been discussed or studied in [START_REF] D'ariano | Detection of the density matrix through optical homodyne tomography without filtered back projection[END_REF][START_REF] D'ariano | Homodyne detection of the density matrix of the radiation field[END_REF][START_REF] Leonhardt | Tomographic reconstruction of the density matrix via pattern functions[END_REF][START_REF] Artiles | An invitation to quantum tomography[END_REF] for η = 1 (no photon loss). Max-likelihood methods have been studied in [START_REF] Banaszek | Maximum-likelihood estimation of the density matrix[END_REF][START_REF] Artiles | An invitation to quantum tomography[END_REF][START_REF] D'ariano | Homodyne tomography and the reconstruction of quantum states of light[END_REF][START_REF] Gut ¸ȃ | Maximum likelihood estimation of the density matrix through quantum tomography[END_REF] and procedure using adaptive tomographic kernels to minimize the variance has been proposed in [START_REF] Paris | Adaptive quantum tomography[END_REF]. The estimation of the density matrix of a quantum state of light in case of efficiency parameter 1 2 < η ≤ 1 has been discussed in [START_REF] D'ariano | Tomographic measurement of the density matrix of the radiation field[END_REF][START_REF] D'ariano | Homodyne tomography and the reconstruction of quantum states of light[END_REF][START_REF] D'ariano | Tomographic methods for universal estimation in quantum optics[END_REF] and considered in [START_REF] Richter | Pattern functions used in tomographic reconstruction of photon statistics revisited[END_REF] via the pattern functions for the diagonal elements.

Outline of the results

The goal of this chapter is to define estimators of both the density matrix and the Wigner function and to compare their performance in L 2 risk. In order to compute estimation risks and to tune the underlying parameters, we define a realistic class of quantum states R(B, r), depending on parameters B > 0 and 0 < r ≤ 2, in which the elements of the density matrix decrease rapidly.

In Section 2, we prove that the fast decay of the elements of the density matrix implies both rapid decay of the Wigner function and of its Fourier transform, allowing us to translate the classes R(B, r) in terms of Wigner functions.

In Section 3, we give estimators of the density matrix ρ. The legend was somehow forged that no estimation of the matrix is possible when 0 < η ≤ 1/2. The physicists argue that their machines actually have high detection efficiency, around 0.8 ; it is nevertheless satisfying to be able to solve this problem in any noise condition. We give here the socalled pattern functions to use for estimating the density matrix in the noisy case with any value of η between 0 and 1. These pattern functions allow us to solve an inverse problem which becomes (severly) ill-posed when 0 < η ≤ 1/2. In this case, we regularize the inverse problem and this introduces a smoothing parameter which we will choose in an optimal way. We compute the upper bounds for the rates achieved by our methods, with L 2 risk measure.

In Section 4, we study a kernel estimator of the Wigner function in L 2 risk, over the same class of Wigner functions. It is a truncated version of the estimator in [START_REF] Butucea | Minimax and adaptive estimation of the Wigner function in quantum homodyne tomography with noisy data[END_REF] and tuned accordingly. We compute upper bounds for the rates of convergence of this estimator in L 2 risk.

To conclude, we may infer that the performances of both estimators are comparable. We obtain nearly polynomial rates for the case r = 2 and intermediate rates for 0 < r < 2 (faster than any logarithm, but slower than any polynomial). It is convenient to have methods to estimate directly both representations of a quantum state. The estimator of the matrix ρ can be more easily projected on the space of proper quantum states. On the other hand, we may capture some features of the quantum states more easily on the Wigner function, for instance when this function has significant negative parts, the fact that the quantum state is non classical.

Decrease and smoothness of the Wigner distribution

We recall that the Wigner distribution W ρ was defined in the introduction. In the Fock basis, we can write W ρ in terms of the density matrix [ρ m,n ] as follows (see Leonhardt [START_REF] Leonhardt | Measuring the Quantum State of Light[END_REF] for the details).

W ρ (q, p) = m,n ρ m,n W m,n (q, p) where W m,n (q, p) = 1 π e 2ipx h m (q -x)h n (q + x)dx. (9)
It can be seen that W m,n (q, p) = W n,m (q, -p) and if m ≥ n,

W m,n (q, p) = (-1) m π n! m! 1 2 e -(q 2 +p 2 ) × √ 2(ip -q) m-n L m-n n 2q 2 + 2p 2 (10)
thus, writing

z := q 2 + p 2 , l m,n (z) := |W m,n (q, p)| = 2 m-n 2 π n! m! 1 2 e -z 2 z m-n L m-n n (2z 2 ) (11)
where L α n (x) := (n!) -1 e x x -α d n dx n (e -x x n+α ) is the Laguerre polynomial of degree n and order α. Concerning the Fourier transforms, we also recall that

W m,n (q, p) = (-i) m+n 2 W m,n q 2 , p 2 . ( 12
)
In this section we show how a decrease condition on the coefficients of the density matrix translates on the corresponding Wigner distribution. First the case r < 2 :

Proposition 1. Assume that 0 < r < 2 and that there exists B > 0 such that, for all

m ≥ n, |ρ m,n | ≤ e -B(m+n) r/2 .
Then for all β < B, there exists z 0 (depending explicitly on r, B, β, see proof ) such that

z := q 2 + p 2 ≥ z 0 implies |W ρ (q, p)| ≤ A(z)e -βz r (13)
as well as

W ρ (q, p) ≤ A(z/2)e -β(z/2) r (14)
where A(z) := 1 π m,n e -B(m+n) r/2 + 4 Br z 4-r .

If r = 2, the result is a little different :

Proposition 2. Suppose that there exists B > 0 such that, for all m ≥ n, m+n) .

|ρ m,n | ≤ e -B(
Then there exists z 0 such that z := q 2 + p 2 ≥ z 0 implies

|W ρ (q, p)| ≤ A(z)e - B (1+ √ B) 2 z 2 (15)
as well as

W ρ (q, p) ≤ A(z/2)e - B (1+ √ B) 2 (z/2) 2 (16)
for

A(z) = 1 π m,n e -B(m+n) + 2e B B(1+ √ B) 2 z 2 .
Note that

B (1+ √ B)
2 < min(B, 1). Even when B is very large, we cannot hope to obtain a faster decrease because e -z 2 is the decrease rate of the basis functions themselves (Lemma 2).

The proof of these propositions is defered to Appendix 5. More general results and converses are studied in [START_REF] Aubry | Ultrarapidly decreasing ultradifferentiable functions, Wigner distributions and density matrices[END_REF]. Let us now state a few general utility lemmata.

Lemma 1. Let y and w be two C 2 functions : [x 0 , +∞) → (0, +∞) such that y ′ (x) → 0, w is bounded, satisfying the differential equations

y ′′ (x) = φ(x)y(x) w ′′ (x) = ψ(x)w(x),
with continuous φ(x) ≤ ψ(x), and initial conditions y(x 0 ) = w(x 0 ). Then for all x ≥ x 0 , w(x) ≤ y(x).

Démonstration. Suppose that there exists x 1 ≥ x 0 where w(x 1 ) > y(x 1 ). Then for some

x 2 ∈ [x 0 , x 1 ] we have w ′ (x 2 ) > y ′ (x 2 ) and w(x 2 ) ≥ y(x 2 ). Consequently, for all x ≥ x 2 , w ′′ (x) -y ′′ (x) ≥ 0, and w ′ (x) -y ′ (x) ≥ w ′ (x 2 ) -y ′ (x 2 ). When x → ∞, lim inf w ′ (x) ≥ w ′ (x 2 ) -y ′ (x 2 ) > 0, which contradicts the boundedness of w.
This lemma is used to prove a bound on the Laguerre functions.

Lemma 2. For all m, n ∈ N and s :=

√ m + n + 1, for all z ≥ 0, l m,n (z) ≤ 1 π    1 if 0 ≤ z ≤ s e -(z-s) 2 if z ≥ s. (17)
Démonstration. When z ≤ s, the result follows from the uniform bound on Wigner functions obtained by applying the Cauchy-Schwarz inequality to [START_REF] D'ariano | Homodyne detection of the density matrix of the radiation field[END_REF]. When z ≥ s, L α n (2z 2 ) doesn't vanish and keeps the same sign as L α n (2s 2 ). Now, as it can be seen from [27, 5.1.2], the function w(z) := √ zl m,n (z) satisfies the differential equation

w ′′ = (4(z 2 -s 2 ) + α 2 -1/4 z 2 )z. On the other hand, y(z) := √ sl m,n (s)e -(z-s) 2 satisfies y ′′ = (4(z -s) 2 -2)y. When z ≥ s, 4(z -s) 2 -2 < 4(z 2 -s 2 ) + α 2 -1/4 z 2 (18)
from which we conclude with Lemma 1 that w(z) ≤ y(z).

Finally, a lemma to bound the tail of a series.

Lemma 3. If ν > 0 and C > 0, there exists a z 0 such that z ≥ z 0 implies m+n≥z e -C(m+n) ν ≤ 2 Cν z 2-ν e -Cz ν . ( 19
)
Démonstration. First notice that m+n≥z e -C(m+n) ν = t≥z (t + 1)e -Ct ν ≤ ∞ z (t + 1)e -Ct ν dt.
When t ≥ z and z is large enough, we have

∞ z (t + 1)e -Ct ν dt ≤ 2 Cν ∞ z Cνt -(2 -ν)t 1-ν e -Ct ν dt ≤ 2 Cν z 2-ν e -Cz ν
which is what we needed to prove.

Density matrix estimation

The aim of this part is to estimate the density matrix ρ in the Fock basis directely from the data (Y i , Φ i ) i=1,...,n . We show that for 0 < η ≤ 1/2 it is still possible to estimate the density matrix with an error of estimation tending to 0 as n tends to infinity (Theorem 3). In both cases (η > 1 2 and η ≤ 1 2 ), we construct an estimator of the density matrix (ρ j,k ) j,k≤N -1 from a sample of QHT data. We give theoretical results for our estimator when the quantum state ρ is in the class of density matrix with decreasing elements defined in (6).

Pattern functions

The matrix elements ρ j,k of the state ρ in the Fock basis (1) can be expressed as kernel integrals : for all j, k ∈ N,

ρ j,k = 1 π π 0 p ρ (x, φ)f j,k (x)e -i(k-j)φ dφdx (20)
where f j,k = f k,j are bounded real functions called pattern functions in quantum homodyne literature. A concrete expression for their Fourier transform using Laguerre polynomials was found in [START_REF] Richter | Realistic pattern functions for optical homodyne tomography and determination of specific expectation values[END_REF] :

for j ≥ k, fk,j (t) = 2π 2 |t| W j,k (t, 0) = π(-i) j-k 2 k-j k! j! |t|t j-k e -t 2 4 L j-k k ( t 2 2 ). (21)
where fk,j denotes the Fourier transform of the Pattern function f k,j .

Let us state the lemmata which are used to prove upper bounds in Propositions 3, 4 and 5.

Lemma 4. There exist constants C 2 , C ∞ such that

N j+k=0 f k,j 2 2 ≤ C 2 N
17 6 and

N j+k=0 f k,j 2 ∞ ≤ C ∞ N 10 3 .
This is a slight improvement over [1, Lemma 1].

Démonstration. By symmetry we can restrict the sum to j ≥ k. For fixed k and j we have fk,j In view of ( 21), the main result in [START_REF] Krasikov | Inequalities for orthonormal Laguerre polynomials[END_REF] can be rewritten as follows : if k ≥ 35 and j -k ≥ 24, then fk,j

2 ∞ ≤ 2888π 2 (j + 1) 1 2 k -1 6 . (22)
In consequence, for these values of k and j, fk,j

2 2 ≤ C(jk -1 6 + j 1 2 k 1 3 ). ( 23
)
On the other hand, a classical bound on Laguerre polynomials found in [START_REF] Szegö | Orthogonal polynomials[END_REF] yields that, for fixed values of j -k, fk,j

2 ∞ ≤ Ck 1 3
, hence for all k ≥ 35 and j -k < 24, fk,j

2 2 ≤ C(j 1 2 k 1 3 + k 5 6 ). ( 24
)
When k < 35, we can use another result in [START_REF] Krasikov | Inequalities for Laguerre polynomials[END_REF] which gives fk,j

2 ∞ ≤ Ck 1 6 j 1 2 independently of j -k, thus fk,j 2 2 ≤ Cj. (25)
Comparing ( 23), ( 24) and ( 25) we see that when N is large enough, in the sum over 0 ≤ j, k ≤ N , the terms k ≥ 35, j -k ≥ 24 dominate and ( 23) yields the first inequality.

The second inequality is obtained by doing a similar computation, starting with f j,k ∞ ≤ fj,k 1 and using [START_REF] Meziani | Nonparametric estimation of the purity of a quantum state in quantum homodyne tomography with noisy data[END_REF] to bound fj,k

2 1 ≤ C(j 3 2 k -1 6 + j 1 2 k 5 6 )
when k ≥ 35 and j -k ≥ 24.

In the presence of noise, it is necessary to adapt the pattern functions as follows. From now on, we shall use the notation γ := 1-η 4η . When 1 2 < η ≤ 1, we denote by f η k,j the function which has the following Fourier transform :

f η k,j (t) := fk,j (t)e γt 2 . (26)
When 0 < η ≤ 1 2 , we introduce a cut-off parameter δ > 0 and define f η,δ k,j via its Fourier transform :

f η,δ k,j (t) := fk,j (t)e γt 2 I |t| ≤ 1 δ . ( 27
)
Then we compute bounds on these pattern functions.

Lemma 5. For 1 > η > 1/2, there exist constants C η 2 and C η ∞ such that N j+k=0 f η k,j 2 2 ≤ C η 2 N 5 6 e 8γN and N j+k=0 f η k,j 2 ∞ ≤ C η ∞ N 1 3 e 8γN .
Démonstration. The proof is similar to the previous one and we skip some details. Once again we assume j ≥ k and write

f η k,j 2
e 2γt 2 dt ≤ C fk,j 2 ∞ s -1 e 8γs 2 .
In the sum we are considering the terms k ≥ 35 and j -k ≥ 24 are dominant and, once again thanks to ( 22), remembering that s =

√ j + k + 1, f η k,j 2 2 ≤ Ck -1 6 e 8γ(j+k)
hence the first inequality.

The second inequality is, in the same fashion, based on

f η k,j 2 ∞ ≤ f η k,j 2 1 ≤ C j 1 4 k -1 12 |t|<2s e γt 2 dt 2 ≤ Cj -1 2 k -1 6 e 8γ(j+k)
when k ≥ 35 and j -k ≥ 24, and the bound on the sum readily follows.

Estimation procedure

For N := N (n) → ∞ and δ := δ(n) → 0, let us define our estimator of ρ j,k for 0

≤ j + k ≤ N -1 by ρη j,k := 1 n n ℓ=1 G j,k Y ℓ √ η , Φ ℓ , (28)
where

G j,k (x, φ) :=    f η j,k (x)e -i(j-k)φ if 1 2 < η ≤ 1 f η,δ j,k (x)e -i(j-k)φ if 0 < η ≤ 1 2
. using the pattern functions defined in [START_REF] Smithey | Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomogra-26 phy : Application to squeezed states and the vacuum[END_REF] and [START_REF] Szegö | Orthogonal polynomials[END_REF]. We assume that the density matrix ρ belongs to the class R(B, r) defined in [START_REF] Cavalier | Efficient estimation of a density in a problem of tomography[END_REF]. In order to evaluate the performance of our estimators we take the L 2 distance on the space of density matrices τ -ρ 2 2 := tr(|τ -

ρ| 2 ) = ∞ j,k=0 |τ j,k -ρ j,k | 2 .
We consider the mean integrated square error (MISE) and split it into a troncature bias term b 2 1 (n), a regularization bias terms b 2 2 (n) and a variance term

σ 2 (n). E   ∞ j,k=0 ρη j,k -ρ j,k 2   = j+k≥N |ρ j,k | 2 + N -1 j+k=0 E[ρ η j,k ] -ρ j,k 2 + N -1 j+k=0 E ρη j,k -E[ρ η j,k ] 2 =: b 2 1 (n) + b 2 2 (n) + σ 2 (n).
The following propositions give upper bounds for b 2 1 (n), b 2 2 (n) and σ 2 (n) in the different cases η = 1, 1/2 < η < 1 or 0 < η ≤ 1/2 and r = 2 or 0 < r < 2. Their proofs are defered to Appendix 5. Proposition 3. Let ρη j,k be the estimator defined by ( 28), for 0 < η < 1, with δ → 0 and N → ∞ as n → ∞, then for all B > 0 and 0

< r ≤ 2, sup ρ∈R(B,r) b 2 1 (n) ≤ c 1 N 2-r/2 e -2BN r/2 (29)
where c 1 is a positive constant depending on B and r.

Proposition 4. Let ρη j,k be the estimator defined by [START_REF] Vardi | A statistical model for positron emission tomography[END_REF], for 0 < η ≤ 1/2, with N → ∞ as n → ∞ and 1/δ ≥ 2 √ N . In the case r = 2, for β := B/(1 + √ B) 2 there exists c 2 , while in the case 0 < r < 2, for any β < B there exists c 2 and n 0 such that for n ≥ n 0 :

sup ρ∈R(B,r) b 2 2 (n) ≤ c 2 N 2 δ 4r-12 e -2β (2δ) r -1 2 ( 1 δ -2 √ N ) 2 . ( 30
)
Note that for 1/2 < η ≤ 1 we have b 2 (n) = 0 for all 0 < r ≤ 2 (ρ η j,k is unbiased). Proposition 5. For ρη j,k the estimator defined by [START_REF] Vardi | A statistical model for positron emission tomography[END_REF],

sup ρ∈R(B,r) σ 2 (n) ≤ c 3 δN 17/6 n e 2γ δ 2 if 0 < η ≤ 1/2 (31) sup ρ∈R(B,r) σ 2 (n) ≤ c ′ 3 N 1/3 n e 8γN if 1/2 < η < 1 (32) sup ρ∈R(B,r) σ 2 (n) ≤ c ′′ 3 N 17 6 n if η = 1 ( 33
)
where c 3 , c ′ 3 are positive constants depending on η.

We measure the accuracy of ρη j,k by the maximal risk over the class R(B, r)

lim sup n→∞ sup ρ∈R(B,r) ϕ -2 n E   ∞ j,k=0 ρη j,k -ρ j,k 2   ≤ C 0 . ( 34
)
where C 0 is a positive constant and ϕ 2 n is a sequence which tends to 0 when n → ∞ and it is the rate of convergence. Cases η = 1 (no noise), 1 2 < η < 1 (weak noise) and 0 < η ≤ 1 2 (strong noise) are studied respectively in Theorems 1, 2 and 3.

Theorem 1. When η = 1, the estimator defined in [START_REF] Vardi | A statistical model for positron emission tomography[END_REF] for the model [START_REF] D'ariano | Tomographic measurement of the density matrix of the radiation field[END_REF], where the unknown state belongs to the class R(B, r), satisfies the upper bound (34) with

ϕ 2 n = log(n) 17 3r n -1 obtained by taking N (n) := log(n) 2B 2 r .
Démonstration. With the proposed N (n) one checks that the bias ( 29) is smaller than the variance (33) which is bounded by a constant times log(n)

17 3r n -1 .
Theorem 2. When 1 2 < η < 1, the estimator defined in [START_REF] Vardi | A statistical model for positron emission tomography[END_REF] for the model [START_REF] D'ariano | Tomographic measurement of the density matrix of the radiation field[END_REF], where the unknown state belongs to the class R(B, r), satisfies the upper bound (34) with -For r = 2,

ϕ 2 n = log(n) 12γ+B 3(4γ+B) n -B 4γ+B with N (n) := log(n) 2(4γ+B) 1 + 2 3 log(log n) log(n)
.

-For 0 < r < 2,

ϕ 2 n = log(n) 2-r/2 e -2BN (n) r/2
where N (n) is the solution of the equation 8γN + 2BN r/2 = log(n).

In that case we have

N (n) = 1 8γ log(n) -2B (8γ) 1+r/2 log(n) r/2 + o(log(n) r/2 ).
Démonstration. When r = 2, the proposed N (n) ensures that the variance (32) is equivalent to the bias [START_REF] Vogel | Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase[END_REF], which is bounded by a constant times log(n)

12γ+B 3(4γ+B) n -B 4γ+B .
When 0 < r < 2, the proposed N (n) makes the variance (32) bounded by a constant times e -2BN (n) r/2 , which is smaller than the bias, the latter being bounded by a constant times

N (n) 2-r/2 e -2BN (n) r/2 .
The asymptotic expansion of N (n) is a standard consequence of its definition by the equation 8γN + 2BN r/2 = log(n).

Theorem 3. When 0 < η ≤ 1 2 , the estimator defined in [START_REF] Vardi | A statistical model for positron emission tomography[END_REF] for the model [START_REF] D'ariano | Tomographic measurement of the density matrix of the radiation field[END_REF], where the unknown state belongs to the class R(B, r), satisfies the upper bound (34) with

ϕ 2 n = N 2-r/2 e -2BN r/2
where N and δ are solutions of the system

   2β (2δ) r + 1 2 ( 1 δ -2 √ N ) 2 + 2γ δ 2 = log(n) 2β (2δ) r + 1 2 ( 1 δ -2 √ N ) 2 -2BN r/2 = (log log(n)) 2 (35)
for arbitrary β < B in the case 0 < r < 2 or    β+4γ 2δ 2 + 1 2 ( 1 δ -2 √ N ) 2 -5 3 log(N ) = log(n) β 2δ 2 + 1 2 ( 1 δ -2 √ N ) 2 -2BN -3 log(N ) = 0 ( 36
)
with β := B (1+ √
B) 2 in the case r = 2. Theses bounds are optimal in the sense that (35) and (36) are obtained by minimizing the sum of the bounds ( 29), ( 30) and (31).

Démonstration. We use the standard notations a

(n) ∼ b(n) if a(n) b(n) → 1 and a(n) ≈ b(n) if there exists a constant M < ∞ such that 1 M ≤ a(n) b(n) ≤ M for all n.
Let us first examine the case 0 < r < 2. Remark that the left-hand term of the second equation in (35) is strictly negative when 1/δ = 2 √ N and increases to ∞ with 1/δ. This proves that the solution satisfies 1/δ > 2 √ N and that Proposition 4 applies. Furthermore, if we suppose that 1/δ √ N is unbounded when n → ∞, then (up to taking a subsequence) by the first equation

1 2 +2γ δ 2
∼ log(n) whereas, by subtracting the two, 2γ δ 2 ∼ log(n), which is contradictory. So 1/δ ≈ √ N and we deduce that N ≈ log(n). Then (30) yields

log b 2 2 (n) N 2-r/2 e -2BN r/2 ≤ (4r -12) log(δ) + r 2 log(N ) -(log log(n)) 2 → -∞ whereas (31) gives log σ 2 (n) N 2-r/2 e -2BN r/2 ≤ log(δ) + ( 5 6 + r 2 ) log(N ) -(log log(n)) 2 → -∞.
We see that the dominant term is the bound (29) on b 2 1 (n), hence the result. When r = 2, the same reasoning as above yields

1/δ > 2 √ N , 1/δ ≈ √ N and N ≈ log(n).
Then the right-hand side of (30) and (31) are of the same order as N e -2BN , which is the bound (29) on b 2 1 (n).

Wigner function estimation 4.1 Kernel estimator

We describe now the direct estimation method for the Wigner function. For the problem of estimating a probability density f : R 2 → R directly from data (X ℓ , Φ ℓ ) with density R[f ] we refer to the literature on X-ray tomography and PET, studied by [START_REF] Vardi | A statistical model for positron emission tomography[END_REF][START_REF] Korostelëv | Minimax theory of image reconstruction[END_REF][START_REF] Iain | Speed of estimation in positron emission tomography and related inverse problems[END_REF][START_REF] Cavalier | Efficient estimation of a density in a problem of tomography[END_REF] and many other references therein. In the context of tomography of bounded objects with noisy observations [START_REF] Goldenshluger | On the shape-from-moments problem and recov-ering edges from noisy Radon data[END_REF] solved the problem of estimating the borders of the object (the support). The estimation of a quadratic functional of the Wigner function has been treated in [START_REF] Meziani | Nonparametric estimation of the purity of a quantum state in quantum homodyne tomography with noisy data[END_REF]. For the problem of Wigner function estimation when no noise is present, we mention the work by [START_REF] Gut | Minimax estimation of the Wigner in quantum homodyne tomography with ideal detectors[END_REF]. They use a kernel estimator and compute sharp minimax results over a class of Wigner functions characterised by their smoothness. In a more recent paper [START_REF] Butucea | Minimax and adaptive estimation of the Wigner function in quantum homodyne tomography with noisy data[END_REF],

Butucea, Gut ¸ȃ and Artiles treated the noisy problem for the pointwise estimation of W ρ ; however the functions needed to prove minimax optimality there do not belong to the class of Wigner functions that we consider here.

In this chapter, as in [START_REF] Butucea | Minimax and adaptive estimation of the Wigner function in quantum homodyne tomography with noisy data[END_REF], we modify the usual tomography kernel in order to take into account the additive noise on the observations and construct a kernel K η h which performs both deconvolution and inverse Radon transform on our data, asymptotically. Let us define the estimator :

W η h (q, p) = 1 πn n ℓ=1 K η h q cos Φ ℓ + p sin Φ ℓ - Y ℓ √ η , ( 37
)
where 0 < η < 1 is a fixed parameter, and the kernel is defined by

K η h (u) = 1 4π 1/h -1/h exp(-iut)|t| N η (t/ √ η) dt, K η h (t) = 1 2 |t| N η (t/ √ η) I(|t| ≤ 1/h), (38)
and h > 0 tends to 0 when n → ∞ in a proper way to be chosen later. For simplicity, let us denote z = (q, p) and [z, φ] = q cos φ + p sin φ, then the estimator can be written :

W η h (z) = 1 πn n ℓ=1 K η h [z, Φ ℓ ] - Y ℓ √ η .
This is a one-step procedure for treating two successive inverse problems. The main difference with the noiseless problem treated by [START_REF] Gut | Minimax estimation of the Wigner in quantum homodyne tomography with ideal detectors[END_REF] is that the deconvolution is more 'difficult' than the inverse Radon transform. In the literature on inverse problems, this problem would be qualified as severely ill-posed, meaning that the noise is dramatically (exponentially) smooth and makes the estimation problem much harder.

L 2 risk estimation

We establish next the rates of estimation of W ρ from i.i.d. observations (Y ℓ , Φ ℓ ), ℓ = 1, . . . , n when the quality of estimation is measured in L 2 distance. In the literature, L 2 tomography is usually performed for boundedly supported functions, see [START_REF] Korostelëv | Minimax theory of image reconstruction[END_REF] and [START_REF] Iain | Speed of estimation in positron emission tomography and related inverse problems[END_REF]. However, most Wigner function do not have a bounded support ! Instead, we use the fact that Wigner functions in the class R(B, r) decrease very fast and show that a properly truncated estimator attains the rates we may expect from the statistical problem of deconvolution in presence of tomography. Thus, we modify the estimator by truncating it over a disc with increasing radius, as n → ∞. Let us denote

D(s n ) = {z = (q, p) ∈ R 2 : z ≤ s n } ,
where s n → ∞ as n → ∞ will be defined in Theorem 4. Let now

W η, * h,n (z) = W η h,n (z)I D(sn) (z). (39)
From now on, we will denote for any function f ,

f 2 D(sn) = D(sn) f 2 (z)dz,
and by D(s n ) the complementary set of D(s n ) in R 2 . Then,

E W η, * h -W ρ 2 2 = E W η h -W ρ 2 D(sn) + W ρ 2 D(sn) = E W η h -E W η h 2 D(sn) + E W η h -W ρ 2 D(sn) + W ρ 2 D(sn) .
When replacing the L 2 norm with the above restricted integral, the upper bound of the bias of the estimator is unchanged, whereas the variance part is infinitely larger than the deconvolution variance in [START_REF] Butucea | Sharp optimality for density deconvolution with dominating bias. i and ii[END_REF]. As the bias is dominating over the variance in this setup, we can still choose a suitable sequence s n so that the same bandwidth is optimal associated to the same optimal rate, provided that W ρ decreases fast enough asymptotically. The following proposition gives upper bounds for the three components of the L 2 risk uniformly over the class R(B, r). Proposition 6. Let (Y ℓ , Φ ℓ ), ℓ = 1, . . . , n be i.i.d. data coming from the model ( 7) and let W η h be an estimator (with h → 0 as n → ∞) of the underlying Wigner function W ρ . We suppose W ρ lies in the class R(B, r), with B > 0 and 0 < r ≤ 2. Then, for s n → ∞ as n → ∞ and n large enough,

sup ρ∈R(B,r) W ρ 2 D(sn) ≤ C 1 s 10-3r n e -2βs r n , sup ρ∈R(B,r) E[ W η h ] -W ρ 2 D(sn) ≤ C 2 h 3r-10 e -2 1-r β h r sup ρ∈R(B,r) E W η h,n -E W η h,n 2
D(sn) ≤ C 3 s 2 n nh exp 2γ h 2 ,
where β < B is defined in Proposition 1 for 0 < r < 2 and

β = B/(1 + √ B) 2 for r = 2, γ = (1 -η)/(4η) > 0, C 1 , C 2 , C 3 are positive constants, C 1 , C 2 , depending on β, B, r and
C 3 depending only on η.
We measure the accuracy of W η, * h by the maximal risk over the class R(B, r)

lim sup n→∞ sup ρ∈R(B,r) E W η, * h -W ρ 2 ϕ -2 n (L 2 ) ≤ C. ( 40
)
where C is a positive constant and ϕ 2 n is a sequence which tends to 0 when n → ∞ and it is the rate of convergence.

In the following Theorem we see the phenomenon which was noticed already : deconvolution with Gaussian type noise is a much harder problem than inverse Radon transform (the tomography part). Theorem 4. Let B > 0, 0 < r ≤ 2 and (Y ℓ , Φ ℓ ), ℓ = 1, . . . , n be i.i.d. data coming from the model [START_REF] D'ariano | Tomographic measurement of the density matrix of the radiation field[END_REF]. Then W η, * h defined in (39) with kernel K η h in (38) satisfies the upper bound (40) with

-For r = 2, put β = B/(1 + √ B) 2 ϕ 2 n = (log n) 16γ+3β 8γ+2β n -β 4γ+β , with s n = (h) -1 and h = 2 4γ+β log n + 1 4γ+β log(log n) -1/2
.

-For 0 < r < 2 and β < B defined in Proposition 1,

ϕ 2 n = h 3r-10 exp - 2 1-r β h r ,
where s n = 1/h and h is the solution of the equation

2 1-r β h r + 2γ h 2 = log n -(log log n) 2 .
Sketch of proof of the upper bounds. By Proposition 6, we get sup

Wρ∈R(B,r) E W η h -W ρ 2 ≤ C 1 s 10-3r n e -2βs r n + C 2 h 3r-10 exp - 2β (2h) r + C 3 s 2 n nh exp 2γ h 2 . =: A 1 + A 2 + A 3
For 0 < r < 2 and by taking derivatives with respect to h and s n , we obtain that the optimal choice verifies the following equations :

2βs r n + 2γ h 2 = log(n) + log(hs 2(4-r) n ) 2 1-r β h r + 2γ h 2 = log(n) + log(h 2r-7 s -2 n ).
We notice therefore that A 2 is dominating over A 3 , which is dominating over A 1 . The proposed (s n , h) ensure that the term A 2 is still the dominating term and gives the rate of convergence.

The case r = 2 is treated similarly, by taking derivatives we notice that the term A 2 and the term A 3 are of the same order and that the term A 1 is smaller than the others. 17) and (41). Combining ( 42) and (43) yields the announced result. The bound on W ρ is then a direct consequence of (12).

Proof of Proposition 2

Let φ(z) := θz 2 -1, where θ := 17) and (41). Combining (44) and (45) yields the announced result. The bound on W ρ is then a direct consequence of (12).

1 (1+ √ B) 2 is the solution in (0, 1) of (1 - √ θ) 2 = Bθ. When m + n ≤ φ(z), then s ≤ √ θz and z -s ≥ z(1 - √ θ) = √ Bθz. By (17), this means that l m,n (z) ≤ 1 π e -Bθz

Proof of Proposition 3

By (6) the term b 2 1 (n) can be bounded as follows

b 2 1 (n) = j+k≥N |ρ j,k | 2 ≤ j+k≥N exp(-2B(j + k) r/2 ).
Compare to the double integral and change to polar coordinates to get

b 2 1 (n) ≤ c 1 N 2-r/2 exp(-2BN r/2 ).

Proof of Proposition 4

To study the term b2 2 (n), we denote

F 1 [p ρ (•|φ)](t) := E ρ [e itX |Φ = φ] = W ρ (t cos φ, t sin φ),
the Fourier transform with respect to the first variable.

E[ρ η j,k ] = E[G j,k ( Y √ η , Φ)] = E[f η,δ j,k ( Y √ η )e -i(j-k)Φ ] = 1 π π 0 e -i(j-k)φ f η,δ j,k (y) √ ηp η ρ (y √ η|φ)dydφ = 1 π π 0 e -i(j-k)φ 1 2π f η,δ j,k (t)F 1 [ √ ηp η ρ (• √ η|φ)](t)dtdφ = 1 π π 0 e -i(j-k)φ 1 2π |t|≤1/δ f j,k (t)e γt 2 F 1 [p ρ (•|φ)](t) N η (t)dtdφ.
As N η (t) = e -γt 2 and by using the Cauchy-Schwarz inequality we have

E[ρ η j,k ] -ρ j,k 2 = 1 π π 0 e -i(j-k)φ 1 2π |t|>1/δ f j,k (t)F 1 [p ρ (•|φ)](t)dtdφ 2 ≤ 1 π π 0 1 2π |t|>1/δ f j,k (t) W ρ (t cos φ, t sin φ) dt 2 dφ. If 1/δ ≥ 2 √ N ≥ 2s with s = √ j + k + 1, then whenever t ≥ 1/δ we get by Lemma 2 | f j,k (t)| = π 2 |t|l j,k (t/2) ≤ π|t|e -1 4 (|t|-2s) 2 .
On the other hand, by Propositions 1 and 2 we have

| W ρ (t cos φ, t sin φ)| ≤ A( |t| 2 )e -β( |t| 2 ) r for β := B (1+ √
B) 2 in the case r = 2, or for arbitrary β < B and t large enough in the case 0 < r < 2. In both cases A is a polynom of degree 4 -r. We deduce the inequality

E[ρ η j,k ] -ρ j,k 2 ≤ C ∞ 1 δ t 5-r e -1 4 (t-2s) 2 -β2 -r t r dt 2 ≤ C( 1 δ ) 12-4r e -1 by Lemma 8 in [5], hence b 2 (n) 2 ≤ CN 2 ( 1 δ ) 12-4r e -1 2 ( 1 δ -2 √ N ) 2 -β2 1-r ( 1 δ ) r
which covers both cases in the proposition.

Proof of Proposition 5

Let us write σ

2 j,k (n) := E ρη j,k -E[ρ η j,k ] 2
. We bound it by

σ 2 j,k (n) = E 1 n n ℓ=1 G j,k ( Y ℓ √ η , Φ ℓ ) -E[G j,k ( Y ℓ √ η , Φ ℓ )] 2 = 1 n E G j,k ( Y √ η , Φ) -E[G j,k ( Y √ η , Φ)] 2 ≤ 1 n E G j,k ( Y √ η , Φ) 2 . ( 46
)
Proof of (31) For 0 < η ≤ 1/2, let us denote by K δ the function with the following

Fourier transform K δ (t) = I(|t| ≤ 1 δ )e γt 2 , then f η,δ j,k = f j,k (t) K δ (t) and we have σ 2 j,k (n) ≤ 1 n E f η,δ j,k ( Y √ η )e -i(j-k)Φ 2 ≤ 1 n E f j,k * K δ ( Y √ η ) 2 ≤ 1 n E f j,k (t)K δ ( Y √ η -t)dt 2 .
By using the Cauchy-Schwarz inequality

σ 2 j,k (n) ≤ 1 n |f j,k (t)| 2 dtE K δ ( Y √ η -t) 2 dt ≤ 1 n |f j,k (t)| 2 dtE 1 2π K δ (u)e -iu Y √ η 2 du ≤ 1 nπ f j,k 2 2 1/δ 0 e 2γu 2 du.
Then,

σ 2 (n) ≤ C nπ N -1 j+k=0 f j,k 2 2 ηδ 1 -η e 2γ δ 2 .
By Lemma 4 we have

N -1 j+k=0 f j,k 2 2 ≤ C 2 N 17/6 thus σ 2 (n) ≤ C 1 ηδN 17/6 nπ(1 -η) e 2γ δ 2 .
Proof of (32) and ( 33) By ( 28), for 1/2 < η ≤ 1,

σ 2 j,k (n) ≤ 1 n E f η j,k ( Y √ η )e -i(j-k)Φ 2 ≤ 1 nπ π 0 f η j,k (y) 2 √ ηp η ρ ( √ ηy|φ)dydφ ≤ 1 nπ f η j,k 2 ∞ For 1/2 < η < 1, by Lemma 5, σ 2 (n) ≤ C ∞ N 1/3 nπ e 8γN .
For η = 1, by Lemma 6

σ 2 j,k (n) ≤ 1 n π 0 |f j,k (x)| 2 p ρ (x, φ)dxdφ ≤ C n f j,k 2 2
hence by Lemma 4,

σ 2 (n) ≤ C C 2 N 17/6 n .

Proof of Proposition 6

It is easy to see that F E[ W η h ] (w) = W ρ (w)I( w ≤ 1/h).

We have, for n large enough s n ≥ z 0 and by ( 13) in the case r = 2. Now we write for the L 2 bias of our estimator :

W ρ 2
E[ W η h ] -W ρ 2 D(sn) ≤ E[ W η h ] -W ρ 2 2 = 1 (2π) 2 F E[ W η h ] -W ρ 2 2
= 1 (2π) 2 W ρ (w) , by the assumption on our class and ( 14), for 0 < r < 2. The case r = 2 is similar.

As for the variance of our estimator :

V W η h = E W η h -E W η h 2 D(sn) = 1 π 2 n E K η h [•, Φ] - Y √ η 2 D(sn) -E K n h [•, Φ] - Y √ η 2 D(sn) . (47)
On the one hand, by using two-dimensional Plancherel formula and the Fourier transform shown above, we get :

E K n h [•, Φ] - Y √ η 2 D(sn) ≤ π 2 |W ρ (w)| 2 dw ≤ π 2 . ( 48
)
In the last inequality we have used the fact that W ρ We replace in the second order moment, then as h → 0

E K η h [•, Φ] - Y √ η 2 D(sn) ≤ M (η)s 2 n 16γh exp 2γ h 2 (1 + o(1)). ( 49
)
The result about the variance of the estimator is obtained from (47)-(49).

2 dt(

 2 with s = √ k + j + 1). Because of Lemma 2, it is clear that the second integral is negligible in front of the first one, which we simply bound by 4s fk,j 2 ∞ .

1 B 2 /r z 2 . 2 ≤ 4 πBr z 4 -

 12224 Let φ(z) := (z -√ βz r/2 ) 2 -1. Since r < 2, for z larger than a certain z 0 (which depends only on β, B and r), it is true that φ(z) ≥ β It follows thate -Bφ(z) r/2 ≤ e -βz r (41) If m + n ≤ φ(z), then s ≤ 1 + φ(z) and z -s ≥ z -1 + φ(z) = √ βz r/2. By[START_REF] Krasikov | Inequalities for Laguerre polynomials[END_REF], this means that l m,n (z) ≤ 1 π e -βz r . So m+n≤φ(z) |ρ m,n |l m,n (z) ≤ Ae -βz r (42) for A := 1 π m,n e -B(m+n) r/2 . On the other hand, using Lemma 3 with ν := r/2, if φ(z) ≥ z 0 , m+n≥φ(z) |ρ m,n |l m,n (z) ≤ 4 πBr φ(z) 2-r/2 e -Bφ(z) r/r e -βz r (43) by (

  2 . So m+n≤φ(z) |ρ m,n |l m,n (z) ≤ Ae -Bθz 2 (44) for A := 1 π m,n e -B(m+n) . On the other hand, by Lemma 3, if φ(z) ≥ z 0 , m+n≥φ(z) |ρ m,n |l m,n (z) ≤ 2 πB φ(z)e -Bφ(z)

9 -

 9 2r exp(-2βt r )dtdφ ≤ C 1 s 10-3r n e -2βs r n , where β < B and for n large enough in the case 0 < r < 2, respectively β = B/(1 + √ B) 2

2 I

 2 r) e -2 1-r β w r dw ≤ C 2 h 3r-10 e -2 1-r β h r

2 2 =2ηh 2 ( 1 +

 221 Tr(ρ 2 ) ≤ 1 where ρ is the density matrix corresponding to the Wigner function W ρ . On the other hand, the dominant term in the variance will be given byE K η h [•, Φ] -•, φ) * N N η ([z, φ] -u)dφdzdu ≤ M (η)πs 2 n (K η h (u)) 2 du,using Lemma 6 below and the constant M (η) > 0 depending only on η, defined therein.Indeed, let us note that√ ηp η ρ (• √ η, φ) is the density of Y / √ η = X + (1 -η)/(2η)ε and let us call N N η the Gaussian density of the noise as normalized in this last equation. Let us first compute, by Plancherel formula, o(1)), as h → 0.

Lemma 6 .p ρ satisfies 0 ≤ π 0 p≤ π 0 p 0 W

 6000 For every ρ ∈ R(B, r) and 0 < η < 1, we have that the corresponding probabilitydensity ρ (•, φ) * N N η (x)dφ ≤ M (η), 0 ρ (x, φ)dφ ≤ Cfor all x ∈ R eventually depending on φ, where M (η) > 0 is a constant depending only on fixed η and C > 0.Démonstration. Indeed, using inverse Fourier transform and the fact that W ρ (w) ≤ 1 we get :π 0 p ρ (•, φ) * N N η (x)dφ ≤ F 1 [p ρ (•, φ)](t) • N N η (t)dtdφ ≤ c(η) π ρ (t cos φ, t sin φ) exp -t 2 (1 -η) 4η dtdφ ≤ c(η) 1 w W ρ (w) exp -w 2 (1 -η) 4η dw ≤ M (η),where c(η), M (η) are positive constants depending only on η ∈ (0, 1).

( 1 δ -2 √ N ) 2 -β2 1-r ( 1 δ ) r