ENAC Ikhlas Selmi

Test of GNSS receiver behavior in presence of Multiple Correlation Peaks Induced by Evil Waveform

Ikhlas Selmi is currently part of the SIGnal processing and NAVigation (SIGNAV) research axis of TELECOM team at ENAC. She graduated as engineer in Telecommunication and Signal Processing from ESIEE. She received her Ph.D. in 2013 from Telecom SudParis on mitigating GNSS interference between indoor and outdoor signals. She is currently working on GNSS signal distortion generated by the satellite payload.

Paul Thevenon graduated as electronic engineer from Ecole Centrale de Lille in 2004 and obtained in 2007 a research master at ISAE in space telecommunications. In 2010, he obtained a PhD degree in the signal processing laboratory of ENAC in Toulouse, France. From 2010 to 2013, he was employed by CNES, the French space agency, to supervise GNSS research activities and measurement campaigns. Since the July 2013, he is employed by ENAC as Assistant Professor. His current activities are GNSS signal processing, GNSS integrity monitoring and hybridization of GNSS with other sensors.

Christophe Macabiau graduated as electronics engineer in 1992 from the ENAC in Toulouse, France. Since 1994, he has been working on the application of satellite navigation techniques to civil aviation.

INTRODUCTION

In 1993, the first important GNSS signal distortion due to a payload failure was observed on GPS L1 C/A [START_REF] Enge | Detecting Anomalous Signals from GPS Satellites[END_REF]. This event prompted the GNSS community to standardize a kind of non-nominal signal distortion, named Evil Waveforms (EWF). A Threat Model (TM) has been proposed by ICAO for GPS L1 C/A to describe the possible distortions that can be observed on GPS signals [START_REF]Annex 10 -Aeronautical Telecommunications[END_REF]. The goal of the TM was to study three possible effects on GNSS correlation function computed by a receiver: dead zones (flat spot in the correlation function), distortions (correlation function asymmetry) and false peaks. The focus of this paper is on a consequence of the later effect, and in particular, on how a receiver will behave in presence of multiple correlation peaks induced by the combination of an Evil Waveform and its RF and tracking characteristics. The present study will focus on the analysis of Galileo E5a results obtained for Dual Frequency Multi-constellation (DFMC) user receiver to assess the effect of multiple peaks on the correlation function (that translate to multiple zero crossings on the discriminator function). But similar results are obtained with DFMC User Space applied to GPS L1, L5, Galileo E1c.

The paper has the following structure. First, the EWF ICAO Threat Models and the considered assumptions when simulating them with the tested parameters of the DFMC User Space are presented in section 1. Then, an example of multiple zero crossing situation in presence of EWF and the reaction of the receiver in this scenario are analyzed. In section 3, the considered receiver behavior models in presence of EWF with multiple zero crossing and the associated parameters (such as the pull-in area and C/N0 loss threshold) are defined. The analysis of the observed multiple zero crossing situations and the effect on the absolute tracking bias for Galileo E5a TM-A, B and C are given in section 4. In section 5, the sensitivity of these effects to the User Space and Threat Space is analyzed. Finally, the conclusion will summarize the observed effects of the multiple zero crossing situation and the recommendations to deal with the possible integrity threat induced by this situation.

SECTION 1 -DESCRIPTION OF THE EWF THREAT SPACE AND THE USER SPACE

EWF Threat Space

ICAO proposed three types of failures that could be related to payload functions, to the observed EWF event of 1993, and that would result in at least one of the three problematic effects on GPS L1 C/A receivers [START_REF]Annex 10 -Aeronautical Telecommunications[END_REF]:

• Threat Model A (TM-A) consists of the normal C/A code except that all the positive chips have a falling edge that leads or lags relative to the correct end-time for that chip. This TM is associated with a failure in the navigation data unit (NDU), the digital partition of a GPS or GLONASS satellite o This type of failure results in the creation of a flat zone at the top of the correlation function and an offset on the entire correlation function.

• Threat Model B (TM-B) introduces amplitude modulation and models the degradations in the analog section of the GPS or GLONASS satellite. More specifically, it consists of the output from a second order system when the nominal C/A code (baseband signal) is the input. TM-B assumes that the degraded satellite subsystem can be described as a linear system dominated by a pair of complex conjugate poles. These poles are located at: 𝜎± 𝑗2𝜋𝑓𝑑 where

• 𝜎 is the damping factor in unit of Mnepers/s

• 𝑓 𝑑 is the resonant frequency in unit of MHz.
o This type of failure results in the creation of distorted and false peaks.

• Threat Model C (TM-C) introduces both lead/lag and amplitude modulation. Specifically, it consists of outputs from a second order system when the C/A code signal at the input suffers from lead or lag. This waveform is a combination of the two effects described above.

These 3 threat models were endorsed by ICAO for GPS L1 C/A. As it can be understood, they depend upon 3 parameters:

• Δ representing the lead or lag relative to the correct end-time of the chip preceding the falling transition

• 𝜎 and 𝑓 𝑑 representing the second order system creating the amplitude modulation of the chip

The range of values of these parameters defines a so-called Treat Space (TS) that is representative of the feared events. The GPS L1 C/A TS is given in Table 1.

The proposed ICAO TSs for Galileo E1c and E5a signals are given in Table 2 [START_REF] Selmi | SBAS message table changes linked to spare definition[END_REF].

TM-A Galileo E1c

[-0.12 -0.12] --Galileo E5a

[-0.1 -0.1]

TM-B Galileo E1c - [0.1-63] [0.1-18] Galileo E5a - [0.1 -23] [0.1 -8] TM-C Galileo E1c [-0.12 -0.12] [0.1-63] [0.1-18] Galileo E5a [-0.1 -0.1] [0.1 -23] [0.1 -8]

User Space and considered assumption for EWF simulations

The considered receiver parameters at the reference station and the aircraft to simulate the EWF effect in the ENAC SW are given in Table 3.

A default assumption that is always considered in literature is that when the EWF appears, it is assumed that the signal instantaneously switches from nominal to an EWF status. An alternative assumption could be to consider the situation where the signal loses its nominal properties at 𝑡 𝐹𝑎𝑖𝑙𝑢𝑟𝑒 , then follows an unknown transient status and finally, takes on the EWF characteristics at 𝑡 𝐸𝑊𝐹 (>𝑡 𝐹𝑎𝑖𝑙𝑢𝑟𝑒 ). In the case of this second assumption, there are more chances that the user/reference station receiver could notice the strange signal change (through C/N0 variation and absence of significant useful signal leading to tracking loss). It seems that the first assumption may lead to worst scenario in term of integrity risk since the receiver may not detect the strange change in the signal caused by the presence of the EWF. The first assumption is considered in this study. It is also assumed that the smoothing filter has reached its steady state when the EWF occurs (at 𝑡 𝐸𝑊𝐹 ). This assumption is usually considered in literature when studying EWF effects.

In addition, when processing EWF, the observed effect (on pseudorange measurements and correlation function shape) remains constant since its instant of occurrence 𝑡 𝐸𝑊𝐹 and during the user and reference station process, notably during their smoothing filter convergence period. This assumption is also considered in this study and has an influence on the methodology used to determine the induced tracking error. The absolute tracking error due to the EWF is determined as the difference between the tracking error in presence of the EWF and the tracking error in nominal conditions. Both tracking errors are determined as the steady state bias (reached after filter convergence) associated with the zero crossing point of the DLL discriminator function. The EWF induced absolute and differential tracking errors are used to determine the probability of missed detection P md [START_REF] Rife | Formulation of a time-varying maximum allowable error for ground-based augmentation systems[END_REF] to respect in order to satisfy the DFMC SBAS integrity risk requirements for each EWF. When the derived P md requirement is different from 1, the processed EWF is considered as hazardous (can induce integrity risks). The worst EWF induced bias corresponds to the highest absolute bias which induces the lowest and most demanding P md .

Note that, in the present study, the tracking errors are determined only in the presence of line-of-sight signal distorted by the EWF. The effect of Multipath (non line of sight signal) combined to EWF is not considered in this study.

SECTION 2 -ILLUSTRATION OF MULTIPLE CORRELATION PEAKS AND THE RESULTING MULTIPLE ZERO CROSSING ON THE DISCRIMINATOR FUNCTION

For some extreme ICAO threats where multiple peaks are visible as shown in Figure 1, the expected behavior of GNSS receiver may change depending on the used multiple peaks mitigation techniques. Figure 1 shows the correlation function on the left and the E-L discriminator (with an E/L chip spacing of 0.1 chip) on the right obtained when Galileo E1c is distorted by TM-C with Δ = 0.09 µs, F d = 11 MHz and σ = 0.2 Mnepers/s and filtered with a 6 th order Butterworth filter at a bandwidth of 24 MHz. The obtained correlation function has multiple peaks corresponding to multiple zero crossings on the E-L discriminator that can be tracked by the receiver. When dealing with such heavily distorted correlation peaks, the tracking loop behavior is uncertain. However, it is certain that in this situation, the real receiver tracking loop behavior will be driven by two main parameters:

the tracking loop pull-in area: it corresponds to the area of the correlation function where the tracking loop will search for a peak and converge towards it. -the distorted peak amplitude: if the distortion creates a peak that is much lower than the nominal one, this will result in a decrease of the estimated signal to noise density ratio (C/N0), a degradation of the tracking accuracy or even a loss of lock.

More details on the considered definition for these parameters in this study are given in the next section.

Here is a tentative description of the chain of events following the occurrence of an EWF.

1. the correlation peak that was tracked is suddenly shifted due to the correlation function offset induced by the EWF. 2. If the shift is in the pull-in area of the delay tracking loop, then, the loop may be able to continue its tracking process. 3. If the shift is larger than the tracking loop pull-in area, then the loop is not able to have a sufficiently stable tracking of the distorted signal, and the tracking loop begins to diverge. The C/N0 estimated by the tracking loop, based on the tracked peak amplitude, will decrease and can trigger a loss of lock. 4. The tracking process may be stopped and a re-acquisition process could begin. 5. If so, the reacquired peak will likely be the one that has the highest amplitude as it is required by REQ-166

("The equipment shall acquire the main correlation peak for each Galileo ranging source used in the navigation solution") in ED259 [START_REF] Eurocae | ED-259-Minimum Operational Performance Standards for Galileo -Global Positioning System -Satellite-Based Augmentation System Airborne Equipment[END_REF], not particularly the one closest to the nominal peak. As can be seen, different possible receiver behaviors can occur after the occurrence of an EWF, and many depend on some design choices of the receiver. This situation may be further complicated by the implementation of peak ambiguity resolution techniques for BOC modulations or multipath mitigation techniques [START_REF] Margaria | A New Strategy to Mitigate Side-Peak Lock Errors in Tracking Multi-Peak Code Correlation Functions[END_REF]- [START_REF] Julien | A New Unambiguous BOC(n,n) Signal Tracking Technique[END_REF], that may react to the presence of multiple peaks or correlation function distortion.

In this paper, we have voluntarily chosen to model the receiver behavior with a simplified model and without taking into account complex tracking techniques related to multiple BOC correlation peak or multipath mitigation.

SECTION 3 -MODELS OF THE RECEIVER BEHAVIOR IN PRESENCE OF MULTIPLE CORRELATION PEAKS

In this study, two models of receiver behavior are considered:

• RBM 1: Model for receiver behavior with uninterrupted tracking • RBM 2: Model for receiver behavior with loss of tracking and re-acquisition process (including steps 4 and 5 depicted in red in Figure 2). It also corresponds to the first acquisition applied when the EWF occurs before the satellite starts to be visible to the user or reference station receiver (called rising satellite case).

The overall process to determine the position of the correlation peak in presence of potential multiple zero crossing is described in Figure 3.

