B Salzmann
email: bjoern.salzmann@unifr.ch

A Pulkkinen

B Hildebrand

T Jaouen

S N Zhang

E Martino

Yazyev

Q Li

G Gu

H Berger

O V Yazyev

A Akrap

C Monney

Nature of

published or not. The documents may come

Over the past decades, investigations of the anomalous low-energy electronic properties of ZrTe5 have reached a wide array of conclusions. An open question is the growth method's impact on the stoichiometry of ZrTe5 samples, especially given the very small density of states near its chemical potential. Here we report on high resolution scanning tunneling microscopy and spectroscopy measurements performed on samples grown via different methods. Using density functional theory calculations, we identify the most prevalent types of atomic defects on the surface of ZrTe5, namely Te vacancies and intercalated Zr atoms. Finally, we precisely quantify their density and outline their role as ionized defects in the anomalous resistivity of this material.

Introduction. Zirconium Pentatelluride (ZrTe 5 ) is a material that first attracted attention for its anomalous transport properties at temperatures around 150 K, namely a peak in resistivity and a simultaneous change of the sign of the Hall coefficient [1,2]. Surprisingly, its chemical potential shifts from the valence band at room temperature to the conduction band at low temperature, suggesting an unknown intrinsic source of charge [3]. Although the material has been investigated for decades [4][5][6][7], no consensus has emerged on either the relevant details of its electronic structure, nor on its exact topological classification. Various recent studies using different measurement and simulation techniques have come to a wide range of conclusions concerning the nature of ZrTe 5 and the resistivity anomaly.

In recent studies, ZrTe 5 has been classified as a semiconductor [8] with a potential transition to a semimetal [9,10], a three-dimensional (3D) Dirac semimetal [11,12], a Weyl-semimetal [13], and additionally as both a weak [3] and a strong [14] topological insulator. Experiments contradicting both the classification as a 3D Dirac semimetal [15] and as a strong topological insulator [16,17] have also been published. One source of confusion appears to be variations between samples due to different growth techniques, namely Chemical Vapor Transport (CVT) and Flux methods. This can be addressed on a microscopic level with scanning tunneling microscopy (STM), as demonstrated recently for the transition metal dichalcogenide TiSe 2 [18]. However, previous STM and scanning tunneling spectroscopy (STS) studies on ZrTe 5 have focused on investigating the presence of surface states as a signature of topological properties [14,16,17], showing small and defect-free topographic images, and are relatively sparse in general due to difficulties in obtaining cleaved surfaces well suited for STM measurements. Although changes of the crystal lattice parameters induced by native defects have already been proposed [3,14], a direct comparison between samples grown with these two methods has so far not been explored with the full array of available techniques [8]. While charge localization at defect sites has been proposed as an origin of the seeming violation of charge conservation with changing temperature, no study of the specific defects present in ZrTe 5 samples has been published thus far [3].

Here, we present comparative STM and STS measurements on samples grown with both Flux and CVT methods. These allow insight into both structural and spectroscopic information up to atomic resolution. With the help of simulated STM images calculated within the density functional theory (DFT), we identify two types of repeating point-defects. Namely, Te vacancies and intercalated Zr atoms, the latter one being unique to CVT samples. We further describe a long-range chemical potential fluctuation in the electronic structure of the CVT samples that we identify as charge puddles. Finally, we propose that the identified native surface defects could explain the anomalous chemical potential of ZrTe 5 , since they potentially represent an important source of intrinsic doping.

Methods. Flux growth of ZrTe 5 occurs within a vacuum sealed ampoule, containing the elementary constituents of ZrTe 5 , with a Zr:Te ratio much smaller than 1:5 [11]. The Te melt begins to dissolve the solid Zr, acting as both the Flux as well as part of the desired product [19]. CVT growth involves an additional vapor component, Iodine, to serve as a transport agent for the constituents. The base components are placed in the ampoule in a ratio that is close or equal to the stoichiometric ratio of the desired product [8,16,20]. The STM and STS measurements were performed using an Omicron low temperature STM system in fixed current mode, with a bias applied to the sample. All measurements were made at a temperature of 4.5 K, with a fixed current of 0.2 nA unless otherwise noted. Samples were cleaved using scotch tape under vacuum (10 -8 mbar), before being transferred into ultra-high vacuum in the STM (<10 -11 mbar) and cooled down to the measurement temperature.

DFT calculations were performed using the VASP code [21][22][23][24][25] within the projector augmented wave method [26]. The exchange and correlation effects were treated within the generalized gradient approximation using the PBE functional [27]. The ZrTe 5 surface was modeled with a 2-layer supercell with 6×2 periodicity in the surface plane, corresponding to simulation cell dimensions 24.28 Å×27.67 Å. The atoms in the bottom ZrTe 5 layer were fixed to their bulk positions, and the top layer atoms were allowed to relax until the forces were converged to 0.02 eV/ Å. In the structural relaxations, the Brillouin zone was sampled using a Monkhorst-Pack grid of 3×3×1 k-points, and the kinetic energy cutoff was set to 400 eV. The simulated STM images were generated using the Tersoff-Hamann approach [28].

Results. Figure 1 shows a comparison between a small size atomically-resolved STM image and the atomic structure of ZrTe 5 , made of layers stacked along the c direction, see Fig. 1 CVT CVT Figure 2 shows large-scale STM images taken on both CVT and Flux samples with U bias = ±0.1 V. We observe two main deviations from a perfectly ordered sample. First, there is a long-range variation in the intensity on the surface of CVT samples (see Fig. 2(a) and (b)), with different regions tens of nanometers in size appearing brighter or darker. These variations are absent on Flux samples (see Fig. 2(c) and (d)). Second, there are two types of small and sharp defects that appear frequently in the STM images. We associate them with different types of atomic defects that occur either in or just below the surface layer. While the dark defects are found in both samples, the bright ones are only present in CVT samples. Similar STM images as those of Fig. 2 have been obtained on different surfaces of the same samples, as well as from other samples from the same batch, confirming the reproducibility of these results.

To identify the nature of these defects, we show in Fig. 3(a) a comparison between small-scale STM images and simulated STM images for different types of defects at U bias = -0.3 V, -0.1 V,0.3 V and 1.0 V [START_REF]STM and STS images, DFT calculations, resistivity measurements)[END_REF]. The first type of defect, shown in Fig. 3(a), appears as a reduction in intensity of the Zr chain at all bias voltages. We identify it as a Te vacancy in the top surface chains. The simulated images exhibit a strong agreement with the experimental data. They were obtained by removing one of the surface Te atoms from the top ZrTe 3 chain, as indicated in the top structural model (see Fig. 3(a)), and subsequently allowing the structure to relax. The relaxation is an important step, as it allows the remaining side Te atom to relax toward the middle of the surface chain, an energetically more favorable position [START_REF]STM and STS images, DFT calculations, resistivity measurements)[END_REF]. This results in a simulated defect with a more symmetric surface electronic density, in agreement with the experiment. The second type of defect, shown in Fig. 3(b), appears as a local increase in the intensity of the Zr chain for negative bias voltages and positive bias voltages up to around +0.4 V. At higher positive bias, the defect instead appears as a lowering of intensity at the same location. This specific behavior with bias voltage allows us to identify the defect as an intercalated Zr atom located directly below one of the Zr atoms in the chains.

We now turn back to the investigation of the long-range intensity variations observed at the surface of CVT samples, see Fig. 2(a) and (b). These intensity variations are clearly visible at U bias = +0.1 V and -0.1 V, but are absent at higher voltages [START_REF]STM and STS images, DFT calculations, resistivity measurements)[END_REF]. This indicates that they are not topographic in nature. In that case, they would appear at all voltages. Rather, they appear to be an electronic phenomenon confined to energies close to the chemical potential (µ). In the images taken at U bias = -0.1 V and +0.1 V, the intensities of the longrange variations invert. Brighter regions at one voltage appear darker at the other voltage, and vice versa. This further supports the interpretation of these variations being purely electronic.

These intensity variations can be understood considering the ZrTe 5 density of states (DOS) shown schematically in Fig. 4(a). The total DOS exhibits a minimum close to µ, and increases for both higher and lower energies within the relevant energy interval of ±0.1 V. Therefore, if the chemical potential shifts locally with respect to that minimum, STM measurements at positive and negative bias voltages will receive unequal contributions from the sample's integrated DOS. Whether the positive or negative voltage shows the region as brighter is then determined by the direction in which the chemical potential is shifted with respect to the local minimum in the DOS.

To support this scenario, we have performed STS measurements on CVT ZrTe 5 , which indeed display local variations in the DOS, as exemplified in Fig. 4(b) and (c). Figure 4(b) shows the spatially averaged dI/dV curves as a function of U bias within two different regions of the same STS map [START_REF]STM and STS images, DFT calculations, resistivity measurements)[END_REF]. Evidently, the DOS within the two regions are shifted with respect to each other, while otherwise showing the same features. Further, the overall shape of the DOS agrees with that seen in previous experiments [17]. Figure 4(c) shows a color map of several STS spectra taken along a line across the same STS map used for the average spectra in (b). It reveals a shift in the location of the plateau in the unoccupied states at U bias = 0.2 V over a distance of 6 nm.

Discussion. The bias voltage behavior of the longrange intensity variation and their profile in STS is very reminiscent of similar observations made in graphene [START_REF] Samaddar | [END_REF]31] and topological insulators like BiSbTeSe 2 [32]. In semimetals with a small density of states near the chemical potential or in small gap semiconductors, the electric potential of charge impurities-like ionized donors or acceptors in topological insulators-is weakly screened and gives rise to the so-called Coulomb disorder, i.e. a strong local band bending with fluctuating potential [33]. This local band bending can bring the valence or conduction band of a compensated semimetal or semiconductor across the chemical potential, leading to a local charge accumulation, called charge puddles. The spatial size of these charge puddles typically depends on the gap size, dielectric constant and density of charge impurities. In the case of ZrTe 5 , these charge puddles have typically a diameter of about 15 nm at 4.5 K in CVT samples. Interestingly, there are no charge puddles in our Flux samples. This can be explained by looking at their average STS spectra (away from charge puddles) in Fig. 4(d), showing a comparison between the spatially averaged DOS of a CVT and a Flux sample. At 4.5 K, while the chemical potential lies in the band gap at the surface of CVT samples, it is shifted inside the valence band at the surface of the Flux samples. Therefore, we expect more charge carriers at the surface of Flux samples, leading to an increased screening that suppresses the occurrence of charge puddles. While at first it seems compelling to associate the appearance of these charge puddles on the CVT samples with the influence of the bright atomic defects, no clear correlation between their location and the charge puddles can be seen in our STM images (see Fig. 2).

We now turn to the quantitative discussion of atomic defects. From a comparison between STM images and DFT calculations, Te vacancies (dark defects) were identified in both Flux and CVT samples, while Zr intercalated atoms (bright defects) were observed only in CVT samples (see Fig. 2). Te vacancies appear at a rate of (4.5±0.4) • 10 -3 and (5.2±0.4) • 10 -3 per unit cell when averaged over large imaged areas (more than 10'000 nm 2 each) for CVT samples and Flux samples, respectively. The average frequency of the intercalated Zr atoms on CVT samples is (2.1±0.3)•10 -3 per unit cell. Altogether, the Zr to Te ratio near the surface for both types of samples amounts to 1:4.98±0.01 and 1:4.995±0.001 in CVT and Flux samples respectively. We stress here that the tiny difference between these two ratios is significant, and is essentially coming from the presence of intercalated Zr atoms in CVT samples only. It is particularly interesting to relate our results to the work of Shahi et al. [8]. They notably studied both types of samples with bulk sensitive methods like x-ray diffraction and energy dispersive xray spectroscopy, finding a deficiency of Te in both CVT and Flux samples with a ratio of about 1:4.60 ±0.20 and 1:4.98 ±0.17, respectively. This is consistent with our conclusions, since this bulk study also evidences a better stoichiometry in Flux samples with our study providing direct evidence of the defects causing the deviation from a perfect stoichiometry.

What is the influence of the atomic defects on the electronic structure of ZrTe 5 ? As discussed previously, our STS measurements (see Fig. 4(d)) show a chemical potential difference between CVT and Flux grown samples. This difference corresponds to additional electron doping in CVT samples, consistent with the presence of intercalated Zr atoms (partially) giving their 4d and 5s electrons, in agreement with our calculations of the defectrelated doping of the two samples above. It is very likely that the atomic defect states common to both types of samples (Te vacancies) play an important role for transport properties, in particular for the anomalous chemical potential shift of ZrTe 5 . Indeed, states close to the chemical potential have mostly Te character [14], therefore such defect states could be easily activated even at low temperatures, leading to the change in dominant charge carrier type [3]. Finally, the electron doping due to Zr intercalated atoms accounts for the shift in temperature of the resistivity peak in CVT samples, in comparison to Flux samples and also the appearance of charge puddles in CVT samples due to a shift of their chemical potential in the gap. [START_REF]STM and STS images, DFT calculations, resistivity measurements)[END_REF].

Conclusion. In this letter, we have identified native defects for both ZrTe 5 samples grown by CVT and Flux methods and demonstrated their relevance to the bulk properties of the material. These include not only pointdefects, but also charge puddles appearing in CVT grown ZrTe 5 . The above results clearly show the importance of taking into account differences between samples produced by different growth methods. They also invite further investigation on the influence of growth method on the quality of the produced samples, and from there on the electronic and topological character of ZrTe 5 .

I. EXPERIMENTAL DATA

Figure SM1 shows a series of STM images at bias voltages from -1 V to +1 V taken at the same location on a CVT grown sample of ZrTe 5 . As mentioned in the main text, the long range intensity variations visible at ±0.1 V do not appear at higher voltages. Additionally, one recognizes two different types of recurring atomic defects from their intensity variation as a function of bias voltage, as described in the main text.

Figure SM2 shows a series of STM images at bias voltages from -0.3 V to +0.3 V taken at the same location on a Flux grown sample. Even at ±0.05 V bias voltage, no long range intensity variation is visible.

Figure SM3(a) and (b) shows two STS maps taken on CVT grown samples. Spectra from these maps were used in figure 4 in the main text. The map in Figure SM3 (a) shows a long range variation in intensity caused by a shifting of the DOS as described in the main text. The locations averaged for the dI/dV spectra in figure 4 Figure SM4 shows resistivity data for both sample types obtained on samples of the same batch as those used for the STM measurements. The resistivity peak for the Flux grown sample appears at lower temperatures than that of the CVT grown sample.

II. DFT CALCULATIONS

In addition to the Te vacancy and Zr intercalated atom presented in the main text, we also investigated a Zr vacancy as a potential defect occurring in ZrTe 5 in the same manner as the other two defects, by comparison with DFT simulations. Figure SM5 shows the comparison between these simulations and the experimental data at the same voltages as in the main text. We find that the correspondence between calculation and experiment is worse for this defect than for the intercalated Zr atom. Namely, the gap in the chain at the location of the defect seen in the simulations at 0.3 V and 1 V in figure SM5 does not appear in the measurements. At the same time, the relative brightness of the defect at -0.3 V and -0.1 V seen in the experiment is not well reproduced by the calculations. Most importantly, the calculations for the Zr vacancy do not reproduce the change from the defect being brighter than the surrounding chain to being darker when moving from 0.3 V to 1 V bias voltage.

Before simulating the STM images of the different defects, the crystal structure was allowed to relax in response to the introduced defect. Figure SM6 shows the resulting relaxed structures for each of the studied defects. The way the structure relaxes gives insight into the range of the lattice effects of the defect. In this context, the comparison between the Zr vacancy and the intercalated Zr again favours the Zr intercalation which shows a strong effect on the atomic position of Te atoms at longer distances away from the defect in the a direction which are absent for the Zr vacancy. This is best seen in the top view in figure SM6. When looking at the size of the defect visible in STM images, there is clearly a significant range of the defect along this direction. The defect size in the a direction at 0.3 V bias voltage where the defect appears largest is around 11 Å, close to three unit cells in the a direction (11.9 Å), which agrees well with the lattice distortion in the simulation of the Zr intercalation whereas the effects of the Zr vacancy are more short ranged.

Figure 1 .

 1 Figure 1. (a), (b) Crystal structure of ZrTe5 as seen from different directions. Te and Zr atoms are shown in gold and green respectively. (c) 8.3×10 nm 2 STM image of the (010) surface of ZrTe5 showing atomic resolution using a tunnelling current of 0.3 nA. (d) Enlarged and rotated region indicated in (c). The overlay shows the locations of atoms close to the surface. Below the image a single layer of ZrTe5 as seen along the chain direction is shown to visualize the height profile of the surface. The pairs of Te atoms aligned perpendicular to the Zr chains are shown in silver both in the overlay and the profile for visual clarity.

  (a) and (b). Each layer is comprised of chains of ZrTe 3 running along the a direction and connected by two Te atoms. Due to the relatively weak interlayer bonding, the material cleaves easily parallel to the a-c plane, exposing the (010) surface. Figures 1(c) and (d) show high-resolution STM images of the occupied states taken on a CVT sample at a bias voltage U bias = -0.3 V with atomic resolution. The comparison with the structural model in (d) highlights a very good correspondence between the expected position of the surface atoms and the image, allowing to identify the brightest parts of the STM image with the pairs of Te atoms atop the Zr chains.

Figure 2 .

 2 Figure 2. Large scale topographic STM images taken on samples of both types. (a), (b) STM images taken at the same location on a CVT ZrTe5 sample at different bias voltages. Both images are 50×50 nm 2 in size. (c), (d) Two STM images taken at the same location on a Flux ZrTe5 sample at the same bias voltages as (a) and (b). Both images are 50×50 nm 2 in size. See Supplementary Materials for STM images over a full set of bias voltages.

Figure 3 .

 3 Figure 3. Comparison between theoretical calculations and experimental images for the bright and dark surface defects. From top to bottom: Model of surface atomic structure in the vicinity of the defect with the location of the defect marked, model of defect as seen along the chain direction, simulated STM images of the defect at four different bias voltages and measured image of the defect at those voltages. (a) Missing Te atom from the top of a chain. The relaxation of the remaining Te atom toward the center of the chain is shown schematically. (b) A Zr atom intercalated below the top layer. The intercalated atom is shown in blue to distinguish it from the Zr atoms of the regular lattice.