Figure 3. Zero crossing research diagram for Model 1 & 2

For both models, the tracking loop pull-in area 𝐴 𝐸-𝐿 on the E-L discriminator function [START_REF] Kaplan | Satellite Signal Acquisition, Tracking, and Data Demodulation[END_REF] has to be defined for each receiver parameters of the User Space. In this study, the pull-in area is defined as the discriminator area with positive slope around the zero-crossing position obtained in the nominal condition (when no EWF is applied) taking into account only the RF/IF filtering effect. It is defined as:

𝐴 𝐸-𝐿 = [𝑃 𝑚𝑎𝑥 - 𝛼 𝐸 𝐶𝑆 2 , 𝑃 𝑚𝑎𝑥 + 𝛼 𝐿 𝐶𝑆 2 ] (𝑐ℎ𝑖𝑝𝑠)
Where

• 𝛼 𝐸 and 𝛼 𝐿 are factors computed for each tested user configuration to define the stable area with a positive slope within the filtered discriminator function.

• 𝑃 𝑚𝑎𝑥 is the position of the maximum peak that is determined:

o on the distorted correlation function for RBM 2
o on the nominal correlation function for RBM 1.

• 𝐶𝑆 is the chip spacing between the Early and Late replica Figure 4 shows an example of the nominal Galileo E1c correlation functions and discriminators filtered by a 6 th order Butterworth filter with 12 and 24 MHz bandwidths. The green areas on the discriminator function (and the red one on the correlation function) correspond to the pull-in areas obtained for each case of receiver configuration with an E/L chip spacing of 0.1 chip. It can be observed that the pull-in area width (proportional to α E + α L ) varies a lot (by a factor larger than 2 in the example) according to the user parameters. Depending on the number of zero crossings found in 𝐴 𝐸-𝐿 , a different process is applied:

• if there is no zero-crossing detected, the distorted signal cannot be tracked, and the associated DLL delay is not taken into account.

• if only one zero-crossing is detected, the associated delay is saved.

• if more than one zero-crossings are detected, the closest one found when starting the research from 𝑃 𝑚𝑎𝑥 and in the DLL direction search (right or left side of 𝑃 𝑚𝑎𝑥 , determined by the value of the discriminator after the introduction of the EWF) is selected. Depending on the tested RBM, the first selected peak corresponds to:

o the closest zero crossing to the nominal one for RBM 1 o the first stable zero crossing for RBM 2

If the associated correlation peak to this zero crossing has a correlation loss lower than the correlation loss threshold Δ(𝐶/𝑁 0 ) 𝑇 , the associated delay is saved as the DLL shift induced by the distortion. Otherwise, it is assumed that this correlation peak cannot be tracked by the DLL, and the next zero crossing is checked to compare its correlation peak amplitude loss to the threshold Δ(𝐶/𝑁 0 ) 𝑇 . This process is repeated until reaching a valid zero crossing (with peak correlation loss lower than the correlation loss threshold Δ(𝐶/𝑁 0 ) 𝑇 ) or going though all detected zero crossings within the pull-in area.

The considered threshold on the 𝐶/𝑁 0 loss, noted Δ(𝐶/𝑁 0 ) 𝑇 , is computed as: o Minimum correlation loss 𝐶 𝑙𝑜𝑠𝑠,𝑚𝑖𝑛 : 1 dB for wide band signals and receiver band larger than 24 MHz [START_REF] Novella | GNSS interference thresholds for L1 and L5 bands[END_REF].