Figure 4 .

 4 Figure 4. (a) Diagram showing the proposed origin of the long-range intensity variation (see text). (b) Comparison between STS spectra averaged over two areas of 3.75×3 nm 2 each, taken on a CVT sample showing opposite intensity variations. The spectra are shifted by 39 meV with respect to each other. (c) Map of STS spectra taken along a line across the transition between two areas with opposite intensity variation from the same STS map used in (b). The periodic pattern also visible in the spectra is due to the differences in the DOS on top of and between the Zr chains. (d) Comparison between STS spectra averaged over large areas on both CVT and Flux samples. The full STS maps used to obtain these data are shown in the Supplementary Materials.

  FigureSM1shows a series of STM images at bias voltages from -1 V to +1 V taken at the same location on a CVT grown sample of ZrTe 5 . As mentioned in the main text, the long range intensity variations visible at ±0.1 V do not appear at higher voltages. Additionally, one recognizes two different types of recurring atomic defects from their intensity variation as a function of bias voltage, as described in the main text.FigureSM2shows a series of STM images at bias voltages from -0.3 V to +0.3 V taken at the same location on a Flux grown sample. Even at ±0.05 V bias voltage, no long range intensity variation is visible.FigureSM3(a) and (b) shows two STS maps taken on CVT grown samples. Spectra from these maps were used in figure4in the main text. The map in FigureSM3(a) shows a long range variation in intensity caused by a shifting of the DOS as described in the main text. The locations averaged for the dI/dV spectra in figure4(b) of the main text are marked with boxes in their respective colors. The position of the line spectrum in figure 4(c) in the main text is shown in black. Figure SM3(b) shows the STS map whose spectra were averaged into the dI/dV curve representing the CVT sample in figure 4(d) of the main text. Figure SM3(c) shows an STS map taken on a Flux grown sample. The locations of the spectra averaged for the dI/dV curve representing the Flux sample in figure 4(d) of the main text are indicated by the green box. Figure SM3(d) contains a comparison between experimental dI/dV curves and calculated DOS from DFT calculations. The general behaviour of the DOS is reproduced correctly and the location of the Fermi level in the calculations is well aligned with the experimental data.FigureSM4shows resistivity data for both sample types obtained on samples of the same batch as those used for the STM measurements. The resistivity peak for the Flux grown sample appears at lower temperatures than that of the CVT grown sample.

Figure SM1 .

 SM1 Figure SM1. Full voltage series measured on a CVT grown sample from -1 V to +1 V taken at the location of the STM images in Figure 2(a) and (b) in the main text. All images show a 50x50 nm 2 area imaged with a fixed tunneling current of 0.2 nA and at 4.5 K.

Figure SM2 .

 SM2 Figure SM2. Full voltage series measured on a Flux grown sample from -0.3 V to +0.3 V taken at the location of the STM images in Figure 2(c) and (d) in the main text. All images show a 50x50 nm 2 area imaged with a fixed tunneling current of 0.2 nA and at 4.5 K.

Figure SM4 .

 SM4 Figure SM4. Resistivity as a function of temperature for samples grown by the two methods, both showing a peak in resistivity at different temperatures.

Figure SM5 .

 SM5 Figure SM5. Comparison between DFT calculation and experiment for a second proposed origin for the bright atomic defects seen on CVT grown samples. From left to right: structural model of the Zr vacancy defect, comparison between calculated and experimental images at different bias voltages.

Figure SM6 .

 SM6 Figure SM6. Relaxed structures in the regions around all three studied defects as seen in, from top to bottom, b, c and a directions. Atoms strongly displaced from their original positions are highlighted in different colors.

Département de Physique and Fribourg Center for Nanomaterials, Université de Fribourg, CH-1700 Fribourg, Switzerland

School of Engineering Science, LUT University, FI-53850, Lappeenranta, Finland

Université Rennes, CNRS, Institut de Physique de Rennes -UMR 6251, F-35000 Rennes, France

Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

National Centre for Computational Design and Discovery of Novel Materials MARVEL, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

Condensed Matter Physics and Materials Science Department, Brookhaven National Laboratory, Upton, New York 11973, USA

Institut de Physique des Nanostructures, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

Acknowledgments. This project was supported by the Swiss National Science Foundation (SNSF) Grant No. P00P2 170597. A.P. acknowledges the Osk. Huttunen Foundation for financial support, and CSC -IT Center for Science, Finland, for computational resources. B.S. acknowledges support from MaNEP through its Summer internship grants. S.N.Z. and O.V.Y. acknowledge support by NCCR Marvel. A.A. acknowledges funding by the SNSF, Grant No. P00P2 170544. Work at BNL is supported by US DOE, DE-SC0012704. We are very grateful to P. Aebi for fruitful discussions and for sharing with us his STM setup. We acknowledge M. Rumo and G. Kremer for help in preparing samples. Skillful technical assistance was provided by F. Bourqui, B. Hediger and O. Raetzo.

Supplemental Material

for Nature of native atomic defects in ZrTe 5 and their impact on the low-energy electronic structure


S Tricot

T Jaouen

D Sébilleau

P Schieffer

Energy dependence of interference phenomena in the forward-scattering regime of photoelectron diffraction

Keywords: X-ray photoelectron diffraction, forward scattering, energy-scan, high-order interference

X-ray photoelectron diffraction is a powerful spectroscopic technique in the direct legacy of C. S. Fadley that combines high sensitivity to the arrangement of atoms in crystals and element specificity providing unique fingerprints of selected atomic sites in matter. When used with kinetic energies between 500 eV and 1500 eV, its interpretation and description is based on the fact that the atomic scattering factors are strongly forward-peaked in such a way that low-angles scattering and backscattering events are respectively dominant and almost irrelevant in the photoemission process. In this paper we aim to demonstrate with the help of multiple-scattering simulations that energy scans of high-energy (500-1500 eV) forward-scattering photoelectron diffraction can provide valuable structural and chemical information about thin epitaxial films or stacking of two-dimensional materials.

Introduction

Angle-resolved photoemission of electrons from nearsurface atoms of solid crystals with kinetic energies ranging from 500 eV to 1500 eV provides efficient, element-specific, structural information about local atomic structures such as lattice parameters, composition and reconstruction of surfaces and interfaces, and adsorbate geometries (1; 2; 3; 4; 5). In this particular kinetic energy regime of the so-called X-ray Photoelectron Diffraction (XPD) technique, forward focusing effects indeed dominate the electron-atom scattering amplitudes in such a way that the photoemission intensity shows zero-order interference enhancements at emission angles corresponding to dense atomic directions [START_REF] Egelhoff | X-Ray photoelectron and auger electron forward scattering: A new tool for surface crystallography[END_REF]. Among the plenty of pioneering works carried out by C. S. Fadley in the field of photoemission, the development of the theoretical description of high-energy forward scattering (FS) XPD has taken a very special place throughout his career. That story begun by his first use of XPS FS to determine the internuclear axis or orientation of bonds of CO adsorbed on transition metal surfaces (6; 7) and continued over the years through his constant implication for the improvement of the simulation and interpretation of XPD for bulkand surface-atomic structures using multiple-scattering methods [START_REF] Sébilleau | Multiple Scattering Theory for Spectroscopies[END_REF].

In general, whether in the FS regime of XPD, in low-energy backscattering studies (<500 eV) [START_REF] Woodruff | Photoelectron diffraction: Early demonstrations and alternative modes[END_REF], or in hard X-ray photoelectron spectroscopy in which photoelectron diffraction patterns are dominated by Kikuchi bands (5-20 keV) [START_REF] Winkelmann | High-energy photoelectron diffraction: model calculations and future possibilities[END_REF], it is essential to perform multiple-scattering (MS) calculations and to vary the structural parameters for optimizing the theoryexperiment agreement (11; 12) and extracting the maximum amount of structural information. Such approach has been previously demonstrated to be highly powerful in so-called angleresolved photoemission extended fine structure studies [START_REF] Barton | Theory of angle-resolved photoemission extended fine structure[END_REF], that focused on energy scans of backscattering at low photoelectron kinetic energies [START_REF] Woodruff | Photoelectron diffraction: Early demonstrations and alternative modes[END_REF] for accessing the geometry of molecules, atoms or monolayers on crystal surfaces. We aim, in this paper, to demonstrate through MS calculations that in the high-energy (500-1500 eV) FS regime of photoelectron diffraction, the little-explored energy-scan data collection mode can specifically provide valuable structural and chemical information about thin epitaxial films or stacking of two-dimensional materials.

In the high-energy FS regime, fine structures related to highorder interferences between various scattering paths can combine to FS-related intensity enhancements in the angle scans ultimately leading to a breakdown of the forward scattering behavior characterized by intensity depressions instead of enhancements along the dense atomic directions [START_REF] Chassé | Breakdown of the forward scattering model in MgO(001)[END_REF]. We will specifically address the energy dependence of these interference phenomena in the FS regime of photoelectron diffraction. Through MS calculations performed on Si(001), a benchmark crystal in solid-state physics, we will demonstrate that exploiting interference phenomena in FS energy scans allow to solve tiny changes of atomic structures with an extremely high sensitivity at the sub-Å length scale.

Computational details

The multiple scattering spherical wave cluster calculations have been performed in the Rehr-Albers framework [START_REF] Rehr | Scattering-matrix formulation of curved-wave multiple-scattering theory: Application to x-ray-absorption fine structure[END_REF] by using the MsSpec program (16; 17; 18). The phase shifts have been calculated by using a complex Hedin-Lundqvist exchange and correlation potential (19; 20; 21) whose imaginary part is able to describe the finite electron mean free path in the final state, with a muffin-tin radius of 1.1 Å. Polar (θ) and azimuth (ϕ) angles are defined with respect to the normal and to the (100) emission plane of the (001) cluster surface, respectively. At angles θ = 0°and ϕ = 0°, the photon source defines an angle of 54.7°with respect to the surface normal aligned with the analyzer direction. Lattice vibrations were described by averaging over T -matrix elements and using isotropic mean-square displacements (MSD). For simulations done at room temperature (RT), the MSD was 6 × 10 -3 Å 2 . We systematically applied MSD values 20% larger for the surface atoms as they vibrate more than the bulk ones. The Rehr-Albers truncation order which controls the degree of sphericity introduced in the description of the basis functions used to expand the wave function around each atomic center was set to 2 for the Si 2s level. It is consistent with the rule of thumb given by the Fadley group [START_REF] Chen | Convergence and reliability of the Rehr-Albers formalism in multiple-scattering calculations of photoelectron diffraction[END_REF] according to which the truncation order is well approximated by l i +1, where l i is the angular momentum of the excited core state. The number of scattering paths to take into account in computing the final cross-section can dramatically increase when considering a large number of atoms and for large scattering orders. To avoid wasting computational resources in calculating scattering paths that weakly contribute to the final signal, we filtered out paths whose total lengths were longer than 30 Å (corresponding to 1.5-2 electron mean free path) and whose angles of any scattering events, except the last one, were greater than 40 • or 60 • .

Since we focus on high-energy FS regime of XPD, the effect of the deviations from bulk geometry of atoms located in the surface region is small [START_REF] Gewinner | X-ray photoelectron diffraction from Si(111): short versus longer range structural sensitivity[END_REF], so that we consider bulkterminated surfaces. To model the signal of a full substrate, a spherical cluster is built with a radius of 20 Å, truncated at the surface and one plane below the emitter's plane. The full signal is the sum of diffraction patterns computed for emitters in all planes from the surface plane down to the 15 th plane beneath the surface. This cluster shape allows to achieve MS calculations for relatively deep planes while keeping a reasonable number of atoms in the cluster (980 atoms). To model thin films, the same method as above was used but with a hemispherical cluster 15 Å in radius containing 363 atoms (9 planes).

Results and discussion

3.1. Energy-dependence of Si 2p XPD polar scans: Experiment vs. MS calculations In order to demonstrate how powerful MS calculations are for accurately describing both the intensity and shape of FS peaks in XPD polar scans, we first show, in Figure 1, a comparison between a high-resolution experimental XPD polar scan of the Si 2p core level of a Si(001) substrate recorded in the (100) emission plane [bottom curve with dotted symbols, Fig. 1b)] and MS calculations performed at 1382.28 eV of kinetic energy [simulated bottom curve in full line, Fig. 1b)]. Examples of the clusters used to compute the theoretical polar scans are shown in Fig. 1a) for emitter atoms in the 6 th (left hand-side) and 15 th planes (right-hand side) below the surface. The complete calculation was done for emitters in planes 1 to 15. The shape of each cluster was chosen to be spherical and truncated at the surface and planes deeper than one plane below the emitter were not considered. The multiple-scattering expansion of Figure 1: a) Side view of the cluster used to compute the full diffraction diagrams shown in (b) for planes 6 and 15. Complete calculations are done for emitters in planes 1 to 15. The shape of each cluster is spherical and truncated at the surface. Planes deeper than the plane below the emitter are not taken into account. b) Multiple scattering (MS) simulations up to scattering order 5 of XPD polar scans of the Si 2p core level for a Si(001) substrate in the (100) emission plane for selected photoelectron kinetic energies between 500 eV and 1500 eV (full lines). The lower curve at 1382.28 eV can be compared to the experimental data taken from Ref. [START_REF] Katayama | Elimination of X-Ray Photoelectron Diffraction Effect of on Si(100) for Accurate Determination of on SiO 2 Overlayer Thickness[END_REF] (symbol + line). The curves are normalized such as the main structure that spreads over 10-15 • around normal emission are normalized to [0, 1]. The curves are shifted in intensity for the sake of readability. Vertical dashes in the upper MS simulated polar scans are a guide to the eye of the energy shift of the high-order interferences-related fine structures around 20 • of polar angle. the photoelectron wave function was carried out up to the fifth order which we checked to be sufficient to achieve convergence. The Rehr-Albers truncation order was set to 2 and its value was decremented by 1 after each scattering event. In the excitation process, only the p → d channel was considered since the radial matrix element for the p → s channel is much smaller and can be dropped out to speed up the calculation. Paths with a scattering angle larger than 40°after a diffusion event were not considered. The calculations have been averaged over 49 directions filling a cone of 1 • half-angle. The MSD of the bulk Si atoms have been fixed to the RT value of 6 × 10 -3 Å 2 [START_REF] Mazur | Anisotropy of the mean-square displacements at the Si(001) -(2 × 1) surface[END_REF], and those of the surface atoms considered to be enhanced by a factor of 1.2. The experimental polar scan has been extracted from Ref. [START_REF] Katayama | Elimination of X-Ray Photoelectron Diffraction Effect of on Si(100) for Accurate Determination of on SiO 2 Overlayer Thickness[END_REF] and was acquired on hydrogen-terminated (001) p-type Si samples using a Physical Electronics 5400 XPS system. An angle resolution down to 1 • was obtained by adjusting the entrance slit of the analyzer. Non-monochromatic Al Kα Xray (1486.6 eV of photon energy) was used for excitation, and the incident light defined a fixed angle of 54.7 • with respect to the analyzer direction. As we can see, the agreement between the experimental polar scan and the MS-calculated one is remarkable. The positions, shapes and relative intensities of the various structures are very well reproduced. As expected in a forward focusing regime, the photoemission intensity shows, in a first approximation, enhancements along the [001] and [101] directions of the Si lattice. The intensity modulations between 0 • and 45 • of polar angles are nevertheless not only related to FS along dense atomic directions. In particular, we notice rather fine structures around 20 • of polar angle [marked by the vertical dashes Fig. 1b)], that gradually shift as the kinetic energy increases from 500 eV to 1500 eV within our theoretical XPD polar scans [full lines, Fig. 1b)]. This behavior, together with characteristic small amplitudes and widths (of the order of a few degrees), is typical of high-order interferences of the emitted and scattered electronic waves.

The energy dependence of the photoemission intensity along the [001] and [101] atomic directions then shows that highorder interferences fine structures are ubiquitous and combine to FS-related intensity enhancements. The most pronounced effect concerns the XPD peak along the [001] direction that shows, in both the experimental and the simulated XPD diagrams, high interferences-related fine structures at +/-3 • and a strong dependence on the kinetic energy of its shape exhibiting the well-known volcano shape at 500 eV, 900 eV, 1100 eV and 1500 eV of kinetic energies. This peculiar shape has been previously observed in other systems such as SiC, MgO, and GaN systems and has been shown to be induced by close neighboring atoms via inter-chain processes and to mainly depend on their scattering potentials and on the energy at which the scattering takes place (26; 27; 28). In the following, we will demonstrate that thanks to high angle resolutions as currently provided by most of photoemission set-up in synchrotron facilities, we can access detailed structural information of ultrathin layers by finely analyzing the energy dependence of interference phenomena acting on normal-emission FS peaks in so-called energy scans. Indeed, focusing on energy-scans of the normal-emission FS peak instead of large-polar angles ones (>30 • ), presents the great advantages to be much less affected by the apparatus function and angle averaging effects induced by possible mosaic twists of the overlayer crystallites resulting from the epitaxial growth process (29).

Calculated energy-scan of the Si(001) normal-emission FS peak

Let us now focus on the kinetic-energy dependence of a full XPD diagram around normal emission calculated for a Si cluster of 9 atomic planes representative of a Si(001) thin film of 1.1 nm in thickness. The MS calculations have been performed for photoelectron emission from the 2s orbital of Si rather than from the 2p one since it allows to significantly reduce the computational time without affecting the drawn conclusions. Scattering paths containing deflections by more than 60°after a scattering event were not into account in this calculation. Figure 2a) shows stereographic projections of (θ, ϕ) MS calculations performed up to the fourth scattering order for the 2s core level of Si at 4 selected kinetic energies in the FS regime of photodiffraction and with polar angles between normal emission and 15 • . No experimental averaging was taken into account at this stage and the MSD of the bulk and surface Si atoms have been fixed to the RT values of 6 × 10 -3 Å 2 and 7.2×10 -3 Å 2 , respectively. As observed for the Si substrate, the normal-emission XPD structure strongly depends on the kinetic energy of the photoelectrons. Around 800 eV and 1250 eV, it consists of a peak of ∼1 • wide, whereas at ∼ 900 eV and ∼1000 eV, a depression surrounded by a ring of enhanced intensity at polar angles of 3 • appears, forming the so-called volcano shape. To get a systematic view of the energy dependence of the interference phenomena acting on the normal-emission FS peak, we plot in Fig. 2b) the modulation function χ(E) of the Si 2s cross-section as a function of the kinetic energy (E), in the whole FS regime. The modulation function χ(E) is defined by:

χ(E) = I 0 (E) -I 0,direct (E) I 0,direct (E) , (1)
where I 0 (E) and I 0,direct are respectively the calculated, total and direct (diffraction-free), normal photoemission intensities at a given kinetic energy E. Interestingly, the modulation function χ(E) shows an overall oscillating behavior over a longenergy range reflecting interference phenomena between scattering paths at short-length scale in real-space. It especially shows strong variations between maxima values of ∼1.2 close to 800 eV and less than 0.1 around 1000 eV where almost no diffraction-induced intensity enhancement occurs. Superimposed to such a large amplitude and long-period oscillation, it exhibits, throughout the energy-scan, successive local maxima which reflect constructive and destructive interferences at larger length scales, to what has already been observed in the backscattering energy range of photodiffraction (100-400 eV) [START_REF] Woodruff | Photoelectron diffraction: Early demonstrations and alternative modes[END_REF]. Note that, whereas the low-energy region of XPD is usually dedicated to the study of adsorption of molecules on surfaces due to the favored 180 • scattering and the rather small mean-free path, in our case, the relevant scattering paths being mainly due to FS or scattering with small angles <90 • , energy scans of high-energy photoelectron diffraction are particularly well suited for probing the structural properties of thin epitaxial films or stacking of two-dimensional materials.

Origin of the long-period oscillating behavior in the normal-emission FS energy scan.