Δ(𝐶/
o Noise power density 𝑁 0 . • (𝐶/𝑁 0 ) 𝑡𝑟𝑎𝑐𝑘𝑖𝑛𝑔 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑 is the minimum C/N0 tracked by the receiver.

The applied threshold Δ(𝐶/𝑁 0 ) 𝑇 for GPS L1 and L5Q and Galileo E1 and E5a are given in Table 4. -150 [START_REF]GPS SARPS JWGs/6-Flimsy 21[END_REF] -150 -200 [START_REF] Rtca | DO 292 -Assessment of Radio Frequency Interference Relevant to the GNSS L5/E5A Frequency Band[END_REF] 53 24 [START_REF] Novella | GNSS interference thresholds for L1 and L5 bands[END_REF], [START_REF] Rtca | DO 292 -Assessment of Radio Frequency Interference Relevant to the GNSS L5/E5A Frequency Band[END_REF] 33 Galileo E1(b+c) -151.45 [START_REF]ICAO NSP JWGs/6-Flimsy 23[END_REF] -154.45 -201.5 [START_REF] Rtca | DO 235B -Assessment of Radio frequency Interference relevant to the GNSS L1 frequency band[END_REF] 53.05 29 [START_REF] Julien | Investigation of Galileo E1 OS/SoL Acquisition, Tracking and Data Demodulation Thresholds for Civil Aviation[END_REF] 28.1 Galileo E5a(I+Q) -149.45 [START_REF]ICAO NSP JWGs/6-Flimsy 23[END_REF] -152.45 -200 [START_REF] Rtca | DO 292 -Assessment of Radio Frequency Interference Relevant to the GNSS L5/E5A Frequency Band[END_REF] 53.55 27 [START_REF] Novella | GNSS interference thresholds for L1 and L5 bands[END_REF], [START_REF] Rtca | DO 292 -Assessment of Radio Frequency Interference Relevant to the GNSS L5/E5A Frequency Band[END_REF] 30.6

SECTION 4 -OBSERVATION OF THE PRESENCE OF MULTIPLE CORRELATION PEAKS OVER THE EWF GALILEO E5A TS AND THE USER SPACE FOR THE TWO RECEIVER BEHAVIOR MODELS.

In this section, the distribution of the observed zero-crossings situations (depending on the number of the zero crossings found in the pull-in area) and the tracking biases associated to the selected zero crossing are presented for both RBM 1 & 2. Each sample of these distributions corresponds to a TM-A, B and C case processed by one tested receiver configuration (defined by the filter type, bandwidth and E/L chip spacing). Figure 5 and Figure 6 show the distribution of the number of obtained zero crossings (divided by the total number of simulated cases) with respect to the two tested receiver behavior models 1 & 2 (RBM 1 & 2) for TM-A, B and C applied to Galileo E5a signal. It appears that one zero crossing (unique zero crossing situation) is obtained for around 80% of the couples {TM-A, B and C distortion / receiver configuration} for both RBM 1 & 2. Among the remaining 20%, most of {TM-A, B and C distortion / receiver configuration} couples give a no-zero crossing situation and only few of them (less than 0.3%) give multiple zero crossing situation for both RBM 1 & 2.

Figure 5. Distribution of the number of the couples TM-A cases and associated receiver configurations with respect to the number of zero crossing obtained with the two receiver models 1 & 2 Figure 6. Distribution of the number of the couples TM-B (right) and TM-C (left) cases and associated receiver configurations with respect to the number of zero crossing obtained with the two receiver models 1 & 2

The effect of each EWF is evaluated by the induced absolute tracking bias estimated as the difference between the observed delays on the distorted and nominal correlation function processed by the same receiver configuration.

The absolute tracking bias difference between the bias obtained with RBM 1 and 2 for each TM-A, B and C and each tested receiver configuration is computed as:

𝛿𝑏 𝑖 = |𝑏_𝑢𝑠𝑒𝑟 RBM1 𝑖 -𝑏_𝑢𝑠𝑒𝑟 RBM2 𝑖 |
Where 𝑏_𝑢𝑠𝑒𝑟 RBM1 𝑖 and 𝑏_𝑢𝑠𝑒𝑟 RBM2 𝑖 are the absolute tracking bias obtained with user receiver configuration 𝑖 (defined in the User Space of Table 3) with RBM 1 and 2 respectively. A large 𝛿𝑏 𝑖 indicates a case where the EWF can potentially lead to an integrity risk for one user and potentially not for the other user. It is therefore an indicator of the sensitivity of an EWF impact to the RBM.

The distribution of the absolute bias difference 𝛿𝑏 𝑖 is given in Figure 7 and Figure 8 for TM-A, B andC. It shows that most of the absolute bias differences are around zero except very few TM-B and C cases where the difference 𝛿𝑏 𝑖 can go up to 85 meters. It means that for most couples {distortion / receiver configuration}, the same peak is tracked in RBM 1 and 2. However, for less than 3% of couples (TM-B and C distortion / receiver configuration), two different peaks (separated by 10 to 80 meters) are tracked with RBM 1 and 2.

5.a. Sensitivity of the Tracking Error obtained from both RBMs to the Threat Space

The question addressed in this section is to determine in which area of the TM-A, B and C TSs, Multiple Peaks situations are more likely to appear and may induce different tracking bias depending on the tested RBM. However, for both TM-B and TM-C, more receiver configurations are inducing multiple zero crossing with RBM2 (uninterrupted tracking). For RBM1, less than 15% of receiver configurations are in multiple zero crossing situation (for the concerned distortions) and it reaches 30% for RBM2. The analysis of Figure 11 shows that a subset of receiver configurations gives δb i > 1 cm for a different TS area from the one inducing multiple zero crossings. This shows that the second origin (b) described above is the cause of the large difference of absolute bias between the two considered RBMs. This observation confirms the importance and the impact of the notion of pull-in area when analyzing the multiple zero crossing situation.

For TM-B, less than 20% of the tested receiver configurations give different tracked peak in RBM 1 & 2. However, for TM-C, this percentage goes up to 24%. The question addressed in this section is which User Space parameters are more likely to experience Multiple Peaks situation that induces different tracking bias depending on the tested RBM. They correspond to the 15 or 30% of receiver configuration experiencing multiple zero crossing situation for some area of TM-B and TM-C TS as defined in section 5.a.

The distribution of the TM-B and TM-C cases inducing multiple zero crossing (divided by the total number of distortions) on Galileo E5a with respect to tested filter bandwidths and chip spacings is analyzed. Contrary to the Threat Space, it seems that all the user space is affected by multiple zero crossing situations. Hence, the 15 or 30% of receiver configurations inducing multiple zero crossing can take any value among the tested filter bandwidths and E/L chip spacings.

Figure 12 shows the distribution of the number of TM-B and TM-C, that induce two different tracked peaks when applying RBM 1 and 2 with a tracking bias difference 𝛿𝑏 𝑖 larger than 1 cm, with respect to the filter bandwidth and chip spacing (of Table 3). It appears that only some of the receiver configurations give two different peaks in RBM 1 & 2 with an associated bias difference larger than 1 cm. For example, with an E/L chip spacing of 1 chip and all filter bandwidths (from 12 to 22 MHz), the same correlation peak is always tracked with RBM 1 & 2 (since separated with less than 1 cm) for TM-B and TM-C (with Δ = ±0.12 µs).

Then, even if all tested receiver configurations can induce multiple zero crossing situation with some TM-B and C cases (for both RBM 1&2), only some of them can give different tracked peaks when applying RBM 1&2.

SECTION 6 -CONCLUSION

In this paper, only Galileo E5a results are presented to illustrate the observed impact of RBM1 and 2 and the multiple zero crossing induced by EWF situation. The results obtained for TM-A, B and C are summarized in Table 5 and Table 6. Table 5 gives the percentage of multiple zero crossing situations observed in the pull-in area for RBM 1 and 2 and Table 6 shows the maximum absolute difference δb i and the percentage of non-zero values for each TM. The presented results for the assessment of multiple zero crossings observed in pull-in area with RBM 1 (uninterrupted tracking) and 2 (reacquisition in case of loss of tracking) for TM-A, B and C applied to Galileo E5a signal show that only few EWF cases (less than 0.3%) are concerned. These cases are located around low σ and F d values for TM-B and TM-C. However, all tested User Space receiver configurations can induce this situation with multiple zero crossing in the pull-in area when applying both RBM 1&2.

The analysis of the effect of the tested receiver behaviors on the obtained tracking bias show that, for some TS area, two different peaks separated by a large shift (that can reach tens of meters for TM-C) are tracked with RBM 1 & 2. However, less than 5% (0% for TM-A) of couples {receiver configuration / distortion} induce two different peaks tracked by RMB1 and 2 and separated by a distance larger than 1 cm. It is then important to define the possible real receiver behaviors in presence of EWF in order to evaluate the worst (highest value inducing the lowest and most demanding P md requirement) absolute and differential bias induced by the TM-A, B and C. The computation of the worst induced absolute or differential bias is needed to characterize the effect of each distortion in order to evaluate the integrity risk that may be caused by its presence (through the evaluation of the associated missed detection probability P md that needs to be respected).

The tested receiver behavior RBM 1 & 2 correspond to two extreme receiver reactions to EWF occurrence. Different receiver behaviors (from RBM 1 & 2) may give different -and possibly worst (, ie with lower associated P md ) -results for EWF induced tracking bias. In order to have a more accurate simulation of the receiver behavior in presence of EWF and the tracked peak in multiple zero crossing situation, one solution would be to agree on the main peak definition for this scenario.

Another solution would be to test a real receiver behavior in presence of simulated EWF generated by an RF signal simulator. However, the obtained results would only reflect the behavior of the tested receiver and cannot be generalized to all user receivers.

Future work on the evaluation of the multipath and noise effect on the receiver behavior in presence of EWF would be necessary to carry out more accurate analysis of the receiver reaction in presence of EWF with multiple zero crossing situation. In addition, the considered assumptions to define the pull-in area need to be further consolidated for a more realistic analysis of multiple zero crossings issue.

Figure 1 .

 1 Figure 1. Example of correlation with multiple peaks (left) and the associated discriminator E-L with multiple zero crossings (right) for Galileo E1c signal distorted by TM-C with 𝜎 = 0.2 Mnepers, 𝐹 𝑑 = 11 MHz and 𝛥 = 0.09 µ𝑠

Figure 2

 2 Figure 2 represents the diagram of the chain of event expected at the receiver.

Figure 2 .

 2 Figure 2. Expected chain of event in the GNSS receiver following EWF occurrence.

Figure 4 .

 4 Figure 4. Examples of pull-in area defined for Galileo E1c E-L discriminator filtered by 6 th order Butterworth filter and bandwidth of 12 (Left) and 24 (right) MHz and an E/L spacing of 0.1 chip.

Figure 7 . 2 Figure 8 . 2 SECTION 5 -

 72825 Figure 7. Distribution of the difference on the absolute tracking bias obtained with TM-A and the tested user receiver configuration obtained with the two receiver models 1 & 2

Figure 9 and

 9 Figure9and Figure10present the distribution of the number of receiver configurations inducing multiple zero crossings (divided by the total number of simulated receiver configurations) with respect to TM-B parameters (F d and σ ) for TM-B and TM-C with Δ = ±0.12 µs applied to Galileo E5a. To illustrate the observed results for TM-C distortions, only one example of |Δ| value is given in this section and in the following ones. It corresponds to the TM-

Figure 9 .

 9 Figure 9. Distribution of the number of receiver configurations inducing multiple zero crossing with respect to TM-B parameters for the two RBM 1 (left) & 2 (right) for TM-B Figure 11 represents the number (divided by the total number of tested receiver configurations) of receiver configurations giving a difference on the absolute bias 𝛿𝑏 𝑖 larger than 1 cm with respect to σ and F d parameters for TM-B and TM-C (with Δ = ±0.12µ𝑠). It shows the TS area where different peaks (located at a distance larger than 1 cm) are tracked with RBM 1 and 2. It appears that the TS areas experiencing multiple zero crossing situation in the pull-in area with RBM 1 & 2 (shown in Figure 9 and Figure 10) are not the one where different peaks are tracked with RBM 1 & 2. The observation of two distinct RBM tracking a different peak can be due to two origins: a. If the same (or almost) pull-in area is considered for RBM1 and RBM2 and the EWF induces multiple peaks in this pull-in area, each RBM can track a different peak b.the EWF causes a loss of tracking lock for RBM 2, which triggers a re-acquisition. The re-acquired peak is then searched at a different location from RBM 1. In this case the pull in areas considered can be distinct since they are centered around 𝑃 𝑚𝑎𝑥 that is potentially different between RBM 1 and 2.

Figure 10 .

 10 Figure 10. Distribution of the number of receiver configurations inducing multiple zero crossing with respect to TM-B parameters for the two RBM 1 (left) & 2 (right) for TM-C with 𝛥 = ±0.12 µ𝑠

Figure 11 .

 11 Figure 11. Distribution of the number of receiver configurations inducing a difference larger than 1 cm on the absolute tracked peaks biases obtained with RBM 1 & 2 with respect to TM-B parameters for TM-B (left) and TM-C with 𝛥 = ±0.12 µ𝑠 (right) 5.b Sensitivity of the Tracking Error induced in the multiple zero crossing situation to the User Space

Figure 12 .

 12 Figure 12. Distribution of the number of receiver configurations inducing a difference larger than 1 cm on the absolute tracked peaks biases obtained with RBM 1 & 2 with respect to US parameters: E/L chip spacing and filter bandwidths for TM-B (left) andTM-C with 𝛥 = -0.12 µ𝑠 (right).

Table 1

 1

		. ICAO-standardized GPS L1 C/A (ICAO 2006) and L5 Threat Space
	ICAO parameters	∆(Chips)	𝝈(Mnepers/s)	𝒇 𝒅 (MHz)
	TM-A	[-0.12 -0.12]	-	-
	TM-B	-	[0.8-8.8]	[4-17]
	TM-C	[-0.12 -0.12]	[0.8-8.8]	[7.3-13]

Table 2 .

 2 ICAO Threat Space for Galileo E1c and E5a

	ICAO parameters	∆(µ𝒔)	𝝈(Mnepers/s)	𝒇 𝒅 (MHz)

Table 3 .

 3 Parameters of Reference station and User Space applied to GPS L1/L5 and Galileo E1c/ E5a for DFMC SBAS users

Galileo E1c signal (𝐶𝐵𝑂C

		reference	user	reference	user
	Tracking technique	EML (BOC(1.1) local replica for Galileo E1c)	EML (BOC(1.1) local replica for Galileo E1c)	EML (𝐵𝑃𝑆𝐾(10) local replica)	EML (𝐵𝑃𝑆𝐾(10) local replica)
	Correlator spacing	0.08, 0.1 and 0.12 chip	0.9, 1 and 1.1 chip
	Pre-correlation bandwidth (double-sided)	24 MHz	12,14,16,18,20, 22,24 MHz	24 MHz	12,14,16,18,20, 22,24 MHz
	Equivalent reception filter	one filter is tested (6th-order Butterworth) to estimate the reference error as the minimum on all the correlator spacing values	4 filters are tested (6th-order Butterworth, 0-group delay resonator, 150 ns maximum differential group delay resonator, 150 ns maximum differential group delay 6th-order Butterworth)	one filter is tested (6th-order Butterworth) to estimate the reference error as the minimum on all the correlator spacing values	4 filters are tested (6th-order Butterworth, 0-group delay resonator, 150 ns maximum differential group delay resonator, 150 ns maximum differential group delay 6th-order Butterworth)

(6.1)) and GPS L1 C/A

Galileo E5a and GPS L5 signal 𝐵𝑃𝑆𝐾(

10

))

  𝑟𝑒𝑐𝑒𝑖𝑣𝑒𝑟 is the margin taken on the minimum 𝐶/𝑁 0 tracked by the receiver. It is set to 2 dB.• 𝛿Δ(𝐶/𝑁 0 ) is the uncertainty margin taken on Δ𝐶/𝑁 0 . It is set to 2 dB.• (𝐶/𝑁 0 ) 𝑚𝑎𝑥 is the maximum 𝐶/𝑁 0 received at the user receiver. It is computed as:(𝐶/𝑁 0 ) 𝑚𝑎𝑥 = 𝑃 𝑚𝑎𝑥 + 𝐺 𝑎𝑛𝑡,𝑚𝑎𝑥 -𝐶 𝑙𝑜𝑠𝑠,𝑚𝑖𝑛 -𝑁 0 o Maximal signal power 𝑃 𝑚𝑎𝑥 received at the receiver front end ▪ For DFMC SBAS user, only maximal signal power at the pilot channel is considered as given in Table4by 𝑃 𝑚𝑎𝑥 𝑝𝑖𝑙𝑜𝑡 .

	𝑁 0 ) 𝑇 = (𝐶/𝑁 0 ) 𝑚𝑎𝑥 -((𝐶/𝑁 0 ) 𝑡𝑟𝑎𝑐𝑘𝑖𝑛𝑔 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑 -𝛿(𝐶/𝑁 0 ) 𝑟𝑒𝑐𝑒𝑖𝑣𝑒𝑟 ) + 𝛿Δ(𝐶/𝑁 0 )
	Where
	• 𝛿(𝐶/𝑁 0 )

o Maximum antenna gain 𝐺 𝑎𝑛𝑡,𝑚𝑎𝑥 : 4 dB for L1 & L5

[START_REF]ICAO NSP JWGs/6-Flimsy 13[END_REF]

.

Table 4

 4

			𝛥𝐶/𝑁 0 threshold for GPS L1/L5Q and Galileo E1 and E5a
	Signal	𝑷 𝒎𝒂𝒙 (dBW)	𝑷 𝒎𝒂𝒙 𝒑𝒊𝒍𝒐𝒕 (dBW)	𝑵 𝟎 (dBW/Hz)	𝑪/𝑵 𝟎 𝒎𝒂𝒙 (dBHz)	(𝑪/𝑵 𝟎 ) 𝒕𝒓𝒂𝒄𝒌𝒊𝒏𝒈 𝒕𝒉𝒓𝒆𝒔𝒉𝒐𝒍𝒅 (dBHz)	𝜟(𝑪/𝑵 𝟎 ) 𝑻 (dB)
	GPS L1	-153 [12]	-	-201.5 [13]	51.5	29 [14]	26.5
	GPS L5Q

Table 5 .

 5 Percentage

		of the couples TM (-A, B and C) and receiver configurations inducing multiple zero crossing in the pull in area
			for RBM 1 & 2 for Galileo E5a
	ICAO TM	TM-A	TM-B	TM-C
	Galileo E5a	RBM 1	0.1%	0.1%	0.19%
	DFMC	RBM 2	0.1%	0.3%	0.02%

Table 6 .

 6 Maximum 𝛿𝑏 𝑖 and percentage of 𝛿𝑏 𝑖 ≠ 0 obtained for TM-A, B and C for Galileo E5a

		ICAO TM	TM-A	TM-B	TM-C
	Galileo E5a	Max(𝛿𝑏 𝑖 )	0 m	80 m	82 m
	DFMC	%age of 𝛿𝑏 𝑖 ≠ 0	0%	4.39%	4.38%

Acknowledgement

The authors would like to thank Mikael Mabilleau (EUSPA: European Union Agency for the Space Program) and Denis Bouvet (Thales AVS France) for providing relevant feedback on the study presented in this paper.

Jaron Samson works as a Navigation Engineer at ESA's Technical Centre (ESTEC), being involved in EGNOS, GNSS Evolutions, and R&D related to navigation. From 2012 until 2021, Jaron has worked as a System Engineer at ESA's EGNOS Project Office in Toulouse (France).


Pin-Cheng Lin

Renan Villarreal
email: renan.villarreal@kuleuven.be

Simona Achilli

Harsh Bana

Maya N Nair

Antonio Tejeda

Ken Verguts

Stefan De Gendt

Manuel Auge

Hans Hofsäss

Steven De Feyter

Giovanni Di Santo

Luca Petaccia

Steven Brems

Guido Fratesi

Lino M C Pereira
email: lino.pereira@kuleuven.be

Doping Graphene with Substitutional Mn

Keywords: graphene, doping, manganese, electronic structure, magnetism, ion implantation

We report the incorporation of substitutional Mn atoms in high-quality, epitaxial 1 graphene on Cu(111), using ultra-low energy ion implantation. We characterize in detail the atomic structure of substitutional Mn in a single carbon vacancy and quantify its concentration. In particular, we are able to determine the position of substitutional Mn atoms with respect to the Moiré superstructure (i.e. local graphene-Cu stacking symmetry) and to the carbon sublattice; in the out-of-plane direction, substitutional Mn atoms are found to be slightly displaced towards the Cu surface, i.e. effectively underneath the graphene layer. Regarding electronic properties, we show that graphene doped with substitutional Mn to a concentration of the order of 0.04%, with negligible structural disorder (other than the Mn substitution), retains the Dirac-like band structure of pristine graphene on Cu(111), making it an ideal system in which to study the interplay between local magnetic moments and Dirac electrons. Our work also establishes that ultra-low energy ion implantation is suited for substitutional magnetic doping of graphene; given the flexibility, reproducibility and scalability inherent to ion implantation, our work creates numerous opportunities for research on magnetic functionalization of graphene and other 2D materials.

layers [START_REF] Ajejas | Unraveling Dzyaloshinskii-Moriya Interaction and Chiral Nature of Graphene/Cobalt Interface[END_REF] and magnetic molecules. [START_REF] Candini | Graphene Spintronic Devices with Molecular Nanomagnets[END_REF] However, these approaches present various limitations in terms of stability against ambient temperature and atmosphere (e.g. segregation and oxidation) and of scalable device fabrication, which constitute major difficulties not only in an application scenario but also for basic research. Substitutional doping, [START_REF] Wang | Doping Monolayer Graphene with Single Atom Substitutions[END_REF][START_REF] Robertson | Dynamics of Single Fe Atoms in Graphene Vacancies[END_REF] on the other hand, offers high stability thanks to the covalent bonding between the dopant atom and the surrounding carbon atoms in graphene.

Substitutional doping with non-magnetic elements such as N and B can be achieved with relatively simple processes, even during growth. [START_REF] Zhao | Visualizing Individual Nitrogen Dopants in Monolayer Graphene[END_REF][START_REF] Terrones | The Role of Defects and Doping in 2D Graphene Sheets and 1D Nanoribbons[END_REF] Substitutional doping with typical impurities such as Si (present during growth or processing) has also been observed. [START_REF] Zhou | Direct Determination of the Chemical Bonding of Individual Impurities in Graphene[END_REF][START_REF] Ramasse | Probing the Bonding and Electronic Structure of Single Atom Dopants in Graphene with Electron Energy Loss Spectroscopy[END_REF][START_REF] Susi | Manipulating Low-Dimensional Materials down to the Level of Single Atoms with Electron Irradiation[END_REF] However, substitutional incorporation of magnetic elements such as transition metals is generally not thermodynamically favored due to their different chemical nature compared to C. Substitutional doping with transition elements [START_REF] Wang | Doping Monolayer Graphene with Single Atom Substitutions[END_REF][START_REF] Robertson | Dynamics of Single Fe Atoms in Graphene Vacancies[END_REF][START_REF] Dyck | Doping Transition-Metal Atoms in Graphene for Carbon[END_REF] has been achieved using non-equilibrium conditions: first, vacancies are formed by bombardment with electrons or ions, and subsequently the magnetic dopant atoms are deposited, filling the vacancies that have not dynamically annealed between the two steps. Such dual-step processes are difficult to control, especially given the complex dynamics of vacancies in graphene, thus strongly limiting their reliability, reproducibility and scalability. Ultra-low energy (ULE) ion implantation has the potential to overcome these limitations and has been successfully applied to the substitutional doping of graphene with N and B, [START_REF] Bangert | Ion Implantation of Graphene -Toward IC Compatible Technologies[END_REF][START_REF] Willke | Short-Range Ordering of Ion-Implanted Nitrogen Atoms in SiC-Graphene[END_REF][START_REF] Kepaptsoglou | Electronic Structure Modification of Ion Implanted Graphene: The Spectroscopic Signatures of p-and n-Type Doping[END_REF][START_REF] Willke | Doping of Graphene by Low-Energy Ion Beam Implantation: Structural, Electronic, and Transport Properties[END_REF][START_REF] Cress | Nitrogen-Doped Graphene and Twisted Bilayer Graphene via Hyperthermal Ion Implantation with Depth Control[END_REF] P [START_REF] Susi | Single-Atom Spectroscopy of Phosphorus Dopants Implanted into Graphene[END_REF] and Ge. [START_REF] Tripathi | Implanting Germanium into Graphene[END_REF] Here, we report the incorporation of substitutional Mn atoms in graphene using ULE ion implantation. We selected Mn as a model case of magnetic doping, since it is expected to exhibit the largest magnetic moment among 3d transition metals in a single C vacancy, [START_REF] Krasheninnikov | Embedding Transition-Metal Atoms in Graphene: Structure, Bonding, and Magnetism[END_REF] which is typically advantageous for magnetic functionality. We provide a detailed characterization of the atomic structure of substitutional Mn in graphene, based on scanning tunneling microscopy (STM) and density functional theory (DFT) calculations. Furthermore, using angle-resolved photoemission spectroscopy (ARPES), we show that graphene doped with substitutional Mn retains the Dirac-like behavior of pristine graphene. These findings mo-tivate further research on magnetic functionalization of 2D materials, strongly benefiting from the high stability of substitutional doping, together with the high degree of control, flexibility and scalability of ultra-low energy ion implantation. Sample preparation and basic characterization. The samples consist of epitaxial graphene grown by chemical vapor deposition (CVD) on epitaxial Cu(111) thin films grown on sapphire(0001) substrates as described in Ref. 27. ULE ion implantation was performed by electrostatic deceleration of a 55 Mn + ion beam from 30 keV to 40 eV, with perpendicular incidence with respect to the sample surface, to a fluence of 1.5 × 10 14 ions per cm 2 (fluence rate of 1.7 × 10 11 cm -2 • s -1 ) measured by integration of the electric current on sample during implantation. [START_REF] Bangert | Ion Implantation of Graphene -Toward IC Compatible Technologies[END_REF][START_REF] Kepaptsoglou | Electronic Structure Modification of Ion Implanted Graphene: The Spectroscopic Signatures of p-and n-Type Doping[END_REF][START_REF] Willke | Doping of Graphene by Low-Energy Ion Beam Implantation: Structural, Electronic, and Transport Properties[END_REF] The samples were subsequently studied in the as-implanted state as well as after thermal annealing under ultra-high vacuum UHV (below 10 -9 mbar) up to 700 • C (20 min). This article focuses on the 700 • C annealed state, where STM measurements provide the clearest identification (compared to the other annealing stages) of Mn defects that, according to the analysis described in the following paragraph, we identify as substitutional Mn atoms in graphene. Moreover, this STM-based identification allows us to quantify its concentration as 0.04±0.02% with respect to C atoms in graphene (Figure 1 (a, b) and Table S1), which corresponds to approximately 1% of the implanted ions (i.e. 1% of the implanted 1.5 × 10 14 ions per cm 2 ) being substitutionally incorporated. The remainder corresponds to Mn ions in other configurations: backscattered during implantation (i.e. back into the vacuum upon collision with C and Cu atoms), adsorbed onto graphene, intercalated between graphene and the Cu surface, or diffused into Cu upon annealing. Nonsubstitutional configurations are also observed in the STM data, although not as clearly identified (Figure S3). Synchrotron X-ray photoelectron spectroscopy (XPS) measurements reveal a decrease in Mn concentration (with respect to C in the graphene layer) upon annealing: from 2.3% and 2.2% after thermal annealing at 200 • C and 425 • C, respectively, to 1.0% after annealing at 700 • C (Figure 1(c), additional details are given in supplementary information). We did not determine the Mn concentration in the as-implanted state, i.e. asmounted in the UHV chamber after exposure to ambient conditions, since in that state the surface is inherently highly contaminated, hindering quantification and possibly modifying the sample surface upon photon irradiation (more details in supplementary information).

Results and Discussion

The Mn concentration at low annealing temperature (2.3% for at 200 • C) corresponds to about half of the Mn fluence (1.5 × 10 14 cm -2 ), indicating that approximately half of the implanted ions are in fact backscattered (lost to the vacuum) already during the implantation process (but still deposit their charge, thus being accounted in the fluence value).

The further decrease upon annealing (down to 1% at 700 • C) can be understood as due to desorption of Mn atoms (i.e. loss of Mn to the vacuum) which had initially adsorbed onto graphene during implantation, and diffusion of intercalated Mn deeper into the Cu layer (i.e.

the Mn atoms remain in the sample but their apparent intensity decreases due to inelastic scattering of the photoelectrons). Thermal annealing also removes adsorbed contaminants resulting from exposure to air between the multiple synthesis and characterization steps (cf.

Table S3 and Table S4). Raman spectroscopy (Figure 1(d-f)) shows onset of significant disorder upon implantation (e), with the emergence of a D peak with a peak-height ratio with respect to the G peak (I D /I G ) of 2.14, whereas in the pristine sample (d) no D peak is observed. A low-disorder state is recovered upon annealing at 700 • C (I D /I G = 0.19). The residual D peak intensity is likely related to the substitutional Mn defects, although other defects (formed upon implantation or annealing) cannot be excluded. In short, the (partial) removal of non-substitutional Mn components, combined with a cleaner graphene surface and low disorder, make the 700 • C annealed state particularly suitable for a more detailed study of substitutionally doped graphene, which constitutes the remainder of this article.