In order to get some insight on the nature of the shortlength scale interference phenomena responsible for the large amplitude and long-period oscillating behavior in the normalemission FS energy scan, we performed single scattering calculations over various simplified Si clusters in which atoms have been removed and looked for the involved scattering paths. Figure 3a) shows the calculated modulation function of the Si 2s normal-emission FS peak for a 5 atoms cluster. Obviously, the local maxima reflecting the long-range interferences are not tracked by our calculations, but, interestingly, the long-period oscillating behavior is already perfectly reproduced, with maxima at 714 eV and 1454 eV and a minimum at 1043 eV. [ 101], [011] and [0 11] directions surrounding by intensity ripples which correspond to high-order interferences with scatter-ing angles of 45 • . Looking at more closely between 0 • and 15 • of polar angles [Fig. 3b), right-hand side], we effectively see that intensity maxima are present at normal emission for 714 eV and 1454 eV of kinetic energy whereas at 1043 eV, we recognize the volcano shape. From these simulations, we can thus conclude that the maxima/minimum of the modulation function at 714 eV (1454 eV)/1043 eV partly originates from third and fourth orders constructive/destructive interferences between electronic waves scattered once by each of the four surrounding atoms [see white lines in the left panel of Fig. 3b)]. It is interesting to note that the normal-emission intensity modulation reflects the combined effect of 4 scattering atoms that are located relatively close to the emitting atom [secondnearest neighbors in the bulk Si( 001) lattice] at a moderated scattering angle of ∼45 • . The additional fact that the modulus of the scattering factor at this scattering is still significant, when compared to the normal-emission one, finally explains the enhancement of the interference effects in the [001] high-symmetry direction.

3.4. Impact of angle-resolution and temperature on normalemission FS energy-scans.

We now discuss the experimental requirements in terms of angle-resolution and temperature for accessing the maximum amount of structural information from the interference phenomena in the normal-emission FS energy scans. Focusing, first, on the impact of the finite size of the analyzer entrance hole on the intrinsic angle averaging of the measurement, we performed MS calculations of the Si 2s normal-emission modulation function as a function of the kinetic energy for half-angles of the detector acceptance cone ranging from 0 • to 3.5 • therefore simulating ultimate to low angle resolutions [Fig. 4a)]. As we can see, the modulation function shows a continuous amplitude decrease as the acceptance angle increases and does not show a strong dependence on kinetic energy anymore for angles higher than 3.5 • . This reflects through the disappearance of the highorder interference features in polar diagrams that only consist of FS-related intensity enhancements along the [001] and [101] and [111] atomic directions (results not shown). We also notice, in Fig. 4a), that even though the long-period oscillating behavior is rather robust with respect to the loss of angle-resolution, the local maxima at 650 eV, 790 eV and 1250 eV reflecting the long-range interferences already disappear for angles higher than 1.5 • . This demonstrates that further prospects of exploiting the energy-dependence of interference phenomena in the FS regime of XPD requires angle-resolution better than 1.5 • . Note that, with their angle-resolutions of tenths or even hundredths of degrees, the current two-dimensional, deflection-based hemispherical, or time-of-flight electron analyzers [START_REF] Ang | Time-resolved X-ray photoelectron diffraction using an angle-resolved time-of-flight electron analyzer[END_REF], now largely permit to perform energy-scan FS XPD with ultimate angle and energy resolutions on modern synchrotron facilities.

Another intrinsic mechanism of angle averaging in XPD comes from scattering of photoelectrons with vibrating atoms at finite temperature. To take into account the effect of such thermal vibrations, we varied, in our MS calculations, the amplitude of the uncorrelated and isotropic atomic motions described by the MSD. All calculations presented so far were performed using the experimental RT value for the bulk Si atoms and considering enhanced MSD at the surface (1.2 factor as compared to the bulk). We now consider the temperature dependence of the MSD that shows a quasi-constant value of 2.5 × 10 -3 Å 2 below 50 K, and non-linear and linear variations with a slope of 1.75 × 10 -5 Å 2 /K, between 50 K and 150 K, and above 150 K, respectively [START_REF] Mazur | Anisotropy of the mean-square displacements at the Si(001) -(2 × 1) surface[END_REF]. Figure 4b) shows the energy-scans of the normal-emission Si 2s modulation function for selected values of the Si MSD (the corresponding lattice temperatures are given in the caption) for a 9 atomic-planes cluster. The calculations were performed up to scattering order 4 with an analyzer acceptance cone of 0.5 • half-angle and for a constant Si lattice parameter of 5.43 Å . At the lowest temperature, the modulation function shows a large number of high-order interferences-related fine structures [see vertical dashes on Fig. 4b)]. As the temperature increases, the amplitude of the modulation function decreases, the finest structures progressively dis-appear, and only the long-period oscillating behavior remains at the highest temperatures. The modulation amplitude decreases by only 20% between 0 K and 282 K demonstrating that the lattice vibrations have only a minor effect on the FS energy-scans, in agreement with the experimental polar scans that show that high-order interference effects around the normal-emission FS peak can be readily observed at RT for a Si(001) surface given a high angle resolution [half angle < 1.5 • ] (24; 26; 27; 28).

3.5. Sensitivity of normal-emission FS energy-scans to lattice deformations. In this last section, we aim to conclude by illustrating how sensitive the normal-emission FS energy scans can be to tiny changes of lattice parameters and atomic structures as induced by strain. Here again, we focus on the energy-dependence of the normal-emission Si 2s modulation function computed for the 9 atomic-planes Si(001) cluster with an analyzer acceptance cone of 0.5 • half-angle and bulk and surfaces MSD values similar to those used in the first two subsections. As for Figure 4, these simulations correspond to scattering order 4 and scattering events included in a 60°half-angle cone. Figure 5a) shows the energy dependencies of normal-emission Si 2s modulation functions for lattice parameters, a 0 , of 5.43 Å , 5.33 Å , and 5.53 Å , corresponding to the true atomic structure of Si(001) at 300 K and to compressed and expanded lattices with relative variations of the lattice parameter, ∆a 0 /a 0 , of ±1.84%. Remarkably, the normal-emission FS energy-scan shows an extremely high sensitivity to such a small change in lattice parameter. Although its shape and amplitude are rather independent of a 0 , it shows energy-offsets, over large energy-scales, towards lower or higher kinetic energies depending on the positive or negative sign of ∆a 0 /a 0 and that depend on the kinetic energy, the offsets being larger for higher energies. In order to quantify the high-sensitivity of the normal-emission FS energy-scan to changes in lattice parameters, we introduce the parameter S lp (E) defined as:

S lp (E) = ∆χ(E) ∆a 0 , (2)
where ∆χ(E) is the variation of the modulation function at the kinetic energy E induced by a change ∆a 0 in the lattice parameter. Doing so for our Si(001) example, we find a maximum sensitivity of -4.1/Å at 900-910 eV of kinetic energy [see the vertical arrow in the Fig. 5a)]. Considering that, from an experimental point, modulation function variations ∆χ(E) of ±0.1 can be measured, the proposed method of energy-scanning normalemission FS peaks should allow to access lattice parameters changes as small as ∼ ±0.03 Å . Figure 5b) finally illustrates the sensitivity of normalemission FS energy-scans to in-plane biaxial lattice strains, as introduced during epitaxy for example. To this end, we compare the energy-dependence of the Si 2s normal-emission modulation function computed using the Si(001) lattice parameter a 0 of 5.43 Å to those calculated for thin films under tensile or compressive strains, therefore exhibiting in-plane biaxial deformations ϵ of ±1.0%. Both films are assumed to be uniformly strained and to be subject to tetragonalization of the Si lattice

) ε = -1 % ε = 0 % ε = + 1 % S i 2 s m o d u l a t i o n f u n c t i o n χ( E ) K i n e t i c
e n e r g y ( e V ) as calculated within the linear elasticity model using a Poisson coefficient of 0.28. As we could have expected, the impact of strain on the normal-emission FS energy-scan is very similar to that of the lattice parameter previously discussed. Here also, the strain leaves almost unchanged the shape and amplitude of the energy variations of the modulation functions, but introduces energy-offsets whose magnitude depends on both kinetic energies and ϵ and whose direction depends on the sign of the latter, solely [compressive (tensile) strains corresponding to negative (positive) energy-shifts]. In the same way as before, we introduce the parameter S strain (E) for quantifying the sensitivity of the normal-emission FS energy-scan to strain. It is defined by:

S strain (E) = ∆χ(E) ϵ , (3)
where ∆χ(E) is the variation of the modulation function at the kinetic energy E induced by a slight film deformation ϵ of ±0.25 %. In that case, we also find a huge sensitivity of normalemission FS energy-scans to strain with a maximum sensitivity of ∼-48 at ∼910-920 eV of kinetic energy [see the vertical arrow in the Fig. 5b)] that corresponds to lattice deformations of the order of ±0.2 % considering accessible modulation function variations ∆χ(E) of ±0.1.

Conclusion

In this paper, through MS calculations performed on a benchmark crystal, we have theoretically addressed the energy dependence of interference phenomena in the specific FS regime of photoelectron diffraction. We have shown that for a proper description of a normal-emission XPD peak measured with a high-angle resolution, both the usual FS-related zero-order and high-order interferences phenomena have to be taken into account. The latter being strongly dependent on kinetic energy, energy-scanning normal-emission XPD peaks in the FS regime gives rise to very-marked oscillations of the modulation function of the photoemission intensity that we demonstrated to be extremely sensitive to tiny changes of the crystal lattice therefore providing a highly powerful spectroscopic probe of atomic structures of ultrathin films and two-dimensional heterostructures. We believe that our study, in direct lineage of the C.S. Fadley legacy, may provide a new tool to the still very active XPD community. Indeed, given the high angle and energy resolutions as well as the broad tunable photon energy provided by most of photoemission set-up in modern synchrotron facilities, they should motivate further experimental prospects of exploiting normal-emission FS energy scans.

Figure 2 :

 2 Figure 2: a) Stereographic projections of the Si 2s photoemission intensity for photoelectron kinetic energies of 800, 900, 1000 and 1250 eV, for polar angles θ between 0 • and 15 • , and for azimuth angles ϕ ranging from 0 • to 360 • . The (100) emission plane corresponds to ϕ=0 • . The reported intensity corresponds to the modulation function of Eq. 1. The intensity distribution exhibits large variations especially at the normal emission. b) Si 2s normal-emission modulation function as a function of the kinetic energy. The dots represent the energy points of the stereographic projections of a).

Figure 3 :

 3 Figure 3: a) Single-scattering calculations of the Si 2s normal-emission intensity as a function of the kinetic energy. The inset shows the 5 atoms cluster used for the simulation. Long-energy range oscillations of the modulation function are well reproduced by this small subset of atoms. The dots mark the kinetic energies at which the stereographic projections of b) have been calculated. b) Left: stereographic projections of the Si 2s intensity calculated within the singlescattering approximation for photoelectron kinetic energies of 714 eV, 1043 eV and 1454 eV (from top to bottom). Polar angles vary between normal emission and 60 • and azimuth lies in the [0°, 360°] range. Right: same as left but for polar angles between normal emission and 15 • . Wavefronts are highlighted with white lines in the left panel of b) as a guide to the eye to identify the different interference orders leading to constructive or destructive intensity at normal emission (right panel).

2 S i 2

 22 s m o d u l a t i o n f u n c t i o n χ( E ) Kinetic energy (eV)

Figure 4 :

 4 Figure 4: a) Si 2s normal-emission modulation function as a function of the kinetic energy for different half-angles of the detector acceptance cone simulating ultimate (0 • ) to low (3.5 • ) angle resolutions. Even though long-range oscillations may still be visible at low resolution, fine structure features disappear above 1.5 • or 2 • . b) Same as a) but for different atoms MSD corresponding to temperatures of 0 K, 168 K, 282 K, 396 K, 511 K, 853 K and 1082 K. The acceptance of the analyzer is kept constant with an acceptance cone of 0.5 • half-angle.

Figure 5 :

 5 Figure 5: a) Si 2s normal-emission modulation function as a function of the kinetic energy for Si(001) lattice parameters a 0 of 5.43 Å and ±0.1 Å around this value. b) Same as a) but for in-plane biaxial tensile or compressive strains of 1% of the cubic cell. Both a) and b) show that energy scans in the 500-1500 eV range may provide an extremely high sensitivity to tiny structural changes.


domain: phys stat

text:
Gianni Jacucci

Louis Delloye

Davide Pierangeli

Mushegh Rafayelyan

Claudio Conti

Sylvain Gigan

Tunable spin-glass optical simulator based on multiple light scattering

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

Non-deterministic polynomial-time problems, also known as NP-problems, are tasks whose solve-time scales exponentially with the size of the input [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF]. Such problems appear in most domains whether it is economy, society, science, and cryptography [START_REF] Aaronson | [END_REF]. Importantly, NP-complete problems can be mapped onto each other, solving one solves all. Among these, finding the ground-state of an Ising spin-system represents one of the most known examples in physics [3][4][START_REF] Nishimori | Statistical Physics of Spin Glasses and In-formation Processing: An Introduction[END_REF][START_REF] Lucas | [END_REF][7].

Interestingly, Ising dynamics have been observed in numerous quantum [8,9] and classical systems such as random lasers [10,11], superconducting networks [12,13], polariton condensates [14,15], nonlinear wave propagation in disordered media [16], and degenerate optical parametric oscillators [17][18][19][20][21][22][23][24][25][START_REF] Hamerly | AI and Optical Data Sciences[END_REF][START_REF] Strinati | [END_REF][28][29]. Therefore, all aforementioned systems represent suitable platforms to implement NP-problem solvers. For example, coherent Ising-machines using degenerate optical parametric oscillators find approximate solutions for the ground-state of spin-systems, thanks to their inherent non-linearity, with great control on their couplings. Although optical parametric oscillators can now implement thousands of spins [24], their scalability remains limited by electronic circuits.

Recently, a second class of optical simulators based on linear optics was demonstrated by exploiting § These two authors contributed equally spatial light modulation. [30][31][32][33][34][35][36] In these Isingmachines-where the spins and their interactions are described by phase and amplitude modulation [30][31][32][33][34] or by a properly engineered gauge-transformation on the optical wavefront [35,36]-the ground-state is found by optimising the detected intensity via a recurrent electro-optical feedback. The advantage of this approach relies on the instantaneous calculation of the energy, as the necessary matrix products are encoded in light propagation and are, therefore, independent of the number of spins.

However, photonic Ising-machines based on wavefront-shaping using free-space propagation are limited to the implementation of a specific class of hamiltonians, known as Mattis models [37], where the couplings between spins are correlated and an exact solution of the energy ground-state exists. Recently, a more general class of hamiltonians, known as Sherrington-Kirkpatrick (SK) spin-glasses [START_REF] Mezard | Spin Glass Theory and Beyond[END_REF][START_REF] Parisi | Statistical field theory[END_REF], where the couplings are all-to-all and random was simulated by combining wavefront-shaping and light propagation inside a random medium [32,34]. Although this optical simulator is scalable and can outperform conventional hardware for a large number of spins [32], its applicability was limited to random all-to-all couplings with zero average, as determined by the transmission matrix of a random medium.

In this article, we present various strategies to tune the couplings of a spin-glass Ising simulator based on multiple light scattering. We experimentally demonstrate the ability to control the evolution of a system in a deterministic fashion by exploiting the knowledge of the transmission matrix of the scattering medium [START_REF] Popoff | [END_REF]. This allows us to observe the phase transition from a disordered to a fully-magnetised ground-state. Moreover, our work showcases the possibility of modifying the simulated Ising Hamiltonian by both introducing an external magnetic-field term and by controlling the number of degenerate ground-states.

II. RESULTS

A system of N spins is described by the following Hamiltonian [START_REF] Mezard | Spin Glass Theory and Beyond[END_REF][START_REF] Parisi | Statistical field theory[END_REF]:

H = - N i,j=1 J ij σ i σ j (1)
where σ {i,j} and J ij are the spins and their couplings, respectively. For all-to-all random couplings, the J ij are drawn from an i.i.d distribution, and finding the ground-state of Equation 1-which is referred to as Sherrington-Kirkpatrick model [START_REF] Nishimori | Statistical Physics of Spin Glasses and In-formation Processing: An Introduction[END_REF]-represents an NPhard problem. This specific problem can be mapped into light-propagation in disordered media, where the spins are encoded on a set of input modes or pixels with a binary phase state of 0 and π, corresponding to ±1 spin states, respectively [32,34]. This equivalence can be made explicit by writing the transmitted intensity after a multiply-scattering medium as [32]:

I T = N i,j=1 J ij σ i σ j = -H (2)
with J ij = -M m Re t im t jm , where m runs on the output modes and t im and t jm are transmission matrix (T ) elements. T is linked to the input modes on the SLM (E in ) and the output modes on the camera (E out ) via E out = T E in [START_REF] Popoff | [END_REF].

Equation 2 shows that maximising the light intensity over the selected output modes allows to retrieve the ground-state of the spin-system (details on the numerical framework used can be found in Appendix B) [32]. In particular, as the distribution of couplings is Gaussian and centred around zero, the ground-state corresponds to an ensemble of randomly oriented spins (mean magnetisation m = N i σ i = 0). As shown in Figure 1b, the J ij distribution can be modified by acting on the initial phases on the SLM. Indeed, such an operation leads to the following effective transmission matrix (T ): where ψ i represents the bias associated to each spin. Equation 3 rewritten in terms of matrix coefficients becomes:

T = T × diag(e iψi ) (3)
t im = t im e iψi (4)
therefore, the resulting couplings-for simplicity still indicated as J ij -where J ij = -M m Re t im t jm , can be modified by tailoring the bias.

An efficient way to tune the couplings is to exploit the knowledge of T . In particular, for a single output mode (M = 1), using as ψ i the phase-conjugation of the corresponding single-row

T (ψ i = -arg(t i ))
gives rise to a J ij distribution of only negative values. This results from the fact that after phase-conjugation T is real-valued. Figure 1b shows also that, consistently with having a negative average coupling J ij , the ground-state of the spins network is not anymore random but it exhibits a fully-magnetised state with degenerate orientations (m = ±1).

This result can also be understood in the framework of wavefront shaping [41]. The introduced bias is the phase mask that optimises the detected intensity, of which the binary spin mask represents an undesired perturbation. Therefore, the optimal intensity is retrieved when the spin mask is reduced to a uniform phase-which in terms of spins corresponds to a fullymagnetised ground-state.

The shift of the J ij can be further controlled by adding a noise mask (mask of random phases) to the phase-conjugation solution. This noise mask is modulated by an amplitude parameter (η) controlling the shift of the couplings, as shown in Figure 1b. A derivation of the effect of noise on the couplings is presented in Appendix C.