Atomic-scale structure of substitutional Mn. By combining STM measurements and DFT simulations, we are able to analyse the atomic-scale structure of substitutional Mn in extreme detail. In the following, we present the main findings regarding the identification and structure of the defects. Further details are given in supplementary information. We must first consider the stacking between the graphene lattice and the Cu(111) surface. As a result of the different lattice parameter (2.46 Å for graphene and 2.56 Å for the Cu surface) and of the angular orientation of graphene with respect to the Cu surface (±3.57 • , determined using low-energy electron diffraction, cf. TOP-HCP (c), HCP-FCC (d). Between these points the stacking has lower symmetry. These high-symmetry points define a periodic 2-dimensional lattice, which is observed in the STM topographies as a Moiré superstructure (Figure 2(e)). The identification of the three highsymmetry stacking regions in the STM topographies (Figure 2(e)) is based on the simulated Cu-graphene distances (2.60 Å for TOP-FCC, 2.70 Å for TOP-HCP, and 2.87 Å for HCP-FCC), which is reflected in the STM data as height differences, i.e. lowest for TOP-FCC and highest for HCP-FCC (Figure S7). For a given high-symmetry stacking, and assuming that only one C vacancy is available, Mn can substitute for one of two non-equivalent C atoms.

For example, for TOP-FCC stacking, Mn can substitute for the C atom that is vertically aligned either with the TOP site or the FCC site of the Cu surface (Figure 3). We performed DFT calculations for all six combinations of three stackings and two sites. Furthermore, by comparing DFT-simulated and experimental STM topographies, we identify the surface structure of Mn substituting for a single C atom, i.e. located at a single C vacancy (Figure 3). This unambiguous identification is based not only on the comparison to simulations of substitutional Mn in a single vacancy, but also to simulations of a wide range of other defects, intrinsic or Mn-related: free single vacancies, free double vacancies, Mn in double vacancies, intercalated Mn (without vacancy), all of which can be distinguished from substitutional Mn in a single vacancy (Figure S1). Some of these other defects are indeed observed in our samples (Figure S3), although their identification is more ambiguous than that of substitutional Mn, which is the focus of the present article. We also note that we observe substitutional Mn defects already in the as-implanted state; however, as mentioned above, given the cleaner surface and lower degree of disorder, here we focus on the 700 • C annealed state. For all the six analyzed substitutional configurations, we find that the STM simulations reproduce the shape and contrast observed in the STM measurements. The TOP-FCC case is illustrated in Figure 3; additional simulations are shown in Figure S2. By referring to our structural model (Figure 3(d-g)) we are able to identify the six protrusions producing the characteristic triangular symmetry as the carbon atoms from the sublattice opposite to the vacancy sublattice (corresponding to the substitutional Mn site), with the three nearest neighbours to the substitutional Mn exhibiting the highest intensity.

In contrast, the Mn site appears as a dark region, consistent with the lower position of the Mn atom with respect to the surrounding C atoms. In Figure 3(a), we exhibit in particular two substitutional Mn defects in opposite graphene sublattices, showing triangular features which are mirror reflected with respect to each other, again in agreement with the symmetry relation obtained in the DFT simulations (Figure 3(b,c)). The site identification of the Mn defects in Figure 2(e) is based on this relation between stacking region, Mn site (i.e. C sublattice), and orientation of the triangular-shaped objects in the STM data; no other differences were observed in the topographic signatures for substitutional Mn in different stacking regions. showing the Moiré superstructure and substitutional Mn defects within them: TOP-FCC stacking with Mn on the fcc site (marked with †); TOP-HCP stacking with Mn on the HCP site (marked with ‡). The identification of the stacking region and of Mn site is described in the main text. Symbols ( , , ) are overlaid to link the identification of the three stacking regions in the various panels.

A particularly important aspect of the atomic structure of substitutional Mn is the vertical displacement with respect to the graphene layer (Figure 3(e,g)). Due to the short bond-lengths that characterize graphene, and to the large effective radius of Mn compared to the free space left by the vacancy, substitutional Mn atoms relax away from the graphene plane (Figure 3(e,g)). This out-of-plane displacement has been proposed to occur for various substitutional dopants in graphene (e.g. P [START_REF] Susi | Single-Atom Spectroscopy of Phosphorus Dopants Implanted into Graphene[END_REF] ). However, previous studies on such effects have focused on suspended graphene, which has mirror symmetry with respect to the graphene plane, before being broken by the out-of-plane displacements upon substitution.

In the present work, the Cu substrate breaks that mirror symmetry, and we observe that Mn is displaced towards the Cu surface, while the carbon atoms surrounding the Mn atom are slightly buckled away from the Cu surface, with an overall corrugation not exceeding 0.1 Å. The Mn-Cu distance is ∼ 2 Å while the equilibrium distance between the surface and the graphene layer is ∼ 3 Å. Interestingly, the C nearest neighbours (to Mn) do not retain, upon relaxation, the reconstruction that is typical of the single vacancy (taken as starting configuration); instead they symmetrically displace in plane in order to accommodate the substitutional Mn atom (the C-C bond distance within the same sublattice is 2.82 Å, corresponding to an increase of 11% compared to the same distance obtained far from the Mn defect). This configuration is obtained as final (relaxed state) when starting with the Mn atom just below the graphene layer. The configuration with the Mn atom above the vacancy site (protruding outwards, into the vacuum) can also be obtained as a final state of the DFT simulation, by placing the Mn atom exactly at the vacancy position and allowing for relaxation. For this final structural configuration (protruding outwards), the STM simulation shows a very distinct signature: bright Mn position versus dark surrounding C atoms (Figure S1(c)), i.e. the opposite of what we observe for the defects that we identify as substitutional Mn. From STM data, which is local in nature, we cannot completely exclude that the outwards-protruding configuration is also present for a minority fraction of the Mn atoms (cf. one rare case in Figure S3(a) that resembles the DFT simulation in Figure S1(c)).

Nevertheless, this outwards configuration is very rare, we only detected a few objects in all our STM data that may be associated with it, i.e. more than two orders of magnitude lower probability compared to the inwards configuration. This absence (or unlikely occurrence) of the outwards configuration is consistent with our DFT calculations, yielding a higher total energy (by 2.25 eV) compared to the inwards configuration (Figure 3(e, g)).

Table 1: Total energy difference ∆E tot for substitutional Mn in different stacking-site configurations, with respect to TOP-FCC stacking with Mn in FCC site. For each configuration, the magnetic moment of Mn (µ Mn ) and of the C atoms (µ C ) surrounding the vacancy (average) is also given.

Stacking

Site 1). For the most stable TOP-FCC stacking the two possible sites, namely FCC and TOP, differ for 0.33 eV. A similar value is also found for the two substitutional configurations of the TOP-HCP stacking while the site dependence is less marked in the HCP-FCC case. On the other hand taking the same site for different stacking we found energy differences ranging between 0.12 (TOP-site) and 0.64 eV (FCC-site). On the basis of these results we can conclude that the substitution in TOP position is unfavoured while the FCC site is preferred and its stability is dependent on the local stacking. We were able to clearly identify the stacking-site configuration (Figure 2(e)) for only a few substitutional defects, since it requires a clear imaging of the defect as well as of the Moiré superstructure (e.g. Figure 2(e)). Although the observed trend is consistent with the calculated total energies, i.e. predominance of the two lowest energy configurations (TOP-FCC stacking with Mn in the FCC site and TOP-HCP stacking, with Mn in the HCP site), the present statistics is too low to be conclusive. Future studies should investigate these effects in more detail. Nevertheless, our DFT calculations yield similar electronic properties for all the six simulated stacking-site configurations of substitutional Mn. The Mn atom is characterized by a reduction of electrons (-0.33 e) with respect to the ideal atomic valence. The n-doping of the graphene layer (+0.62 e) is slightly smaller than in the bare vacancy case (+0.75 e). The electrons donated by the Mn atom provide about half of the charge that is necessary for charge stability of the graphene+Mn overlayer, resulting in less charge donated by the substrate (-0.30 e) compared to the bare vacancy case (-0.75 e). The charge transfer also influences the magnetic properties of the system (Table 1), in particular in the close surrounding of the substitutional atom. We obtain calculated magnetic moments for Mn (µ Mn ) around 3.4 µ B , and for the C nearest neighbours (µ C ) around 0.2 µ B , i.e. lower in absolute value compared to the bare vacancy case (0.49 µ B ). Moreover, µ C is anti-parallel with respect to µ Mn (the reason for the negative values of µ C in Table 1).

While so far we discussed essentially local electronic properties, we now turn our attention to the effect on Mn doping on the overall electronic properties (band structure) of graphene.

Our DFT calculations yield a magnetic ground state for substitutional Mn in graphene/Cu, as can be deduced from the density of states (Figure 4(a,b)) projected on different species (PDOS). In addition to the spin-polarized PDOS of Mn we find a 0.15 eV shift between majority and minority spin components at the Dirac cone of graphene (Figure 4(a, b) and Figure S10), confirming the interaction between the Mn atom and the carbon overlayer. In that the Cu-graphene coupling is not significantly affected; no change in Fermi velocity can be resolved (0.91(3)×10 6 ms -1 and 0.88(4)×10 6 ms -1 for pristine and implanted, respectively).

The observed n-doping effect in the pristine sample is well-understood as due to charge transfer from the Cu surface; 29 the increase upon Mn implantation and annealing is unlikely to be a result of substitutional Mn doping, but instead an effect of the intercalated Mn component, for the following two reasons. First, our DFT calculations for substitutionally doped graphene show no increase of n-doping effect, which is instead found for intercalated Mn (Figure S9). Second, after 700 • C annealing, our STM data shows intercalated objects which are likely to correspond to intercalated Mn islands (Figure S3(e)). Based on the understanding that charge transfer effects in graphene-metal interface are mostly dependent on the work functions, [START_REF] Khomyakov | First-Principles Study of the Interaction and Charge Transfer between Graphene and Metals[END_REF] and considering that the work function of Mn [4.1(2) eV] is lower than that of Cu [4.65 (5) eV], [START_REF] Drummond | Work Functions of the transition Metals and Metal Silicides[END_REF] we can therefore expect a measurable n-doping effect from the Mn islands. Moreover, if intercalated Mn (as isolated Mn atoms) are still present after 700 • C (i.e. have not segregated into islands and have not diffused in the the Cu layer), they are also expected to contribute to n-doping on the basis of our DFT calculations. The band gap in the pristine sample is a well-known effect for graphene on Cu(111) [START_REF] Avila | Exploring Electronic Structure of One-Atom Thick Polycrystalline Graphene Films: A Nano Angle Resolved Photoemission Study[END_REF][START_REF] Vita | Understanding the Origin of Band Gap Formation in Graphene on Metals: Graphene on Cu/Ir (111)[END_REF] (resulting from the hybridization between Cu 3d and graphene π states), which is overestimated due to directional broadening when the ARPES beam probes multiple grains which are not perfectly aligned with respect to each other, [START_REF] Avila | Exploring Electronic Structure of One-Atom Thick Polycrystalline Graphene Films: A Nano Angle Resolved Photoemission Study[END_REF] as is the present case. The apparent decrease in band gap in the Mn-doped sample is likely due to Mn DOS (mostly d z 2 orbitals) emerging near the Dirac point (cf. arrows in Figure 4 (a, b, e, j), and hybridizing with the p z states of graphene, as predicted by our DFT calculations (Figure 4 (a, b)). Indeed, in the Mn-doped case, while an energy gap can still be extracted by taking the energy difference between local DOS maxima, (arrow in Figure 4(j)), the DOS is nearly constant in that region, i.e.

there is effectively no gap. In short, the Dirac-like character of graphene is preserved upon substitutional Mn doping, with strong indication of hybridization between Mn and graphene states near the Dirac point. This makes Mn-doped graphene an ideal system in which to study strong coupling between local magnetic moments and Dirac electrons, as we discuss below in the context of future research.

Substitutional Mn formation and stability. There are essentially two possible paths for the formation of substitutional Mn in graphene by ULE ion implantation: (1) direct substitution in which the impinging Mn ion displaces a C atom and takes its position; (2) a diffusion-assisted mechanism in which an implanted Mn atom (adatom or intercalated) diffuses (in plane) and becomes trapped in a C vacancy which had been created by an arbitrary Mn ion. Both mechanisms are likely to take place in the present case and both require the formation of a C vacancy, which by simple kinematics (conservation of kinetic energy and momentum) [START_REF] Alford | Fundamentals of Nanoscale Film Analysis[END_REF] requires a minimum initial energy of the impinging ion given by:

E min = E t (M Mn + M C ) 2 4M Mn M C (1)
where E t is the C threshold displacement energy (i.e. minimum energy that must be transferred to the C atom to create a stable C vacancy), and M Mn and M C are the atomic mass of Mn and C, respectively. With E t ≈ 22 eV, 34 equation 1 gives E min ≈ 37 eV. This minimum energy corresponds to a head-on collision (zero impact parameter). The larger the impact parameter, the smaller the transferred energy becomes, i.e. Mn atoms impinging with larger impact parameter with respect to the closest C atom must have an energy higher than 37

eV in order to still displace a C atom. We chose an implantation energy of 40 eV, i.e. above the threshold of 37 eV, so that a non-vanishing fraction of the implanted Mn ions transfer sufficient energy to create a vacancy, but not too high, to minimize the Mn backscattering energy (leading to loss of Mn to the vacuum).

From our DFT calculation we can also evaluate the stability of substitutional Mn. The binding energy between an implanted Mn atom and a vacancy can be estimated as:

E b = (E Mn-sub -E gr@Cu ) -(E Vac-gr@Cu -E gr@Cu ) -(E Mn-interc -E gr@Cu ) = -5.3 eV (2)
where E Mn-sub and E Mn-interc are the total energies of the system with Mn in substitutional and intercalated positions (i.e. with and without an associated carbon vacancy), E gr@Cu is the total energy of pristine graphene on Cu(111), and E Vac-gr@Cu of that system with a single vacancy. Such a strong binding (-5.3 eV) confirms the tendency of the implanted Mn atoms to combine with existing vacancies in the graphene lattice, as well as the high thermal stability of the substitutional configuration (at least up to 700 • C, as reported here).

While in the present case (ULE ion implantation), both direct and diffusion-assisted substitution is likely to occur, two-step doping approaches in which the vacancy formation and the incorporation of transition metals (e.g. Fe, Co, Ni, Mo, Pt) correspond to two separate processes, [START_REF] Wang | Doping Monolayer Graphene with Single Atom Substitutions[END_REF][START_REF] Robertson | Dynamics of Single Fe Atoms in Graphene Vacancies[END_REF][START_REF] Dyck | Doping Transition-Metal Atoms in Graphene for Carbon[END_REF][START_REF] Rodriguez-Manzo | Trapping of Metal Atoms in Vacancies of Carbon Nanotubes and Graphene[END_REF] are essentially based on the diffusion-assisted mechanism. For example, vacancies are formed by electron irradiation, followed by thermal deposition of the metal atoms, which diffuse on the surface until becoming trapped in the previously formed vacancies. Diffusion-assisted substitution has been observed to be thermally activated (i.e.

occurring only upon annealing above room temperature); [START_REF] Rodriguez-Manzo | Trapping of Metal Atoms in Vacancies of Carbon Nanotubes and Graphene[END_REF] this would imply that, if part of substitutional Mn in our samples was indeed formed via a diffusion assisted-mechanism (i.e.

not exclusively by direct substitution), the post-implantation annealing played a critical role.

Since its atomic number is similar to Mn (Z = 25), it is interesting to further compare our results (in terms of defect structure, electronic properties and stability) to those previously obtained for Fe (Z = 26) doping via two-step processes. The atomic structures of the two types of defects are similar (Mn here and Fe in Ref. 12, in a single C vacancy), characterized by strong relaxation effects resulting from the large effective radius of Mn/Fe compared to the C atom that they substitute: in-plane displacement of the neighboring C atoms away from the metal atom; out-of-plane displacement of the metal atom and (although less pronounced) of the neighboring C atoms. There is, however, a difference in the latter (outof-plane displacement), already mentioned above, resulting from the fact that in the present work the graphene layer is supported on a substrate while in Ref. 12 it is suspended. While in suspended graphene, it is generally accepted that the out-of-plane displacement of the metal dopant and of neighboring C atoms occurs in the same direction, in our case, as mentioned above, the Cu substrate breaks the mirror symmetry with respect to the graphene plane.

Consequently, the metal dopant (Mn) is observed to displace towards the Cu surface whereas the neighboring C atoms are displaced away from it, resulting in a total energy that is 2.25 eV lower than for a configuration in which both the Mn and C atoms are displaced away from the Cu surface. In terms of electronic structure there is a striking difference between substitutional Mn and Fe: while for Mn the electrons in (different) non-bonding d orbitals add up to a local moment of approximately 3 µ B , in Fe the 3d spins pair up resulting in a vanishing local moment. [START_REF] Krasheninnikov | Embedding Transition-Metal Atoms in Graphene: Structure, Bonding, and Magnetism[END_REF] In terms of stability, the fact that we observe no signs of instability up to 700 • C supports the calculated high binding energies of metal atoms such as Mn and Fe in single C vacancies: ∼ -5 eV for substitutional Mn in graphene/Cu(111) as reported here, and ∼ -6 eV and ∼ -7 eV for Mn and Fe, respectively, in suspended graphene. [START_REF] Krasheninnikov | Embedding Transition-Metal Atoms in Graphene: Structure, Bonding, and Magnetism[END_REF] The fact that we do not observe signs of clustering of substitutional Mn atoms up to 700 • C is also consistent with calculated migration barriers of ∼ 3 eV for 3d metals. [START_REF] Krasheninnikov | Embedding Transition-Metal Atoms in Graphene: Structure, Bonding, and Magnetism[END_REF] As mentioned in Ref. 35, the lower stability observed in that work is likely a consequence of a different defect structure, i.e. not an ideal substitution into a single C vacancy.

Future studies should evaluate the roles of direct and diffusion-assisted substitution as well as the effect of thermal annealing and varying ion energy on the substitution efficiency, in particular, towards higher substitutional concentrations. In the present work, with an implantation energy of 40 eV, the 0.04% substitutional concentration (determined from STM) corresponds to about 1% efficiency with respect to the implanted fluence. Further optimization, in particular towards slightly higher energies, is likely to further increase the substitutional efficiency, by increasing the vacancy formation rate, without significantly increasing the Mn backscattering energy. In any case, the substitutional concentration can in principle be increased by increasing the fluence. For example, a substitutional concentration of about 1% (with a substitution efficiency of a few %) was obtained with ULE ion implantation (25 eV) of Nitrogen 20 (for N, equation 1 gives E min ≈ 22 eV, i.e. also slightly below the implantation energy, as in the present case).

Prospects for magnetic functionality. As mentioned above, ARPES measurements at 80 K do not provide evidence of the spin-splitting (ferromagnetic state) predicted by DFT (corresponding to an effective temperature of 0 K). Various reasons may explain this. First, the ordering temperature may simply be below 80 K. Second, the local magnetic moment of substitutional Mn may in fact be strongly decreased or completely screened by the conduction electrons (Kondo effect, which has been reported for magnetic defects in graphene at temperatures as high as 90 K 36 ), which would inhibit a ferromagnetic state. Future studies should therefore specifically address the magnetic properties of substitutionally doped graphene.

On the theoretical side, more complex approaches including many-body effects related to electronic correlation can be adopted, such as tight binding model with Hubbard term [START_REF] Liang | Modeling the Kondo Effect of a Magnetic Atom Adsorbed on Graphene[END_REF] or DFT+dynamical mean-field theory (DMFT) in the Anderson Impurity Model 38 which allow to describe Kondo physics. On the experimental side, various advanced techniques can be applied that provide insight at various scales, for example: local studies based on scanning tunneling microscopy and spectroscopy, including with spin-polarized capabilities; [START_REF] Ren | Effect of Cobalt Adatoms on a Graphene Monolayer Controlled by Substrate-Induced Ripples[END_REF][START_REF] Ren | Interatomic Spin Coupling in Manganese Clusters Registered on Graphene[END_REF][START_REF] Choi | Studies of Single Atom Magnets via Scanning Tunneling Microscopy[END_REF] macroscopic studies based on X-ray magnetic circular dichroism (XMCD); [START_REF] Girovsky | Long-Range Ferrimagnetic Order in a Two-Dimensional Supramolecular Kondo Lattice[END_REF][START_REF] Joly | Screening of the Spin and Orbital Magnetic Moments of Fe Impurities in Cu[END_REF] spin-resolved band structure characterization based on spin-resolved ARPES; 42 magnetotransport studies (e.g. Hall effect and magnetoresistance measurements). [START_REF] Chen | Tunable Kondo Effect in Graphene with Defects[END_REF] Magnetic doping of 2D materials by ULE ion implantation. Our approach to magnetic substitutional doping, based on ULE ion implantation, offers various advantages.

For example, contrary to two-step approaches, [START_REF] Wang | Doping Monolayer Graphene with Single Atom Substitutions[END_REF][START_REF] Robertson | Dynamics of Single Fe Atoms in Graphene Vacancies[END_REF][START_REF] Dyck | Doping Transition-Metal Atoms in Graphene for Carbon[END_REF] we can control both processes of vacancy creation and of dopant incorporation with a single set of implantation parameters, which can be expected to yield better reproducibility. All the implantation parameters (energy, fluence, fluence rate, and temperature) can be precisely monitored and controlled, which in general provide high degree of control and reliability. Moreover, the resulting atomic configuration, in which the Mn atom is underneath (protected by) the graphene layer, can be expected to provide more stability with respect to exposure to the atmosphere and post-processing steps (e.g. deposition of a polymer layer for transfer or of an oxide layer for top gating).

The approach is also extremely flexible, being compatible with virtually any element of the periodic table (as dopant) and any 2D material (as host). Although one can expect some limitations regarding which elements can effectively be substitutionally incorporated into a given 2D material, those should be less restrictive compared to processes that occur closer to equilibrium (e.g. doping during growth). Therefore, the present work motivates further research, even beyond magnetic doping of graphene, for example: magnetic doping of MoS 2 and related 2D semiconductors for spintronics and valleytronics applications; 43 creation of spin defects in h-BN and related 2D insulators for quantum technologies. [START_REF] Tran | Quantum Emission from Hexagonal Boron Nitride Monolayers[END_REF][START_REF] Castelletto | Hexagonal Boron Nitride: A Review of the Emerging Material Platform for Single-Photon Sources and the Spin-Photon Interface[END_REF] Finally, ion implantation is widely used for semiconductor processing in microelectronics; if ULE ion implantation technology is further developed to the level of its higher energy counterpart, one can expect it to be compatible with wafer-scale fabrication.

Conclusion

To conclude, we demonstrate the incorporation of substitutional Mn in graphene using ULE ion implantation and provide a detailed characterization of its local atomic structure. We also show that graphene doped with substitutional Mn retains the Dirac-like behavior of pristine graphene, making it an ideal system in which to study the electronic and magnetic interplay between local magnetic moments and Dirac electrons. These results establish a solid basis for such future studies, for example, dealing with ferromagnetism and Kondo physics in magnetically doped 2D materials. Although this work deals with Mn doping of graphene, our approach is compatible with virtually any dopant element and any 2D material host. Scanning tunneling microscopy. The STM micrographs reported here were obtained in UHV (base pressure ∼ 10 -11 mbar) at low temperatures (2 K, 5 K and 78 K) using an

Methods

Omicron LT STM and an Unisoku USM-1500. We used electrochemically etched W tips and cut Pt-Ir tips. The tip oxide from the W tips is first removed by flash annealing. Both types of tips were characterized by scanning the Au(111) surface. All topographies were acquired in constant I tun. . X-ray photoelectron spectroscopy. XPS measurements were performed at the Su-perESCA beamline at Elettra synchrotron in Trieste. All the spectra were acquired in UHV (base pressure ∼ 10 -10 mbar) and at room temperature (RT), the samples were placed in the normal emission (NE) configuration, with respect to the hemispherical analyzer. The samples were annealed in-situ up to 700 • C, with a ramp rate of 1 • C/sec to monitor the change of Mn concentration. The C 1s, O 1s, and Mn 2p core-level spectra were acquired at photon energies of 400 eV, 650 eV, and 780 eV, respectively. The core-level spectra were fitted with Doniach-Šunjić functions convoluted with a Gaussian distribution and a linear background (apart from Mn 2p core-level, where a Shirley background is used).

Raman spectroscopy. Raman spectra were measured using a confocal Raman microscope (Monovista CRS+, S&I GmbH) equipped with a 532 nm Nd:YAG laser. The laser was directed onto the sample surface through an objective (OLYMPUS, BX43 100×, N.A. 0.7), with the maximum laser power remaining below 1 mW in order to avoid laser-induced modification. All the measurements were obtained in ambient conditions, at RT.

Angle-resolved photoemission spectroscopy. ARPES measurements were performed at the BaDElPh beamline at Elettra synchrotron in Trieste, [START_REF] Petaccia | ElPh: A 4 m Normal-Incidence Monochromator Beamline at Elettra[END_REF] with a photon energy of 34 eV. Before the measurements, the sample was pre-annealed at 700 • C. The ARPES spectra were obtained in UHV (base pressure ∼ 5 × 10 -11 mbar), at 80 K. The presented spectra correspond to a summation of data that were separately acquired using linear p-and spolarized light. The energy and angular resolution were set to 10 meV and 0.1 • , respectively.

Low-energy electron diffraction. As-grown graphene samples were characterized using low-energy electron diffraction (LEED) measured in an UHV chamber with background pressure better than 4 × 10 -9 mbar using a rear-view AEP 8011, VG Scientific instrument.

LEED images were recorded using the screen voltage and an electron beam (spot size ca.

1mm) energy of 1 kV and 52 eV, respectively.

Density functional theory calculations. We performed Density Functional Theory calculations with generalized gradient approximation (GGA) for the exchange-correlation potential, using a localized orbital basis set and a pseudopotential description of the core elec-trons, as implemented in the SIESTA code. In the adopted model graphene is matched to the underlying Cu(111) surface, by taking the theoretical lattice constant (a T heo = 3.59 Å), resulting in a 3% expansion of the overlayer. We considered one defect related to Mnimplantation in a (7 × 7) supercell, corresponding to having a Mn-C ratio of 1/97 ∼ 1%.

Instead of dealing with Moiré reconstruction that would require a very large periodicity, we separately analyze different the possible stacking configurations, corresponding to different high-symmetry regions of the Moiré, placing the Mn atom in different positions with respect to the graphene lattice and the underlying substrate. This choice is supported by the local nature of the graphene-Mn-Cu interaction in proximity of the defect, leading to electronic, magnetic and structural properties that are not dependent on the long range arrangement typical of the Moiré. We performed structural relaxation of the overlayer with the Mn atom while the Cu is taken as fixed to the relaxed coordinate of the surface. The threshold for the force tolerance was set to 0.04 eV/ Å. We adopted a mesh cutoff of 450 Ry and a 5 × 5 × 1 sampling of the Brillouin zone. Along the direction perpendicular to the surface three layers of Cu were considered for structural relaxation while a six layer Cu slab was adopted for the evaluation of the electronic properties. The periodic replica of the system along the same direction were separated by a 40 Å wide vacuum region.

The STM simulations were performed in a Tersoff-Hamman approach, 47 assuming a constant density of states for the tip. The integration range for the electronic density, corresponding to the applied bias, was varied in an energy interval of 1.0 eV around the Fermi

Quantification of surface Mn, C and O

To evaluate the effect of thermal annealing on the amount of Mn and contaminants present at the surface, XPS quantification analysis was performed on three core-levels: Mn 2p 3/2 , C 1s, O 1s (spectra presented in Figure 1(c) and Figure S5). In order to estimate the relative surface concentrations, the following three parameters are typically considered: backgroundsubtracted peak area (A), photoionization cross-section (σ) (Table S2), and photoelectron escape depth. Given the similar photoelectron kinetic energies (E k ) under our measurement conditions, we assumed the escape depth to be the same for the photoelectrons from the core-levels considered here; the error induced by this approximation is further decreased by the fact that all species are distributed within a very narrow depth range. As reference, we took the C signal in the pristine graphene sample, after annealing 700 • C, as one monolayer equivalent (MLE) of carbon in graphene. All other concentrations were taken as relative to this C monolayer. For example, the surface concentration of Mn (c M n ), with respect to one C MLE, was estimated as:

c M n = A M n /σ M n A C,700 /σ C,700
Table S3 compiles the concentration values obtained for Mn, C and O, for different annealing temperatures. We avoided measurements before thermal annealing, i.e. as-mounted in the UHV chamber after exposure to ambient conditions, since in that state the surface is inherently highly contaminated. Such adsorbed contaminants generally hinder the quantification of elements which are present underneath. Moreover, irradiation with high-flux X-rays may irreversibly modify the surface when such high concentrations of contaminants are present. The most striking feature is the decrease in Mn concentration upon 700 • C, as discussed in the main text (Figure 1(c)). Moreover, the significant decrease in surface O in the pristine sample, from 425 • C to 700 • C, confirms that surface contamination is significantly decreased. In the Mn implanted sample, the surface O content remains around 20%, even after annealing at 700 • C, indicating that surface oxidation may be enhanced by the Mn implantation and subsequent exposure to ambient, before the characterization steps.

Most importantly, the C signal of the Mn-implanted sample (c), in particular after annealing at 700 • C, preserves the signature of pristine graphene (a), consistent with the low level of disorder observed with Raman spectroscopy and STM. In short, and as also discussed in the main text, the (partial) removal of non-substitutional Mn components, combined with a cleaner graphene surface and low disorder, make the 700 • C annealed state particularly suitable for a more detailed study of Mn substitutionally doped graphene, which constitutes the focus of this article.

Graphene domains

The LEED image (Figure S6(a)) acquired at 52 eV shows the hexagonal principle diffraction spots of Cu(111) substrate (innermost spot from the set of three spots, see inset) surrounded (on both sides) by the diffraction spots of graphene (outermost spots), forming slightly larger hexagons owing to a smaller real space lattice parameter of graphene (a Gr = 2.46 Å) as compared to the surface lattice parameter of Cu(111) (a Cu = 2.56 Å). By applying trigonometric functions to the reciprocal lattice parameters (indicated on the image), we can identify that the graphene lattice is not perfectly aligned with that of Cu(111). In particular, we observe the presence of two rotated sets of graphene domains (rotated by ca. 5.62 ± 0.12 • with respect to each other), giving rise to two spots, i.e. rotated by 2.81 ± 0.06 • with respect to Cu in opposite directions. The presence of these two sets of domains is also observed in the

Identification of high-symmetry stacking regions

As mentioned in the main text, DFT calculations yield Cu-graphene distances of 2.60 Å for TOP-FCC, 2.70 Å for TOP-HCP, and 2.87 Å for HCP-FCC. We illustrate in Figure S7 the identification of the stacking regions within the Moiré superstructure, by extracting height profiles from the STM micrograph. On both profiles, height differences are observed which can be attributed to the three different stacking regions by comparison to the trend obtained in the DFT simulated heights, i.e. lower to higher from TOP-FCC, to TOP-HCP, to HCP-FCC.

Additional results from the DFT calculations

In this section, we provide additional results concerning the DOS and band structure determined from the DFT calculations. We note that the structural model used for these calculations focuses on the local structural properties of the different stacking regions, aiming in particular at the identification of the defects observed with STM. In particular, the long range Moiré pattern observed for graphene/Cu(111) is not accounted for in the calculations, since it would require much more computation-intensive calculations, with no additional insight of relevance for the present work. tional Mn (inwards configuration). The Dirac cone energy in the substitutional case (-0.39 eV and -0.54 eV from the Fermi level for spin up and spin down components, respectively) is slightly higher than the pristine case (-0.62 eV). This result further confirms that the additional n-doping observed experimentally is not due to substitutional Mn. A small gap at the Dirac cone is observed which is related to the carbon vacancy. The midgap state characteristic of the single C vacancy in graphene disappears upon incorporation of a Mn atom at the vacancy site.

In both cases, intercalated and substitutional, Mn states are superimposed to the graphene bands; the interaction between graphene and Mn states originates discontinuities in the Dirac cone (Figure S9 and S10). In particular, in the substitutional case (Figure S10), these are clearly visible just below the Fermi level, i.e. the region probed by ARPES. However, they are not observed in our experimental ARPES data, most likely due to the fact that the DFT calculations imply a periodic Mn arrangement, while the actual distribution in the sample is random.

Figure 1 :

 1 Figure 1: STM micrographs of Mn implanted graphene on Cu(111) after 700 • C annealing: (a) large-scale micrograph showing substitutional Mn defects pointed by white arrows (V sample = 200 mV, I tun. = 300 pA); (b) small-scale micrograph showing one substitutional Mn defect (V sample = -100 mV, I tun. = 5 nA). (c) Mn 2p core-level spectrum of Mn implanted graphene/Cu(111), annealed at the following temperatures after implantation: 200 • C (black), 425 • C (blue) and 700 • C (pink). (d-f) Raman spectra (laser wavelength = 532 nm) on graphene/Cu(111): (d) pristine sample annealed at 700 • C; (e) Mn implanted, in the as-implanted state; (f) Mn implanted and annealed at 700 • C after implantation. The positions of the D, G and 2D bands of graphene are indicated in panel (e).

  Figure S6(b)), the graphene-Cu stacking varies continuously across the surface (Figure 2(a)) leading to a Moiré superstructure (e). Three points with high-symmetry stacking can be defined indicating the position of the two C atoms in the graphene unit cell with respect to the underlying substrate: TOP-FCC (Figure 2(b))