We tested experimentally our approach for controlling the couplings of the spin simulator (Figure 2a and Appendix A). A laser (Coherent Sapphire SF 532, λ = 532nm) is directed onto a reflective phase-only, liquid-crystal SLM (Meadowlark Optics HSP192-532, 1920×1152 pixels) divided into N macro-pixels (spins). The modulated light is projected on the objective back focal-plane (OBJ1, 10×, NA = 0.1) and focused on a scattering medium made of Teflon (DIFF) with 0.5mm thickness. The scattered light is collected by a second objective (OBJ2, 20×, NA = 0.4) and the transmitted intensity is detected by a CCD camera (Basler acA2040-55µm, 2048×1536 pixels). The spins and the desired bias are encoded by a spatial light modulator (SLM) in a phase pattern whose binary part is sequentially updated until the ground-state is reached. Note that for the optimisation any algorithm can be used, i.e., the setup is algorithm agnostic, as the advantage of the presented simulators resides in the parallel measurements of the energy. [32] Figure 2b shows the results for a single output mode (M = 1) on the detection camera (CCD). The groundstate agrees with the numerical predictions of Figure 1b, showing a high magnetisation. Figure 1b shows a typical ground-state averaged over thermal fluctuations [31]-see Figure A.1a for single realisations.

The difference in numerical (cf. Figure 1b, |m| = 1) and experimental results (cf. Figure 2b, |m| 0.8) can be explained by the contributions that each SLM pixel bears. Indeed, as shown in Figure 2c, the spins not aligned with the dominant orientation correspond to SLM pixels with very small T amplitudes-i.e., they marginally contribute to the output intensity and therefore cannot be detected over the optimisation. The limited sensitivity of the detection camera introduces an experimental noise term, which allows the spins with small T amplitudes to have any possible orientation. This is equivalent to state that the experimental Ising simulator has an effective, finite temperature [32]. showing that the light intensity retrieved after optimisation is close to that obtained from phase-conjugation even if the magnetisation does not reach its theoretical maximum.

Moreover, building upon T it is possible to drive the optical spin-system over a phase transition. By adding an artificial noise-controlled by the amplitude η (cf. Appendix C)-to the phase-conjugated mask enables to control the degree of magnetisation reached in the ground-state (cf. Figure 2d). This can also be understood in terms of couplings: as predicted in Figure 1b, introducing artificial noise effectively tunes the J ij distribution.

Figure 2 summarises the case of M = 1, which simplifies Equation 1 to the class of Mattis Hamiltonian where the J ij are correlated [32,34,37]. However, such systems are described by an Ising simulator where light propagates in free space [30,31,35,36]. To take advantage of the all-to-all random couplings introduced by the scattering medium is necessary to focus on a higher number of output modes (M > 1). Maximising the optical intensity on an area of multiple pixels represents a more complex wavefront shaping task as T has now a size of M × N and the number of accessible input is N. The corresponding spin couplings satisfy rank(J ij ) = N , which implies a NP optimisation problem for the ground-state search [32].

In the following, we exploited two approaches to choose the optimal bias for M > 1 modes: i) doing the phase-conjugation (PC) of the sum of the T rows; ii) performing the singular value decomposition (SVD) of T and use the argument of the first right singular vec- tor-i.e., the vector associated to the largest singular value [42][43][44]. More details on i) and ii) in Appendix A. The SVD results are more effective in inducing a magnetisation as the T size increases (cf. Figure 3 and Figure S3a for experimental and numerical results, respectively). Indeed, the SVD mask leads to a higher focus enhancement than the PC one, with a twofold enhancement for a square T (N = M = 256), and it produces a more marked shift of the J ij distribution (cf. inset of Figure 3 A decrease of the magnetisation in function of M can again be understood in terms of T amplitudes. Figure A.5b shows that for M = 256 a correlation between the alignment of the spins and their T amplitudes is observed. However, compared to M = 1 (c.f. Figure 2), the T amplitudes are more homogeneously distributed and decrease in value-meaning that more spins have a marginal, experimentally undetectable, contribution to the intensity.

In Figure 2 and Figure 3 we demonstrated that, by measuring T , a magnetised ground-state can be promoted in a spin-system with all-to-all couplings. Im- A constant external magnetic-field is implemented by fixing a section of the SLM. The angle associated to the fixed area determines the final state of the spins system. b) Tuning β allows transitioning between the two magnetised-states with opposite orientations. c) Adding a binary pattern (γi) to the bias is possible to determine the spatial distribution of the spins in the ground state. In particular, by using a chequerboard with binary values of 0 and π results in an antiferromagnetic-like ground state -whose fidelity to the bias was quantified in terms of staggered magnetisation (|m |). d) Changing the value of Γ) in the chequerboard part controls the probability of the obtaining a ferromagnetic (|m| = 1) or an antiferromagnetic-like (m = 0) ground state. Remarkably, when Γ = π/2, the Ising Hamiltonian presents three ground states -with a degenerate probability for the ferro-and antiferro-magnetic states.

portantly, the presented approach can also be used to generalise Equation 1 accounting for an external magnetic-field. In particular, as derived in Appendix C, by selecting a given portion of the SLM which will stay fixed to a specific phase (β) defined before the optimisation, i.e., by defining m SLM modes with phases b 1 , b 2 , ...b N = b, Equation 2 becomes:

I T = N i,j=1 J ij σ i σ j + B i σ i (5)
where B ∝ Re e iβ represents the added external field, defined from a fixed phase β ∈ [0, 2π).

This term represents a constant external magneticfield, driving the spins towards a specific sign of magnetisation. In the two limiting cases β = 0, π the machine is expected to settle a non-degenerate groundstate that depends on the value of β (cf. Figure 4a). Note that the bias is applied also onto the area of the SLM dedicated to the external field to effectively drive the evolution of the spins. Moreover, β tunes the intensity of the magnetic-field and therefore the transition between the two degenerate spins' orientations (cf. Figure 4b). The sharpness of which is defined by the ratio between the number of SLM pixels allocated to the spins and those describing the magnetic-field (see Figure A.2).

Remarkably, the ground-state found by the spinglass simulator can be further modified by adding a tailored, binary shift to the T bias (γ). By doing so, the system will converge to the defined shift pattern (cf. Figure 4c). Figure 4c shows that a ground-state resembling an antiferromagnetic system is obtained by defining γ k = (1 -e iπk )Γ, with Γ ∈ [0, 2π]. The fidelity of the ground-state-quantified in terms of staggered magnetisation (m ) [45]-is close to what was observed in Figure 1b. Moreover, Figure 4d shows that by varying the angle associated with the additional shift mask (Γ), the energy landscape of the system can be affected. Without any additional shift, the ground-state is fully-magnetised-as already shown in Figure 1b. When Γ = π 2 , the system has an additional equiprobable ground-state which becomes dominant when Γ = π. The form this shift and corresponding ground-state take can be extended to any arbitrary mask Figure A.4.

III. CONCLUSION

In conclusion, we reported a novel approach to tailor the Hamiltonian of all-to-all coupled spin-systems implemented on a spatial photonic Ising-machine. In contrast to previous works, our method controls the couplings by shifting their distribution. Such tuneability was obtained by exploiting the complex coefficients of the transmission matrix and tuning the number of output modes. In particular, knowing the transmission matrix allows to map a Sherrington-Kirkpatrick spin-glass model with tuneable couplings and an external field term through a co-propagating reference on the input wavefront. Finally, we demonstrated that greater control over the energy landscape-and hence the possible ground-states and their likelihood-can be obtained by introducing a tailored binary phase shift.

The accuracy of the Ising-machine, quantified in terms of magnetisation and transmission matrix amplitudes, can be improved by reducing the experimental noise of the setup. However, to target specific couplings a different approach is required. In particular, we foresee two possible routes: i) develop algorithms that takes advantage of the huge number of pixels of both the SLM and the camera to choose or group them to re-define a spin-system with the target couplings; ii) replace the random medium with an engineered system that maps the desired couplings. Although computationally and experimentally challenging, these strategies may unlock the ability to have a fully-tuneable spin-glass Hamiltonian. Moreover, the performance of the Ising-machine could be further improved by replac-ing the electro-optical feedback in the ground-state search with an all-optical, wavefront-shaping approach [46] The transmission matrix of the scattering medium was estimated as in [START_REF] Popoff | [END_REF]. In detail, each row of the TM can be reconstructed by monitoring how the intensity on a given CCD pixel changes when a phase modulation is applied to the input patterns on the SLM. The TM is sensitive to translations and rotations of the scattering medium as well as to the input and detection hard-ware. In this work, we define the stability as a variation within 10% of it original value. The time is long enough to run our experiments but for larger systems one would need more stable architectures.

The ground-state search is conducted sequentially by means of the recurrent digital feedback. Computation starts from a random configuration of N binary macro-pixels (spins) on the SLM. The measured intensity distribution determines the feedback signal. At iteration, a arbitrary batch of spins is randomly flipped if it increases the intensity at a chosen output mode. The batch size decreases over the optimisation procedure, starting from 12% of the pixels to a single pixel for the last ∼ 600 iterations. Note that our simulator is agnostic to the optimisation algorithm used.

Focusing on multiple output modes: phase conjugation and singular value decomposition

The ideal option to perfectly focus through a medium would be to be able to invert the transmission matrix. But this is extremely hard to do in practice we only use the phase-conjugated matrix. Each element of this matrix is defined as a complex element that has a magnitude of 1 and an argument equal to the opposite of that specific element in the initial transmission matrix, i.e. if θ is the argument of an element in the TM then -θ is the argument of the corresponding conjugated matrix.

For several output modes, one first needs to sum over the rows of the initial TM and then extract the arguments. Indeed, we only control N degrees of freedom thanks to the SLM but for several output modes the TM is N × m. Once this is done one can actually use the opposite angle to counteract the disorder medium. The larger the matrix the less easy it is to actually get a proper focus.

The second method is based on the singular value decomposition (SVD). It is defined as follows: any given matrix A ∈ C of dimensions m × n can be factorised as U ΣV † = A, where the columns of U (resp V ) represents the left (resp right) singular vectors and the diagonal matrix Σ represents the singular values ranked by descending order -i.e., the most representative vectors of A are the first ones of each matrix. In this particular nomenclature, the matrix V (resp U ) encodes for the output (resp input) modes. One can therefore use the opposite of the arguments of the first right singular vector (V here) to act as a phase conjugation substitute. Indeed, the first singular vector is the one with the biggest singular value and therefore is the one that represents the most the actual matrix.