Figure 2 :

 2 Figure 2: (a) Moiré superstructure of graphene/Cu(111), resulting from the different graphene-Cu stacking regions: (b) TOP-FCC, (c) TOP-HCP, (d) HCP-FCC. (e) STM micrograph (V sample = 100 mV, I tun. = 100 pA)showing the Moiré superstructure and substitutional Mn defects within them: TOP-FCC stacking with Mn on the fcc site (marked with †); TOP-HCP stacking with Mn on the HCP site (marked with ‡). The identification of the stacking region and of Mn site is described in the main text. Symbols ( , , ) are overlaid to link the identification of the three stacking regions in the various panels.

Figure 3 :

 3 Figure 3: (a) STM micrograph (V sample = -100 mV, I tun. = 5 nA) exhibiting two substitutional Mn defects in opposite graphene sublattices, showing triangular features which are mirror reflected with respect to each other. DFT simulated STM for substitutional Mn in a TOP-FCC stacking region (energy window [-0.25 eV,E F ]): (b) for Mn in FCC site; (c) in TOP site. Relaxed model structures obtained from DFT: (d, e) top and side view, respectively, for Mn in FCC site; (f, g) top and side view, respectively, for Mn in TOP site.

Figure 4 :

 4 Figure 4: (a, b) Density of states (PDOS) projected on the different elements (a) as well as on the different orbital components (i.e. magnetic quantum number m) of Mn and p z states of graphene (b), for a substitutional Mn concentration of 1%. The majority and minority spin components are plotted as positive and negative PDOS, respectively. (c-l) ARPES spectra (hν=34 eV) measured after 700 • C annealing, on pristine (c-g) and Mn implanted (h-l) samples. Panels (c, h) show 2D spectra measured near the K-point along the Γ → K direction; panels (d, i) display the second derivative of the data in (c, h), highlighting the Dirac cone dispersion of the π and π * bands of graphene. (e, j) Energy distribution curves (EDC) extracted from the p-polarization branch of the data in (c, h), with a wavenumber (k) step of 0.01 Å-1 . (f, k) 2D spectra measured in the region where the graphene π-band and the Cu 3d-band intersect; the observed discontinuity in the graphene π-band results from hybridization with the Cu 3d-band. (g,l) Cu(111) surface state measured at the Γ point. In panels (e, j), the dots mark the local maxima. The arrows in panels (a, b, e, j) point to the gap region, indicating that Mn d z states emerge in j according to the PDOS in a, b. In the ARPES panels (c-l), binding energy is defined as BE = E F -E.

Figure 4 (

 4 Figure4(c-l) compares ARPES data for pristine (c-g) and Mn-implanted graphene (h-l), both after annealing at 700 • C. The first observation is that the data (measured at 80 K) do not provide clear indication of the spin-splitting predicted by DFT (corresponding to an effective temperature of 0 K), i.e. a single set of Dirac cones is observed, without indication of 0.15 eV spin-splitting predicted by DFT (FigureS10). The DOS near the Dirac point does increase as predicted for substitutional Mn (cf. arrows in Figure4, panels a, b, e, j), but finer details cannot be resolved. In fact, the ARPES measurements (Figure4) show only a subtle effect of Mn doping in the band structure. We observe an increase in the n-doped character of the implanted sample (Figure4 (c, d) and (h, i)): the Fermi level position with respect to the Dirac point changes from 0.33(1) eV in the pristine sample to 0.43(1) eV in the implanted sample. A decrease in the band gap is also observed, from 0.40(2) eV to 0.23(2) eV. Other key features are essentially unaffected: the Cu surface state (g and l) is clearly visible, indicating that no significant oxidation of the interface has occurred; 28 the hybridization between Cu 3d and graphene π states (f and k) is also clearly visible, indicating

  CVD growth of graphene. Pre-cut and polished 2 Czochralski grown sapphire wafers (Roditi International Corporation) diced along the c-plane are used to prepare Al 2 O 3 template wafers. The monocrystalline wafers have a surface misorientation smaller than 0.3 • . The wafers were cleaned prior to Cu sputtering in a 3:1 concentrated acid mixture of H 2 O 4 :H 3 PO 4 at 200 • C for 20 min followed by a UPW rinse for 3 min. Cleaned sapphire wafers are placed in a Nimbus 310 sputtering setup with a base pressure of 4 × 10 -6 mbar. Substrates are placed using a 200 mm Si pocket wafer and sputtering is performed at room temperature for 173 s (21 passes under target) at 6 × 10 -3 mbar Ar pressure. The Cu deposition rate is calibrated and corresponds to a rate of about 29 Å s -1 . The applied power was 3000 W and the throw distance is approximately 50 mm. The graphene growth process has been optimized taking into account oriented monolayer growth and low Raman D peak. The growth conditions are 5000 sccm Ar: 125 sccm H 2 : 0.3 sccm CH 4 at a pressure of 750 mbar close to the melting point of Cu. The growth time is 30 min, while ramp-up and down is performed in a mixture of Ar and H 2 .

  level. We simulated a constant-height (3 Å ) STM image applying a Gaussian spatial broadening of 1 Å to the electronic density to mimic the finite experimental resolution associated with the tip. from the EU Horizon 2020 Framework (MagDirac project, ID:796940). K.V. also acknowledges funding from the Agency for Innovation & Entrepreneurship (VLAIO). The authors acknowledge Elettra Sincrotrone Trieste for providing access to its synchrotron radiation facilities and, in particular, to the BaDElPh and SuperESCA beamlines, as well as the Red Espanola de Supercomputacion (RES) for the availability of high performance computing resources and support under the project FI-2020-2-0034. The authors thank P. Lacovig and S. Lizzit for the support during the XPS beam time at SuperESCA, Elettra. The authors also thank P. Walke for the support with the Raman spectroscopy measurements.

Figure S2 :

 S2 Figure S2: STM simulations for substitutional Mn in the three high-symmetry stacking configurations (TOP-FCC, TOP-HCP and HCP-FCC) and, for each stacking configuration, in each of the two possible graphene sublattices. All the simulations are performed with the same tip-surface distance (3 Å), with the same tip size (1 Å), and are plotted in the same contrast scale. The integration energy ranges are indicated (0.5 eV below and 0.5 above the Fermi level E F ).

  Figure S3(a), DFT simulation in Figure S1(c)). Moreover, four different intercalated forms could be discerned. Two single Mn intercalated configurations, in the HCP position (STM in Figure S3(b), DFT in Figure S1(d)) and in the FCC position (STM in Figure S3(c), DFT in Figure S1(e)) were observed, as well as a possible intercalated Mn cluster (Figure S3(d)). Our data further indicates that these three forms of intercalated Mn aggregate into intercalated Mn islands (Figure S3(e)) upon annealing at ≥600 • C.Carbon single vacancyAlthough rare, objects resembling the DFT simulation for the carbon single vacancy (Figure S1(b)) appear in some micrographs (FigureS4), after annealing at 700 • C, suggesting that a negligible concentration of vacancies may be formed upon high temperature annealing

Figure S3 :

 S3 Figure S3: STM micrographs showing objects which can be attributed to other forms of Mn incorporation in graphene/Cu(111), i.e. other than substitutional in the inwards configuration: (a) substitutional Mn in the outwards configuration (V sample = 200 mV, I tun. = 900 pA); (b) single intercalated Mn in HCP position (V sample = -400 mV, I tun. = 100 pA); (c) single intercalated Mn in the FCC position (V sample = -400 mV, I tun. = 100 pA); (d) intercalated Mn cluster (V sample = -400 mV, I tun. = 100 pA); (e) intercalated (Mn or Mn/Cu alloy) islands (V sample = -200 mV, I tun. = 80 pA). Micrographs (a) and (e) were measured after annealing at 600 • C. Micrographs (b-d) were measured in the as-implanted state.

Figure S4 :

 S4 Figure S4: STM micrograph showing a possible carbon single vacancy, measured after 700 • C annealing (V sample = 100 mV, I tun. = 1.5 nA)

Figure S5 :

 S5 Figure S5: XPS measurements of C 1s (a, c) and O 1s (b, d) core-levels on pristine (a, b) and Mn-implanted (c, d) graphene/Cu(111), measured in the as-implanted state (black) and after annealing at 700 • C (red).

Figure S6 :

 S6 Figure S6: (a) LEED pattern of pristine graphene/Cu(111), observed at an electron beam energy of 52 eV. Inset: enlarged LEED spots to illustrate the two graphene spots and the sandwiched copper spot. (b) ARPES Fermi surface map of pristine graphene/Cu(111) near the K point, obtained with hν = 34 eV, at a constant binding energy of 0.33 eV (i.e. Dirac point).

Figure S7 :

 S7 Figure S7: (a) STM micrograph showing the Moiré superstructure with two arrows displaying extracted profiles (V sample = 100 mV, I tun. = 100 pA), (b) extracted profile from HCP-FCC region to HCP-FCC region and (c) second extracted profile for statistics.

Figure S8 :

 S8 Figure S8: Momentum distribution curves (MDC) extracted from the p-polarization branch of the data in Figure 4(c, h), corresponding to pristine (a) and Mn-implanted (b) graphene/Cu(111), with a binding energy (BE) step of 0.04 eV. The dots mark the local maxima.

Figure

  FigureS9shows the spin-resolved band structure of graphene/Cu(111) with intercalated Mn (same configuration as in FigureS1(d)). The band structure shows a downshift of the Dirac cone with respect to the pristine graphene/Cu(111) case due to charge transfer. The Dirac cone in pristine graphene/Cu is found at -0.62 eV from the Fermi level (b) while in Mn intercalated graphene/Cu it is at -0.86 and -0.64 for spin up and spin down components, respectively (c, d). As mentioned in the main text, it is plausible that the additional n-doping effect observed in the implanted sample (compared to pristine) is due to a non-vanishing concentration of such intercalated Mn atoms.

Figure

  FigureS10shows the spin-resolved band structure of graphene/Cu(111) with substitu-

Figure S10 :

 S10 Figure S10: Spin-resolved projected density of states (a) and band structure of graphene/Cu(111) with substitutional Mn atoms (c, d), compared to the band structure of pristine graphene/Cu(111) (b). Green dashed lines in the band structures mark the position of the Dirac point. The Mn/C ratio is 1/97 and the Mn/Cu ratio 1/295. The band structure has been unfolded onto high symmetry paths of the 1 × 1 Brillouin zone. The position of Mn states, as deduced by the PDOS, is marked.

Table S1 :

 S1 Estimation of the substitutional Mn concentration.

	Micrograph Area (nm 2 ) # of defects Concentration
	1	400	1	0.0071%
	2	100	1	0.0284%
	3	400	5	0.0355%
	4	100	1	0.0284%
	5	100	2	0.0568%
	6	400	3	0.0213%
	7	100	2	0.0568%
	8	100	2	0.0568%
	9	100	4	0.1136%
	10	100	2	0.0568%
	11	400	8	0.0568%
	12	400	6	0.0426%
	13	400	7	0.0497%
	14	400	4	0.0284%
	15	400	8	0.0568%
	16	100	3	0.0852%
	17	400	5	0.0355%
	18	400	8	0.0568%
	19	400	8	0.0568%
	based on comparison with the DFT simulated STM topographies. The outwards configu-
	ration of substitutional Mn was rarely observed (STM micrograph in

Table S2 :

 S2 Photoionization cross-sections.

	C 1s	O 1s	Mn 2p 3/2
	0.4247 (hν=400eV) 0.3383 (hν=650eV) 0.6746 (hν=780eV)

Table S3 :

 S3 Temperature dependence of O and C concentration on pristine graphene/Cu(111).

	Element	425 • C annealed 700 • C annealed
	O (%MLE)	51.9	9.4
	C (%MLE)	77.6	100

Table S4 :

 S4 Temperature dependence of Mn, O and C concentration on Mn implanted graphene/Cu(111).

	Element	200 • C annealed 425 • C annealed 700 • C annealed
	Mn (%MLE)	2.3	2.2	1
	O (%MLE)	25.4	22.4	22.4
	C (%MLE)	81.5	86	84.6

Acknowledgement

This work was funded by KU Leuven and FWO Vlaanderen. M.N.N. acknowledges funding

Supporting Information Available

The Supporting Information is available free of charge on the ACS Publications website at ... It provides additional results from the DFT calculations and from the experimental studies.

Supplementary information for: Doping Graphene with

Substitutional Mn STM simulations (DFT) for various defects

As mentioned in the main text, the unambiguous identification of substitutional Mn is based not only on the agreement with DFT simulations for that configuration, but also on the lack of agreement with simulations (Figure S1) of other point defects that may also be present: free single vacancies, free double vacancies, Mn in double vacancies, substitutional Mn in the outwards-protruding configuration, and intercalated Mn (without vacancy).

STM simulations (DFT) for substitutional Mn

As mentioned in the main text, for all the six high-symmetry substitutional configurations, the STM simulations (Figure S2) show generally similar features, and are in good agreement with the shape and contrast observed in the STM measurements and identified as substitutional Mn in the inwards configuration.

Quantification of substitutional Mn

The weighted average (0.04%) and standard deviation (0.02%) of the concentration of substitutional Mn (with respect to C atoms) were determined from the estimated concentration in 19 micrographs from different surface areas on the same sample (Table S1). The weight was determined based on the area of each micrograph with respect to the total area of all the considered micrographs.

Other forms of Mn incorporation

Different types of defects were observed in the STM micrographs that can be attributed to other forms of Mn incorporation (other than substitutional Mn in the inwards configuration),

ARPES analysis: energy and momentum distribution curves

The ARPES maps (p-polarization) were analyzed by extracting energy distribution curves (EDC) at constant k || , with a step of 0.01 Å-1 (Figure 4(e,j)) and momentum distribution curves (MDC) at constant binding energy, with a step of 0.04 eV (Figure S8). The extracted line profiles were fitted by single Lorentzian functions in order to determine their peak positions. The data points corresponding to those peaks in both the EDC and MDC line-profiles were then fitted with linear functions to find the two intercepts (valence and conduction bands) with the K point (wavenumber=1.7 Å-1 ). The obtained gap value is defined as the difference between the two intercepts, resulting from the non-colinear valence and conduction bands. For both MDC and EDC analysis, we obtained the position of Dirac point (with respect to the Fermi level), defined as the midpoint in the gap. The values for the gap and Dirac point position quoted in the main text are based on the EDC analysis, as more commonly adopted. Here, in Table S5, we provide the results from the MDC analysis, for completeness. The Fermi velocity values quoted in the main text were obtained from the MDC analysis as:

1 ∂E ∂k where = 6.58 × 10 -16 eV • s, and ∂E ∂k is taken as the slope of the conduction band, determined from the linear fit of the MDC local maxima near the Fermi level.