Appendix B: Numerical details

The numerical model used in this work is a generalisation of what described in [32]. The optical SG is numerically simulated by forming N pixel blocks (SLM plane). The initial optical field has a constant amplitude, and its phase is a random configuration of N binary phases, φ i = 0, π. A transmission matrix T with random complex numbers is generated. At each iteration, a randomly selected single spin is flipped. The input phase updates if the intensity increased after the linear propagation of the field. The bias in the numerical framework is calculated as in the experiment by starting from the knowledge of T.

Numerical evaluation of I T corresponds to a measurement with a detector in a noiseless system. In general, within this scheme, ∼ 10N iterations are sufficient for a good convergence, i.e., when focus intensity reaches a plateau. All codes are implemented in MATLAB on an Intel processor with 14 cores running at 3.7 GHz and supported by 64 GB ram.

Appendix C: Details of the analytical calculations

Noise bias

Additionally to the context of having the spins superimposed with the SVD retrieved angles to compensate for the medium, one can define a random mask superimposed on top of this which reads as: e iθη , where θ is a random value in [0, 2π] and η ∈ [0, 1] is a global parameter which defines how strong this noise is compared to the SVD angles (η = 1 meaning the noise is a strong and therefore the situation is equivalent to not having retrieved any angles at all). Consequently, the intensity of one output mode m reads as:

I T =| i t im E i | 2 =   j t jm E j   i t im E i = i,j t im t jm E i E j = 1 2   2 i,j t im t jm E i E j  
where t im = A i e -i ψi are the TM coefficients and E i = e iφi e iψi e iηθi . The intensity and resulting couplings can therefore be rewritten as:

I T = i, j J ij σ i σ j (C1) J ij = -Re t im e iψi e iηθi t jm e -iψj e -iηθj (

C2)

External magnetic field

To derive the effective external field, one needs to separate the SLM into two parts, one portion encodes for the spins -optimised by the algorithm -and one portion encodes for the external field -fixed during optimisation. Therefore, considering the case where half the SLM is allocated to the spins σ i = e iφi as well as to the external field h j = e iβj , one can divide the usual sum into two terms :

I T = | i t im E i | 2 =   N/2 j=1 t jm E j + N j=N/2+1 t jm E j   ×   N/2 i=1 t im E i + N i=N/2+1 t im E i   = i,j t jm t im E j E i + i,j t jm t im E j E i + N/2 j=1 N i=N/2+1 t jm t im E j E i + N/2 i=1 N j=N/2+1 t jm t im E j E i
where E i = e φi+ψi or e βi+ψi depending on whether it corresponds to the spin area or the field area, respectively.

The first term is the one derived earlier which represents the classical Hamiltonian. The second one is simply a constant term based on the field value which will not change during optimisation -which will be neglected in the coming lines. Acting on the index order we obtain:

I T = i,j J ij σ i σ j + N/2 j=1 N i=N/2+1 t jm t im E j E i + N/2 i=1 N j=N/2+1 t jm t im E j E i = -H(σ) + N/2 i=1 N j=N/2+1 (t jm t im E j E i + t im t jm E i E j ) = -H(σ) + N/2 i=1 N j=N/2+1 (t jm t im h j e -iψj σ i e iψi + t im t jm h j e iψj σ i e -iψi )
By definition, σ i = e iφi and since the spins are binaryφ i = 0, π -then σ i = 1, -1 and therefore σ i = -σ i . Writing out the TM coefficients as t im = A i e i ψi and the field as h j = e iβj one gets :

I T = -H(σ) + N/2 i=1 N j=N/2+1
(A j e i ψj A i e -i ψi e -iβj e -iψj e iψi σ i + A j e -i ψj A i e i ψi e iβj e iψj e -iψi σ i ) = -H(σ)

+ N/2 i=1 N j=N/2+1
Re (e -i(ψj -ψj ) e i(ψi-ψi) e iβj A i A j σ i ( * )

Now, the ψ i,j are, by definition, the phase conjugated angles of the TM which implies that in case of an ideal phase conjugation ψ i = ψi . Considering the field is uniform, meaning β j = β then this last term can therefore be rewritten as:

( * ) = 2 i j Re e iβj A i A j σ i = 2 Re e iβ j A j i A i σ i = 2 Re e iβ   j A j   i A i σ i
We therefore obtain:

I T = i,j J ij σ i σ j + B i A i σ i
where B = 2 Re e iβ j A j describes the effective external field.

Figure 1 .

 1 Figure 1. Concept of a tuneable optical Ising simulator. a) Schematics of the Ising-machine based on multiple scattering. The spins are encoded in a binary phase on the spatial light modulator. The ground-state of the Ising Hamiltonian corresponds to the optimised binary phase mask maximising the transmitted intensity recorded by a camera.The all-to-all, random couplings of the spins are induced by multiple light scattering during the propagation in the disordered medium-described by the transmission matrix (T ). b) Measuring T determines the correct bias (ψi) to add to the initial phase mask to tune the distribution of couplings. In the case of fully-biased couplings (η = 0), the final state corresponds to a fully-magnetised degenerate ground-state.

Figure 2 .

 2 Figure 2. Tuneable experimental spin-glass simulator (M=1). a) Schematics of the experimental setup. Recurrent feedback from the measured intensity on a CCD updates the SLM configuration to reach the SG ground-state. b) Adding a bias to the spin-system (N = 256 spins) results in a more magnetised ground-state. Inset: intensity focus on the CCD plane as the ground-state is reached. c) The relatively low magnetisation compared to the simulation value can be understood in terms of the finite detection sensibility of the CCD. The spins which are not aligned with the dominant orientation (green (light gray) bars) are indeed those whose T-amplitudes, and therefore intensity-contribution, are the lowest. d) The degree of magnetisation of the ground-state can be readily tuned by adding an artificial noise-a mask of random phases of amplitude η-to the bias.

  Figure 2c is corroborated by Figure A.1b,

Figure 3 .

 3 Figure3. Tuneable experimental spin-glass simulator (M > 1). Magnetisation as a function of the number of output modes M on the CCD. Using the first singular value (SVD) of the T yields a magnetisation notably larger than with phase-conjugation (PC). This can be understood by comparing the relative focus intensity (δI) of the two methods: numerical simulations show that the SVD outperforms the PC when increasing M , with a two-fold enhancement for a square T (N = M = 256). The top inset displays the intensity distribution on the CCD for the two approaches-for a square T the SVD yields a less homogeneous intensity distribution compared to PC.

  and Figure A.3, respectively).

Figure 4 .

 4 Figure 4. Modifying the hamiltonian of an Ising simulator to control the ground state of a spin system. a) A constant external magnetic-field is implemented by fixing a section of the SLM. The angle associated to the fixed area determines the final state of the spins system. b) Tuning β allows transitioning between the two magnetised-states with opposite orientations. c) Adding a binary pattern (γi) to the bias is possible to determine the spatial distribution of the spins in the ground state. In particular, by using a chequerboard with binary values of 0 and π results in an antiferromagnetic-like ground state -whose fidelity to the bias was quantified in terms of staggered magnetisation (|m |). d) Changing the value of Γ) in the chequerboard part controls the probability of the obtaining a ferromagnetic (|m| = 1) or an antiferromagnetic-like (m = 0) ground state. Remarkably, when Γ = π/2, the Ising Hamiltonian presents three ground states -with a degenerate probability for the ferro-and antiferro-magnetic states.

2 I

 2 Figure A.1. a) Magnetisation values of different experimental realisations.Averaging the spins optimal mask allows to account for thermal fluctuations and obtain a reliable ground-state (as shown in Fig1b, main text) to compare to the corresponding theoretical zero-temperature solution. A noise-driven annealing approach does not improve the performance of the optimisation. b) Comparison between the focus intensity obtained via optimisation and from phaseconjugation of the TM for M=1. The relative focus intensity (δI) is close to unity, showing that the optimisation correctly reaches the highest focus intensity.

Figure A. 2 .Figure A. 3 .Figure A. 4 .Figure A. 5 .

 2345 Figure A.2. a) Simulation of the expected magnetisation in function of the number of output modes for both SVD and PC biases. b) TM amplitudes associated to the SLM spins when considering M=256. As in Figure 2, a correlation between the alignment of the spins and their TM amplitudes is observed. Compared to M=1 the TM amplitudes are more homogeneously distributed and significantly decrease in value, leading to a lower experimental magnetisation.

ACKNOWLEDGMENTS

We thank Dr. R. Pandya and Dr. J. Dong for useful discussions. G.J. and L.D. contributed equally to this work. G.J., L.D., D.P., M.R., C.C., and S.G. designed the project. L.D. and G.J. carried out experiments, numerical simulations and data analysis. G.J. and L.D. wrote the paper with contributions from all the authors. This project was funded by the European Research Council under the grant agreement No. 724473 (SMARTIES). D.P. and C.C. acknowledge funding from SAPIExcellence BE-FOR-ERC 2020 Grant, the QuantERA European Research Area Network (ERA-NET) Cofund [Grant No. 731473, the "Quantum Information Processing with Complex Media" (Quom-pleX) project], the Progetti di Ricerca di Interesse Nazionale (PRIN) 2017 program (Photonic Extreme Learning Machine -PELM project). M.R. acknowledges funding from Science Committee of Armenia under the grant agreements No. 20TTSG-1C008 and No. 21APP-1C016.

Appendix A: Experimental details

Transmission matrix calculation and ground-state search