Lei Pan

Jean Dumoncel

Arnaud Mazurier

Clément Zanolli
email: clement.zanolli@gmail.com

Hominin diversity in East Asia during the Middle Pleistocene: A premolar endostructural perspective

Following the recent studies of East Asian mid-Middle to early Late Pleistocene hominin material, a large spectrum of morphological diversity has been recognized and the coexistence of archaic ('Homo erectus-like') and derived ('modern-like') dental morphological patterns has been highlighted. In fact, for most of these Chinese fossils, generally categorized as 'archaic Homo sapiens' or 'post-H. erectus Homo', the taxonomic attribution is a matter of contention. With the help of µCT techniques and a deformation-based 3D geometric morphometric (GM) approach, we focused on the morphological variation in the enamel-dentine junction (EDJ) of 18 upper and lower premolars from Chinese Middle Pleistocene hominins. We then compared our results with a number of fossil and modern human groups, including Early Pleistocene H. erectus from Sangiran; late Early Pleistocene hominins from Tighenif, Algeria; classic Neanderthals; and modern humans. Our results highlight an evolutionary/chronological trend of crown base reduction, elevation of EDJ topography, and EDJ surface simplification in the hominin groups studied here. Moreover, this study brings insights to the taxonomy/phylogeny of six late Middle Pleistocene specimens whose evolutionary placement has been debated for decades. Among these specimens, Changyang premolars show features that can be aligned with the Asian H. erectus hypodigm, while Panxian Dadong and Tongzi premolars are more similar to Late Pleistocene Homo. Compared with early to mid-Middle Pleistocene hominins in East Asia, late Middle Pleistocene hominins evince an enlarged morphological variation. A persistence of archaic morphotypes and possible admixture among populations during the late Middle Pleistocene are discussed.

Introduction

Dental evidence suggests that Eurasian Early and Middle Pleistocene hominins, including those from Central Asia [START_REF] Martinón-Torres | Dental remains from Dmanisi (Republic of Georgia): morphological analysis and comparative study[END_REF], continental East Asia [START_REF] Bailey | A comparative dental metrical and morphological analysis of a Middle Pleistocene hominin maxilla from Chaoxian (Chaohu), China[END_REF][START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF][START_REF] Zanolli | Inner tooth morphology of Homo erectus from Zhoukoudian. New evidence from an old collection housed at Uppsala University, Sweden[END_REF] and insular Southeast Asia, had a relatively independent evolutionary course compared with African hominins [START_REF] Martinón-Torres | Dental evidence on the hominin dispersals during the Pleistocene[END_REF][START_REF] Bermúdez De Castro | A new model for the evolution of the human Pleistocene populations of Europe[END_REF][START_REF] Kaifu | Archaic hominin populations in Asia before the arrival of modern humans: their phylogeny and implications for the "Southern Denisovans[END_REF]. From the early to mid-Middle Pleistocene, hominin remains found in the Asian continent are mostly assigned to Homo erectus [START_REF] Kaifu | Taxonomic affinities and evolutionary history of the Early Pleistocene hominids of Java: Dentognathic evidence[END_REF][START_REF] Martinón-Torres | Dental remains from Dmanisi (Republic of Georgia): morphological analysis and comparative study[END_REF][START_REF] Zanolli | Additional evidence for morpho-dimensional tooth crown variation in a New Indonesian H. erectus sample from the Sangiran Dome (Central Java)[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF].

However, the high-degree of morphological variation of the Asian dental remains [START_REF] Zanolli | Brief communication: Two human fossil deciduous molars from the Sangiran Dome (Java, Indonesia): Outer and inner morphology[END_REF](Zanolli et al., , 2018a;;[START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF][START_REF] Xing | Late Middle Pleistocene hominin teeth from Tongzi, southern China[END_REF] suggests the (co)existence of derived (modern H. sapiens-like) and archaic morphotypes in the Middle Pleistocene. This complex scenario gives rise to a number of questions regarding the routes of dispersal, possible interactions and gene flow between lineages in the mid-Middle to early Late Pleistocene, before the attested presence of anatomically modern humans [START_REF] Liu | The earliest unequivocally modern humans in southern China[END_REF].

In the late Pleistocene Chinese hominin fossil record, the taxonomic status and phylogenetic relationships with earlier, contemporaneous and later human lineages remains unclear. Part of this is probably due to the limited sample size and the heterogeneity of expression of archaic and derived traits in the mid-to late Middle Pleistocene 'non-H. erectus' human record. For instance, the late Middle Pleistocene dental remains from Panxian Dadong, even if limited to four isolated teeth, retain some archaic traits in the anterior dentition, including a weakly expressed finger-like tuberculum dentale on the upper central incisor and robust incisor and canine crowns, in association with derived features in most teeth (symmetrical canine crown without cingulum, modern-like premolar crown morphology; [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF]. Similarly, the late Middle Pleistocene teeth from Tongzi, even if penecontemporaneous with late representatives of H. erectus, differ from the typical condition seen in the latter group and show a series of derived traits that are more commonly seen in Late Pleistocene Homo, including: the absence of tuberculum dentale on the upper incisor, simplified and symmetrical premolar crowns without grooves and ridges on the buccal aspect, and mesiodistally compressed upper molars [START_REF] Wu | New discoveries of human fossil in Tongzi, Guizhou[END_REF][START_REF] Xing | Late Middle Pleistocene hominin teeth from Tongzi, southern China[END_REF]. Conversely, dental remains from other late Middle to early Late Pleistocene populations such as Chaoxian, Changyang and Xujiayao express a degree of archaic, H. erectus-like features-see external crown morphology in [START_REF] Xing | Hominin teeth from the early Late Pleistocene site of Xujiayao, Northern China[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF] and [START_REF] Bailey | A comparative dental metrical and morphological analysis of a Middle Pleistocene hominin maxilla from Chaoxian (Chaohu), China[END_REF], and root canal configuration in [START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF]. Taphonomic factors and occlusal wear have also erased parts of the external tooth morphology, thus adding to the confusion concerning the evolutionary interpretations of these specimens [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Xing | Hominin teeth from the early Late Pleistocene site of Xujiayao, Northern China[END_REF][START_REF] Xing | Late Middle Pleistocene hominin teeth from Tongzi, southern China[END_REF][START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF]. However, even if the external morphology is altered, diagnostic taxonomic and phylogenetic information may still be preserved in the internal structure of teeth, and their internal signature can still be extracted using the methods developed in virtual paleoanthropology (Skinner et al., 2008b;[START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF][START_REF] Weber | Virtual anthropology[END_REF][START_REF] Zanolli | Inner tooth morphology of Homo erectus from Zhoukoudian. New evidence from an old collection housed at Uppsala University, Sweden[END_REF][START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF] .

We used X-ray µCT scanning and 3D image processing techniques to investigate the inner structure of six late Middle Pleistocene premolar specimens from the sites of Changyang (PA76, PA81), Tongzi (PA873), Panxian Dadong (PA 1577 andPA 1578) and Chaoxian (Chaoxian), whose taxonomic attribution is controversial. In order to characterize these specimens, we focus here on the enamel-dentine junction (EDJ) morphology. This interface results from the membrana praeformativa that configures the outer enamel surface (OES) during tooth growth [START_REF] Corruccini | The dentinoenamel junction in primates[END_REF]. In mature teeth, it preserves diagnostic morphological information, even when the occlusal enamel is worn. It is thus recognized as a reliable taxonomic proxy to distinguish between hominid taxa (even at the subspecies level) and can also bring hominid phylogenetic insights [START_REF] Weber | Virtual anthropology (VA): A call for Glasnost in paleoanthropology[END_REF][START_REF] Skinner | Enamel-dentine junction morphology of extant hominoid and fossil hominin lower molars[END_REF]Skinner et al., 2008a;[START_REF] Braga | The enamel-dentine junction in the postcanine dentition of Australopithecus africanus: intra-individual metameric and antimeric variation[END_REF][START_REF] Macchiarelli | From outer to inner structural morphology in dental anthropology: integration of the third dimension in the visualization and quantitative analysis of fossil remains[END_REF][START_REF] Pan | Further morphological evidence on South African of[END_REF][START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF].

Comparing with a number of fossil and modern human specimens/samples, including Chinese early to mid-Middle Pleistocene H. erectus from the sites of Zhoukoudian, Hexian, Xichuan and Yiyuan, the aims of this study are: (1) to examine the ability of premolar EDJ morphological information to discriminate among a selection of fossil and modern human premolar samples; (2) to document the morphological variability of mid-to late Middle Pleistocene human premolars in East Asia; and (3) to re-evaluate the taxonomic affinities of six late Middle Pleistocene specimens.

Materials and methods

Study sample

The Chinese Middle Pleistocene hominin sample consists of two groups; a map of their location is shown in Supplementary Online Material (SOM) Figure S1. The first group includes early to mid-Middle Pleistocene H. erectus (HEC, n = 12) from the sites of Zhoukoudian, Hexian, Xichuan and Yiyuan, with a chronological range from 800 to 282 ka (Table 1). Note that the specimens from Xichuan County, Hubei Provence lack accurate stratigraphic context, but their morphological signature a priori places them close to Zhoukoudian H. erectus [START_REF] Wu | Human fossil teeth from Xichuan, Henan[END_REF][START_REF] Wu | Human evolution in China: a metric description of the fossils and a review of the sites[END_REF]. The allocation of Xichuan hominins to East Asian H. erectus is further supported by a recent comparative study [START_REF] Xing | The fossil teeth of the Peking Man[END_REF] on the dental inner structure of Chinese Middle Pleistocene H. erectus. The second group is represented by six premolars from three mid-to late Middle Pleistocene sites (Changyang, Tongzi, Panxian Dadong and Chaoxian). These specimens do not have a clear taxonomic attribution, as in previous studies they were classified as archaic or premodern H. sapiens, or were categorized as non-H. erectus Homo [START_REF] Xing | Late Middle Pleistocene hominin teeth from Tongzi, southern China[END_REF] (and are thus regarded here a priori as Homo sp. from China, HSPC). We compared the Chinese specimens with 56 fossil and extant human specimens (Table 1), including Javanese Early Pleistocene H. erectus (HEJ) specimens from Sangiran [START_REF] Grine | Fossil hominid teeth from the Sangiran Dome (Java, Indonesia)[END_REF][START_REF] Kaifu | Taxonomic affinities and evolutionary history of the Early Pleistocene hominids of Java: Dentognathic evidence[END_REF]; North African late Early Pleistocene hominins (HNA) from Tighenif [START_REF] Zanolli | Endostructural characterization of the H. heidelbergensis dental remains from the early Middle Pleistocene site of Tighenif, Algeria[END_REF]; Neanderthals (NEA) from Krapina and La Chaise [START_REF] Macchiarelli | How Neanderthal molar teeth grew[END_REF]NESPOS database, 2020), and modern H. sapiens (MHS) from Africa and Eurasia (Pan et al., 2017).

Methods

Tooth wear stages were scored according to the standard proposed by [START_REF] Smith | Patterns of molar wear in hunter-gatherers and agriculturalists[END_REF].

For the Neanderthal and modern H. sapiens specimens, we selected only premolars without dentine exposure on the occlusal surface. However, for the Chinese Middle Pleistocene specimens, because of the limited available sample, we included all specimens with wear degree less than stage 3 ("full cusp removal and/or moderate dentine patches"; Smith, 1984:45).

Information on the crown outer and inner morphology-non-metric features, 2D geometric morphometrics (GM) on the OES-has been reported for most of the HEC and HSPC specimens in previous studies [START_REF] Black | Fossil Man in China. The Choukoutien cave deposits with a synopsis of our present knowledge of the late Cenozoic of China[END_REF][START_REF] Chia | Notes on the human and some other mammalian remains from Changyang, Hupei[END_REF][START_REF] Wu | Human fossil teeth from Xichuan, Henan[END_REF][START_REF] Chen | Uranium series dating of fossil bones from Hexian and Chaoxian fossil human sites[END_REF][START_REF] Bailey | A comparative dental metrical and morphological analysis of a Middle Pleistocene hominin maxilla from Chaoxian (Chaohu), China[END_REF][START_REF] Xing | Morphological variation of Zhoukoudian H. erectus teeth[END_REF][START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF][START_REF] Xing | Late Middle Pleistocene hominin teeth from Tongzi, southern China[END_REF][START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF]Zanolli et al., 2018a). The dental remains from Xichuan (PA524, PA525, PA526, PA527, PA528) and Changyang (PA76, PA81) are reported and their internal structure (as well as that of the Chaoxian specimen) analyzed here for the first time. The terms employed for the morphological descriptions are from [START_REF] Weidenreich | The dentition of Sinanthropus pekinensis: a comparative odontography of the hominids[END_REF], [START_REF] Turner | Scoring procedures for key morphological traits of the permanent dentition: the Arizona State University dental anthropology system[END_REF] and [START_REF] Scott | The Anthropology of Modern Human Teeth: Dental Morphology and its Variation in Recent Human Populations[END_REF].

Microcomputed tomography and processing of the 3D data All Chinese fossil specimens (except for two Zhoukoudian H. erectus specimens housed at Uppsala University) were scanned using a 225 kV micro-XCT scanner at the Institute of Vertebrate Paleontology and Paleoanthropology, China according to the following parameters: 120-140 kV voltage, 100-120 mA current, angular step of 0.5° over a scan angle of 360°. The final volumes were reconstructed with an isotropic voxel size ranging from 40 to 70 µm. The fossil premolars PMU M3887 and M3549, housed at Uppsala University, were scanned at the Multidisciplinary Laboratory of the The Abdus Salam International Centre for Theoretical Physics (ICTP), Trieste, Italy (Zanolli et al., 2018a).

Image stacks were imported into Avizo 9.0 (FEI Visualization Sciences Group), and the images were segmented using semiautomatic procedures and an adaptation of the half-maximum height method [START_REF] Fajardo | Assessing the accuracy of high-resolution X-ray computed tomography of primate trabecular bone by comparisons with histological sections[END_REF][START_REF] Coleman | CT thresholding protocols for taking measurements on three-dimensional models[END_REF]; the EDJ surfaces were generated using the 'unconstrained smoothing' option. For 11 out of 18 Chinese human fossil specimens (marked in Table 1), some cusp tips were worn, so we manually reconstructed the dentine horn tips based on the curvature and orientation of the EDJ at each side of the dentine horns (Smith et al., 2012;Zanolli et al., 2018b). Our reconstructions are depicted in SOM Figure S2.

Geometric morphometric analysis We used a recently developed diffeomorphic surface matching (DSM) method to compare the EDJ shapes among hominin groups. This landmark-free, mesh-based approach relies on the construction of average surface models, and the difference between surfaces is interpreted as the amount of deformation needed to align them by using diffeomorphic shape matching [START_REF] Glaunes | Template estimation form unlabeled point set data and surfaces for computational anatomy[END_REF][START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. Then deformations between surfaces are mathematically modeled as smooth and invertible functions (i.e., diffeomorphisms; [START_REF] Glaunes | Template estimation form unlabeled point set data and surfaces for computational anatomy[END_REF][START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. From a set of surfaces, an atlas of surfaces is created.

The method estimates an average object configuration or mean shape from a collection of object sets, here the EDJ surfaces, and computes the deformations from the mean shape to each specimen. In addition, a set of initial control points located near its most variable parts, and a set of momenta parameterizing the deformations of the mean shape to each individual are estimated (Durrleman et al., 2012a, b;Beaudet et al., 2016a).

This method has been successfully employed in a number of paleontological/paleoanthropological investigations (Beaudet et al., 2016b[START_REF] Beaudet | The endocranial shape of Australopithecus africanus: surface analysis of the endocasts of Sts 5 and Sts 60[END_REF][START_REF] Zanolli | Inner tooth morphology of Homo erectus from Zhoukoudian. New evidence from an old collection housed at Uppsala University, Sweden[END_REF][START_REF] Braga | Efficacy of diffeomorphic surface matching and 3D geometric morphometrics for taxonomic discrimination of Early Pleistocene hominin mandibular molars[END_REF][START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF][START_REF] Urciuoli | The evolution of the vestibular apparatus in apes and humans[END_REF]. A recent study comparing the discriminative power of 3D landmark-based GM and DSM demonstrated that DSM has the potential to better capture the geometric details that describe tooth shape [START_REF] Braga | Efficacy of diffeomorphic surface matching and 3D geometric morphometrics for taxonomic discrimination of Early Pleistocene hominin mandibular molars[END_REF].

To perform this analysis, as a previous step, at each dental position the EDJ surfaces were manually oriented, then superimposed using the rigid and uniform scale option of the 'align surface' module in Avizo 9.0 (Visualization Sciences Group Inc.). This was done by minimizing the root mean square distance between the points of each specimen to corresponding points on the reference surface using an iterative closest point algorithm.

We used the Deformetrica software (https://www.deformetrica.org; [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF] to generate a global mean shape (GMS) with a set of diffeomorphisms from the GMS to each individual and the output (control points and deformation momenta) used to perform the statistical analyses to explore the EDJ shape variation. An example of the original EDJ surface, the simplified EDJ surface imported in Deformetrica v. 4.3 and the output virtual models generated by DSM process is shown in SOM Figure S3 (see also SOM Fig. S4 for the output surfaces of all the Chinese fossil specimens after deformations). Using the packages ade4 1.7-13 (Dray and Dufour, 2007) and Morpho v. 2.6 [START_REF] Schlager | Retrodeformation of fossil specimens based on 3D bilateral semi-landmarks: Implementation in the R package "Morpho[END_REF] for R v. 3.5.1 (R Development Core Team, 2012), we checked the shape variation of the fossil and recent human groups in a principal component analysis (PCA). Even if the first two components of the PCA account for approximately 40 to 65% of the total variance in these analyses, there is still some information in the following components. Howerver, visualizing and interpreting a complex space in multiple dimensions is difficult. We computed between-group principal component analyses (bgPCA) to evaluate the degree of shape variation between hominin groups. We used both PCA and bgPCA as the former is mostly a means to explore the morphospace, with the variance spread over multiple axes, whereas bgPCA, through projection of the data onto the principal components of the group averages, represents a viable alternative to discriminant analyses (Mitteroecker and Bookstein, 2011). Nonetheless, when sample size is small and groups are few, the results of bgPCA must be interpreted with care and a number of tests are needed to confirm that group separation is not spurious [START_REF] Bookstein | Pathologies of between-groups principal components analysis in geometric morphometrics[END_REF][START_REF] Cardini | Seeing distinct groups where there are none: spurious patterns from between-group PCA[END_REF].

As a result, we checked and confirmed that group separation preexisted in the PCA and bgPCA was used to assess group affinity of target specimens. Following the recommendations of [START_REF] Cardini | Cross-validated between group PCA scatterplots: A solution to spurious group separation?[END_REF], we also computed the cross-validated bgPCA (CV-bgPCA) to confirm that the group distribution is not spurious. The six Chinese late Middle Pleistocene Homo sp. specimens were then projected a posteriori in the bgPCA plot to assess their morphological affinities.

The residual of allometric signals were tested using the coefficient of determination (R 2 ) of a multiple regression [START_REF] Bookstein | Morphometric Tools for Landmark Data: Geometry and Biology[END_REF] in which the explicative variable is the EDJ area as proxy for size and the dependent variables are the (bg)PC scores.

Results

Description of enamel-dentine junction morphologies

Here we describe the morphological features on the EDJ of previously unreported material from Xichuan, Changyang and Chaoxian and compared it with a number of Pleistocene and Holocene hominin specimens (see Figs. 1A-V and 2A-S for EDJ surfaces in occlusal view). The Chinese fossil specimens are ordered and described below with regard to dental position.

Left maxillary third premolar (P 3 , PA524, Xichuan) See Figure 1B. Viewed from the occlusal surface, the crown contour is oval and slightly asymmetrical, the width of the buccal aspect is slightly larger than the lingual aspect, the transverse crest is pronounced and connects the buccal and lingual cusps, forming a small and deep anterior fovea and a larger posterior fovea, and a small tubercle is visible at the mesial aspect of the buccal cusp base. Along the marginal ridge, no accessory ridge is identified.

Left maxillary third premolar (P 3 , PA76, Changyang) See Figure 1D. The occlusal outline of the marginal ridges is elliptical and slightly asymmetrical, two faint accessory ridges (one of which is bifurcated at the end) can be seen at the base of the buccal cusp but interrupted in the middle of the small occlusal basin, and no continuous transverse crest is formed. One mesial accessory ridge is identified. Two weakly developed tubercles are visible at the base of the lingual cusp.

Right maxillary fourth premolar (P 4 , Chaoxian) See Figure 1O. This P 4 is attached to a maxilla, whose right P 3 -M 1 and left P 4 -M 2 are preserved [START_REF] Bailey | A comparative dental metrical and morphological analysis of a Middle Pleistocene hominin maxilla from Chaoxian (Chaohu), China[END_REF]. The occlusal basin of this P 4 is relatively small compared with the crown base, similar to the condition in H. erectus (Fig. 1M, N, P, Q). The crown contour is ovorectangular.

The transverse crest is not well expressed and is interrupted by the occlusal basin. In all, compared with the highly crenulated EDJ surface that is expressed at East Asian H. erectus (e.g., [START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF] this study), the EDJ of Chaoxian P 4 is more similar to the less complex morphology of Indonesian H. erectus from Sangiran. Weaklyexpressed accessory tubercles are found at the mesial and distal part of the occlusal basin.

Left mandibular third premolar (P3, PA526, Xichuan) See Figure 2C. The crown outline of PA526 is asymmetrical and mesiodistally elongated in occlusal view. The transverse crest is well developed and continuous. The mesial fovea is slightly smaller than the distal one, as the latter is born by a wide talonid. An accessory ridge can be found at the mesial aspect of the lingual cusp, but it does not reach the mesial marginal ridge. At the buccal surface, the mesial vertical groove is better developed than the distal one.

Right mandibular third premolar (P3, PA527, Xichuan) See Figure 2D. The morphology of PA527 is very similar to that of PA526; however, instead of an accessory ridge, a fainted wrinkle is present at the mesial side of the lingual cusp. Like PA526, the mesial vertical groove at the buccal surface is better expressed than the distal one.

Right mandibular fourth premolar (P4, PA525, Xichuan) See Figure 2L. The outline of the marginal ridge is ovorectangular with a protruding mesiobuccal angle. The transverse crest is moderately expressed and interrupted. At the occlusal basin, two small, weakly expressed wrinkles are presented at the mesial marginal ridge and the distal side of the lingual cusp. Moreover, a weakly expressed accessory cusp is seen at the distal side of the buccal cusp. On the buccal surface, two vertical grooves are clearly detected.

Left mandibular fourth premolar (P4, PA528, Xichuan) See Fig. 2M. The morphology of PA528 is slightly different from that of PA525, but resembles that of the H. erectus specimen M3887 from Zhoukoudian. The occlusal outline of its marginal ridge is also ovorectangular. The transverse crest is well developed and uninterrupted, and the mesial fovea is clearly smaller than the distal fovea. Several accessory cusps can be seen at the distal marginal ridge, and a moderately expressed accessory crest is present at the distal side of the lingual cusp. A moderately expressed wrinkle is visible at the mesial marginal ridge. On the buccal surface, two vertical grooves are well expressed.

Left mandibular fourth premolar (P4, PA81, Changyang) See Figure 2N. The crown contour is ovorectangular, the distal part of the crown being wider than the mesial one. The transverse crest is interrupted by a deep groove and is weakly developed.

Mesial and distal accessory ridges for the lingual cusp can be identified, three accessory cusps are seen at the distal part of the fovea posterior, and an accessory ridge is present distally to the buccal cusp and deflects at the center of the occlusal basin. On the buccal surface, the distal vertical groove is more pronounced than the mesial one.

The East Asian mid-Middle Pleistocene H. erectus (PA67,68,110,524,525,526,527,528,832,M3887 and Shy071) consistently show a complex morphology, with vertical ridges and crests on the buccal aspect of the EDJ, together with a wrinkled occlusal basin and accessory cusps and crests. In the latter half of the Middle Pleistocene, more variation is observed.

For the upper premolars, the P 3 PA76 (Fig. 1D) from Changyang retains archaic traits, as it exhibits multiple accessory crests running in the occlusal basin. In comparison with PA76, the P 3 PA873 from Tongzi (Fig. 1E) displays a simple morphology, with a proportionally larger and smoother occlusal basin and three faint occlusal crests running on the internal aspect of the paracone. The mesial and distal marginal ridges of PA873 are more protruding, with numerous small cuspule-like features. A similarly simplified and enlarged occlusal basin than in PA873 is also shown in the P 3 PA1577 (Fig. 1F) from Panxian Dadong. The P 4 from Chaoxian (Fig. 1O) has a smaller occlusal basin compared with the crown base, similar to East Asian H. erectus from Zhoukoudian (PA68, Fig. 1M; Shy071, Fig. 1N). However, unlike the 'dendritelike' EDJ that is shown in H. erectus, the occlusal surface relief of Chaoxian is less complex (but is not comparable to PA873 and 1577, which are more simple and smooth) with a weakly-expressed, interrupted transverse crest.

As for the lower premolars, a previous study on the crown morphology [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF] showed that the P3 PA1578 (Fig. 2E) from Panxian Dadong exhibits weakly expressed archaic traits and a derived general configuration with respect to the condition in H. erectus. The EDJ outline is more symmetrical in PA1578 than in Early to mid-Middle Pleistocene specimens from Africa and Asia (including the North African specimen Tighenif 2, Indonesian and Chinese H. erectus), but is more asymmetrical than in recent modern humans. Previous analyses of crown geometric morphometrics place the Panxian Dadong premolars close to modern humans [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF]. Studies of the root structure show that, compared with modern humans, PA 1578 has a robust morphology, although it is more gracile than the root complex usually found in H. erectus [START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF][START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF]. On the other hand, the P4 PA81 (Fig. 2N) from Changyang exhibits a wrinkled occlusal basin, with vertical ridges expressed on the buccal aspect of the EDJ (similar to H. erectus; Fig. 2K-M). The talonid of PA81 is mesiodistally elongated with a minute hypoconid dentine horn (Fig. 2N). In sum, from the six premolar specimens dated to late Middle Pleistocene, despite a small sample size, greater morphological variability (i.e., primitive, H. erectus like and derived, modern-human like) is observed compared with East Asian H. erectus.

Geometric morphometric analyses

The results of GM analyses using DSM are shown in Figures 34(see also SOM Fig. S4, which illustrates the output surfaces after DSM analyses). Taking into account that some authors have issued cautions about the interpretation of bgPCA results [START_REF] Bookstein | Pathologies of between-groups principal components analysis in geometric morphometrics[END_REF][START_REF] Cardini | Seeing distinct groups where there are none: spurious patterns from between-group PCA[END_REF], we present in Figures 3 and4 the results from both PCA and bgPCA. Results of CV-bgPCA are illustrated in the SOM Figure S5, showing that the groups are separated in a similar way than in the bgPCA (even if inversion of the distribution of points between the negative and positive space along the axes may occur) and that group separation is not spurious. A mild allometric signal is observed along some of the (bg)PCs (SOM Table S1). While size-related component is present, the variation observed in the PCA and bgPCA thus mostly represents shape differences.

PCA and bgPCA yielded consistent results with H. erectus, Neanderthals and modern humans roughly distributed in a similar way along the first two axes, even if as expected the bgPCA better separates them (Figs. 3 and4). In any case, even if slightly more overlapping, the groups identified when observing the morphospace obtained with bgPCA are already present in the PCA. The extreme shape changes along each (bg)PC are shown at the extremities of each axis. Even if comparative late Early to mid-Middle Pleistocene, Late Pleistocene and Holocene Homo groups are closer in the PCA of P 4 (Fig. 3C), H. erectus (including the North African specimen Tighenif 2, Indonesian and Chinese H. erectus) is discriminated from later groups along (bg)PC1 (Figs. 3 and4). (bg)PC1 reflects a chronoevolutionary trend, with lower EDJ in H. erectus than in the two other groups, as shown by the extreme shapes at the end of the axes (Figs 3 and4; see also Figs. 1 and2). For P 3 and P 4 , H. erectus shows a narrower occlusal basin (buccolingually compressed and small with respect to the crown base) and more internally-positioned dentine horns than in Neanderthals and modern humans (Fig. 3 A-D). As seen in the bgPCA, Neanderthal upper premolars differs from those of modern humans by their mesiodistally wider EDJ and deeper occlusal basin (Fig. 3B,D). For the P3 EDJ, as seen along (bg)PC1, H. erectus exhibits a proportionally more expanded talonid (resulting in a protruding distolingual angle of the crown) than in Neanderthals and modern humans that show more equally developed mesial and distal foveae. Along bgPC2, Neanderthal P3 are separated from those of modern humans due to their smaller occlusal basin (Fig. 4A-B). PC1 and bgPC1 of the P4 EDJ analyses show that H. erectus has low lateral aspects and short dentine horns, with a more buccolingually compressed occlusal basin (and more internally-set metaconid) than in Neanderthals and modern humans (Fig. 4 C,D). bgPC2 discriminates Neanderthals from modern humans, the latter having a lingually-positioned metaconid and buccolingually expanded occlusal basin and more equally developed mesial and distal foveae than in the former (Fig. D). In this context, the Changyang specimens (PA76 P 3 and PA81 P4) and the P 4 from Chaoxian unambiguously fall close to the late Early Pleistocene to mid-Middle Pleistocene human groups, showing low EDJ, large crown base and complex surface relief with relatively small occlusal basin-as supported by previous observations on the EDJ surfaces (Figs. 1D, O and2N). Conversely, the P 3 PA873 from Tongzi, P 3 PA1577 and P3 PA1578 from Panxian Dadong are placed within or very close to the variability of Neanderthals and modern humans (Figs. 3A, B and4A, B).

Results of the cross-validation tests performed on the bgPC scores shows that predictive accuracy is robust (Table 2). Correct classification of 88-100% is obtained in most analyses, except for the P 4 analysis, which exhibits slightly more confusion for Neanderthals and modern humans (80% of correct classification; but H. erectus teeth show 100% discrimination). The P 3 PA76, the P 4 Chaoxian and the P4 PA81 are classified as H. erectus. The P 3 PA873 and the P3 PA1578 are both classified nearly equally as belonging to either modern humans or Neanderthals, while the P 3 PA1577 is classified as a Neanderthal.

Discussion

Increasing evidence from the morphological study of fossil hominins from China suggests remarkable morphological variability toward the end of the Middle Pleistocene [START_REF] Bailey | A comparative dental metrical and morphological analysis of a Middle Pleistocene hominin maxilla from Chaoxian (Chaohu), China[END_REF][START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Chang | The first archaic Homo from Taiwan[END_REF][START_REF] Martinón-Torres | A "source and sink" model for East Asia? Preliminary approach through the dental evidence[END_REF]. Similarly, recent advances in the fields of paleoproteomics and paleogenetics confirm the diversity of hominin groups that coexisted in this chronogeographic interval, including a sister group of Neanderthals: the Denisovans [START_REF] Chen | A late Middle Pleistocene Denisovan mandible from the Tibetan Plateau[END_REF]. However, taxonomic attributions of late Middle to early Late Pleistocene specimens and their interpretation in an evolutionary perspective remain unclear [START_REF] Bae | The late Middle Pleistocene hominin fossil record of eastern Asia: synthesis and review[END_REF][START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Kaifu | Archaic hominin populations in Asia before the arrival of modern humans: their phylogeny and implications for the "Southern Denisovans[END_REF]. For the first time, through a novel geometric morphometric approach we analyzed the EDJ shape of a sample of hominin premolars from China spanning the Middle Pleistocene and compared this sample with other fossil and recent human groups.

At a large scale, our study reported at least two morphotypes in EDJ shape expressed in this time period: one displaying a H. erectus-like condition and the other showing a morphology that is closer to Late Pleistocene humans, sharing traits with both Neanderthals and modern humans. Irrespective of the taxon, in this timeframe, we also document a general trend towards a reduction of the crown base in relation to the occlusal basin, narrowing of the lingual aspect, and increasing of EDJ relief.

From the late Early to mid-Middle Pleistocene, we found a proportionally broadened crown base associated with a smaller occlusal basin and low EDJ relief in hominin premolars from Indonesia (Sangiran), Algeria (Tighenif) and China (Zhoukoudian, Hexian, Xichuan and Yiyuan). This is also compatible with previous studies suggesting a possible retention of Early Pleistocene hominin features (as seen in Sangiran 4 for example) in Chinese Middle Pleistocene H. erectus [START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF][START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Zanolli | Inner tooth morphology of Homo erectus from Zhoukoudian. New evidence from an old collection housed at Uppsala University, Sweden[END_REF]. Thus, until the mid-Middle Pleistocene, only one human group, likely related to H. erectus can be identified in China from the available evidence. Toward the late Middle Pleistocene, the Chinese fossils express at least two morphotypes. The first one is characterized by large crown base, centrallypositioned dentine horn tips and low EDJ topography, together with a set of crowncomplexifying non-metric traits as seen in the specimens from Changyang (PA76 and 81) and Chaoxian. It is noteworthy that the Chaoxian P 4 expresses H. erectus-like features at the EDJ level that are in agreement with precious observations made on the outer crown morphology [START_REF] Bailey | A comparative dental metrical and morphological analysis of a Middle Pleistocene hominin maxilla from Chaoxian (Chaohu), China[END_REF]. In light of these results, together with the similarities in EDJ shape of these specimens with H. erectus (Figs. 1 and2), we interpret the Changyang and Chaoxian hominins as a late representative of an East Asian H. erectus lineage.

The second morphotype, represented by the fossils from Panxian Dadong and Tongzi, is characterized by a series of derived features generally found in Late Pleistocene hominins, like Neanderthals and modern humans, as already suggested by previous studies of the external tooth morphology [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF][START_REF] Zanolli | Evidence for increased hominid diversity in the Early-Middle Pleistocene of Java, Indonesia[END_REF]. The EDJ shape of PA873 and PA1578, even if intermediate between the Neanderthal and modern human conditions, tends to show more affinities with the Neanderthal general pattern. Even if there is no evidence for such an early onset of dispersion into Asia, the recent finding from Misliya in Israel [START_REF] Hershkovitz | The earliest modern humans outside Africa[END_REF] indicates that we cannot rule out the possibility that they may represent a biogeographical variation of early modern humans in China. It is also possible that these non-H. erectus Homo may sample a distinct lineage, possibly Denisovans. In fact, predictions from paleogenetics that Denisovans were a group of Neanderthal-like humans that lived from the Altai to East and South Asia [START_REF] Meyer | A high-coverage genome sequence from an archaic Denisovan individual[END_REF][START_REF] Meyer | Nuclear DNA sequences from the Middle Pleistocene Sima de los Huesos hominins[END_REF][START_REF] Jacobs | Timing of archaic hominin occupation of Denisova Cave in southern Siberia[END_REF][START_REF] Mondal | Approximate Bayesian computation with deep learning supports a third archaic introgression in Asia and Oceania[END_REF], were recently reaffirmed by the combined morphological and molecular analyses of a mandible fragment from Xiahe, on the Tibetan Plateau [START_REF] Chen | A late Middle Pleistocene Denisovan mandible from the Tibetan Plateau[END_REF]. As archaic Denisovans probably occupied a broad geographical range, admixture between Denisovans and early modern humans probably occurred not only in Siberia and Southeast Asia but also in Southern China (Reich et al., 2011;[START_REF] Cooper | Did the Denisovans cross Wallace's line?[END_REF]. Even if no Denisovan premolar is known yet, given the current evidence, it is possible that the Panxian Dadong and Tongzi hominins could belong to this taxon, as it was also proposed for other late Middle Pleistocene specimens from Asia, like the human remains from Xujiayao [START_REF] Xing | Hominin teeth from the early Late Pleistocene site of Xujiayao, Northern China[END_REF][START_REF] Chen | A late Middle Pleistocene Denisovan mandible from the Tibetan Plateau[END_REF] and Penghu [START_REF] Chang | The first archaic Homo from Taiwan[END_REF]. Another possibility is that these specimens sample representative from both Denisovans and Pleistocene modern humans.

Conclusions

The present study suggests the existence of time-related evolutionary changes in the premolar EDJ morphology in a number of Early and Middle Pleistocene hominin groups, with special reference to mid-and late Middle Pleistocene populations in China.

A general trend towards a reduction of crown base in relation to occlusal aspects, narrowed lingual aspects, heightened EDJ relief and a surface simplification is found through time. The current analysis expands our knowledge on the dental inner morphology of Middle Pleistocene hominins, supports the coexistence of archaic ('Homo erectus-like') and derived ('Late Pleistocene Homo-like') dental morphological patterns and, through a dental perspective, suggest the preservation of plesiomorphic, H. erectus-like morphology during the Middle Pleistocene. New fossil discoveries and/or molecular analyses will be critical to explore further the past diversity and to shed light on the dispersal and interaction of hominin populations in East Asia.

earliest Homo lower postcanine dentition: enamel thickness and enamel dentine junction. Journal of Human Evolution 96, 82-96. Pan, L., Thackeray, J.F., Dumoncel, J., Zanolli, C., Oettlé, A., De Beer, F., Hoffman, J., Duployer, B., Tenailleau, C., Braga, J., 2017. Intra-individual c Worn crown, reconstructed dentine horn tips. For the reconstruction see SOM Figure S1.

Table 2

Results of the cross-validation analyses performed on the between-group principal component scores.

The posterior probabilities that the Chinese specimens are closer to any group were then calculated. The highest probabilities are indicated in bold (the greater the number, the higher the probability).

Abbreviations

  metameric variation expressed at the enamel-dentine junction of lower post-canine dentition of South African fossil hominins and modern humans. American Journal of Physical Anthropology 163, 806-815. Pan, L., Zanolli, C., 2019. Comparative observations on the premolar root and pulp canal configurations of Middle Pleistocene Homo in China. American Journal of Physical Anthropology 168, 637-646. Pan, L., Dumoncel, J., Mazurierd, A., Zanolli, C., 2019. Structural analysis of premolar roots in Middle Pleistocene hominins from China. Journal of Human Evolution 136, 102669. R Development Core Team, 2012. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. Radovčić, J., Smith, F.H., Trinkaus, E., Wolpoff, M.H., 1988. The Krapina Hominids: An Illustrated Catalog of Skeletal Collection. Mladost and the Croatian Natural History Museum, Zagreb. Reich, D., Patterson, N., Kircher, M., Delfin, F., Nandineni, M.R., Pugach, I., Ko, A.M.-S., Ko, Y.-C., Jinam, T.A., Phipps, M.E., 2011. Denisova admixture and the first modern human dispersals into Southeast Asia and Oceania. The American Journal

20 a

 20 Age references: A = Yuan et al. (1986); B = Grün et al. (1998); C = Grün et al. (1997); D = Huang et al. (1993); E = Shen et al. (2001); F = Han et al. (2016); G = Shen and Jin (1991); H = Karkanas et al. (2008); I = Chen et al. (1987); J = Shen et al. (2010); K = Huffman (2001); L = Matsu'ura et al. (2020); M = Sahnouni and Van der Made (2009); N = Rink et al. (1995); O = Macchiarelli et al. (2006). b Citations reporting descriptions of the morphology of the specimens : 1 = Xing (2012); 2 = Xing et al. (2018); 3 = Black et al. (1933); 4 = Zanolli et al. (2018a); 5 = Xing et al. (2014); 6= Liu et al. (2017); 7 = Chen et al. (1987); 8= Wu and Wu (1982); 9 = Xing et al. (2016); 10 = Wu and Wu (1982); 11 = Chia (1957); 12 = Xing et al. (2019); 13 =Liu et al. (2013); 14 = Bailey and Liu (2010); 15 = Grine (1994); 16 = Zanolli and Mazurier (2013); 17 = Radovčić et al. (1988); 18 = NESPOS database (2020); 19 = Macchiarelli et al. (2006); 20 = Pan et al. (2017).

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. Volume rendering of lower premolar EDJ of Homo erectus from China (HEC),

Figure 3 .

 3 Figure 3. Principal component analysis (PCA) and between-group principal component

Figure 4 .

 4 Figure 4. Principal component analysis (PCA) and between-group principal component

Table 1

 1 Composition of the study specimens from China and of the comparative samples.

	Chinese Pleistocene specimens	P 3 P 4 P 3 P 4	Provenance	Chronological range a	Citations b
	Middle Pleistocene Homo erectus from China (HEC)
	PA67	1	Zhoukoudian Loc. 1 (Layer 11)	800-770 (A)	1-3
	PA68	1	Zhoukoudian Loc. 1 (Layer 11)	800-770 (A)	1-3
	PMU M3549	1	Zhoukoudian Loc. 1 (Layer 5)	550-500 Ka (B)	2, 4
	PMU M3887	1	Zhoukoudian Loc. 1 (Layer 5)	550-500 Ka (B)	2, 4
	PA110 c	1	Zhoukoudian Loc. 1 (Layer 3)	320-282 ka (C, D); 500 ka (E)	1-3
	PA832 c	1	Hexian	437-387 ka (B)	5, 6
	Shy071 c	1	Yiyuan	420-320 ka (F)	7-9
	PA524 c	1	Xichuan	Unknown	10
	PA525 c	1	Xichuan	Unknown	10
	PA526 c	1	Xichuan	Unknown	10
	PA527 c	1	Xichuan	Unknown	10
	PA528	1	Xichuan	Unknown	10

  : HE = Homo erectus from Asia and late Early Pleistocene Homo from North Africa; NEA = Neanderthals; MHS = modern Homo sapiens.

	Dental position Specimen	Predictive accuracy of comparative specimens	Classification of the investigated Chinese fossils
	P 3	PA76	HE: 100%; NEA: 100%, MHS: 100%	HE: 0.061; NEA: <0.001; MHS <0.001
		PA873	HE: 100%; NEA: 100%, MHS: 100%	HE: <0.001; NEA: 0.001;MHS: 0.002
		PA1577	HE: 100%; NEA: 100%, MHS: 100%	HE: <0.01; NEA: 0.104; MHS: <0.001
	P 4	Chaoxian	HE: 100%; NEA: 80%, MHS: 80%	HE: 0.334; NEA: 0.001; MHS: <0.001
	P 3	PA1578	HE: 100%; NEA: 100%, MHS: 91%	HE: <0.001; NEA: 0.004; MHS: 0.026
	P 4	PA81	HE: 100%; NEA: 100%, MHS: 88%	HE: 0.808; NEA: <0.001; MHS: <0.001

Acknowledgements

This study is dedicated to the memory of our dear friend and coworker Pan Lei, who passed away while this paper was being peer-reviewed. For access to specimens, we sincerely thank C. Argot, O.

Kullmer, E.N. L'Abbé, H. Lelièvre, W. Liu, A. Oettlé, F. Schrenk, C. Thèves, J.F. Tournepiche, S. Xing, and C.K. Sun. The Musée of Angoulême kindly granted access for study and scanning to the specimens from La Chaise Abri Bourgeois-Delaunay; for this, we are indebted to J.F. Tournepiche. These specimens were scanned by A. Mazurier and R Macchiarelli at the ESRF beamline ID 17 (Grenoble) in collaboration with A. Bravin, C. Nemoz and P. Tafforeau, and we thank them for sharing the microtomographic records. For technical assistance, we are indebted to F. Bernardini, F. de Beer, B. Duployer, J.W. Hoffman, and C. Tuniz. For scientific collaboration and discussion, we sincerely thank Y.M. Hou, M. Kundrát, R. Macchiarelli, and S. Xing. We thank Y.X. Hao for his insightful comments which helped us to improve the manuscript. We also acknowledge the Centre de Microtomographie at the Université de Poitiers, the European Synchrotron Radiation Facility (ESRF) beamline ID17 at Grenoble, the Neanderthal Studies Professional Online Service (NESPOS) Society and the Senckenberg Research Institute. This study granted access to the HPC resources of CALMIP supercomputing center under the allocation 2016-(P1440] attributed to the AMIS Laboratory, Toulouse. This work was supported by the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDA19050102 and XDB26000000), the National Natural Science Foundation of China (Grant No. 41702026 and 41872030), State Key Laboratory of Palaeobiology and Stratigraphy (Nanjing Institute of Geology and Palaeontology, CAS; grant number 173119) and The French Centre National de la

Recherche Scientifique (CNRS). We thank the Editor (David Alba), the Associate Editor, and three anonymous reviewers for their comments, which helped to improve the paper.


P Hyot

Benjamin Pochet

Grevin

Y J Dappe

Y Almadori

M T Dau

C Vergnaud

M Jamet

C Paillet

B Hyot

P Pochet

Benjamin Grévin
email: benjamin.grevin@cea.fr

Charge transfers and charged defects in WSe 2 /graphene-SiC interfaces

Keywords: non-contact atomic force microscopy, Kelvin probe force microscopy, transition metal dichalcogenide, van der Waals heterostructures, graphene, charge transfer, Schottky-Mott model, charge defect ORCID

Charge transfers and charged defects in WSe 2

/graphene-SiC interfaces

Introduction

Surface electrostatics represents a wide field of science due to the large diversity of potential applications. Hence, electrostatic fields and charge transfers at surfaces and interfaces have been widely studied for electronic properties modulation [1], nanostructuration [2], molecular structures on surfaces [3], or applications for nanoelectronic devices [4]. The richness of this area arises from the specific surface properties compared to the bulk aspects. On either metallic, semiconducting or insulating surfaces, these effects are of paramount importance to understand new fundamental physics and develop technological applications. In that respect, the recent discovery of a multitude of twodimensionals (2D) materials following the rise of graphene allows considering almost perfect systems for surface and interface electrostatics studies. Indeed, van der Waals heterostuctures that combine the properties of several constituents, associated to the constant increase of the amount Figure 1 Atomic representation of the different unit cells used for the DFT calculations: (a) SLG/ZLG/SiC, (b) BLG/ZLG/SiC, (c) WSe2/SLG/ZLG/SiC, (d) WSe2/BLG/SiC and (e) WSe2-180/WSe2/SLG/ZLG/SiC. In the last structure, the two WSe2 layers present a mutual angle of 180o in order to avoid any particular symmetry aspect. NB: in this work, the graphene buffer layer is referred to as Zero Layer Graphene (ZLG). SLG: single (or mono-) layer graphene. BLG: bi-layer graphene. The calculated amounts of transferred charges are given for each atomic layer, positive or negative signs indicate an excess or a deficit of electrons, respectively. of available 2D materials of different types and nature (metals, semiconductors and insulators) open fantastic perspectives for the exploration of novel 2D solid states physics [5]. This new fundamental research also paves the way toward novel low dimensional and high performance optoelectronic devices [6].

Among the multitude of potential vertical 2D heterostructures, coupling mono-or bi-layer graphene and transition metal dichalcogenides (TMDs) appears to be a promising approach to develop new optoelectronic systems, according to their complementary properties. As already well known, mono-layer graphene is an almost ideal semi-metal exhibiting extremely high carrier mobility and electrical conductivity thanks to the relativistic character of particles at the K-point of the reciprocal space [7]. On the other hand, group VI-B mono-layer transition metal dichalcogenides (TMDs) of the form MX2 (M = Mo, W ; X = S, Se, Te) are direct bandgap semiconductors which exhibit exceptionally strong light matter interactions regarding to their thickness [8].

Several works have demonstrated that photo-induced interlayer charge and/or energy transfer [9,10] can occur at graphene-TMD interfaces. Various kinds of hybrid graphene-TMD photodetectors have been successfully realized, featuring ultrahigh gains [11], broad-band photo-responsivity [START_REF] Mudd | High Broad-Band Photoresponsivity of Mechanically Formed InSe-Graphene van der Waals Heterostructures[END_REF], fast photo-responsivities [START_REF] Massicotte | Picosecond Photoresponse in van der Waals Heterostructures[END_REF], and efficient near-infrared photo-detection [START_REF] Yu | Near-Infrared Photodetectors Based on MoTe2/Graphene Heterostructure with High Responsivity and Flexibility[END_REF]. We refer the reader to review articles [START_REF] Geng | Heterostructures for High-Performance Photodetectors[END_REF] for a more comprehensive introduction to the field of hybrid graphene/TMD photo-detection. However, the interlayer charge/energy transfer processes remain the subject of intensive investigations [10,[START_REF] Yuan | Photocarrier Generation from Interlayer Charge-Transfer Transitions in WS2-Graphene[END_REF].

In that frame, it remains crucial to improve our understanding of the band alignment and charge exchange in 2D TMD/graphene heterojunctions [START_REF] Larentis | Band Offset and Negative Compressibility in Graphene-MoS2[END_REF][START_REF] Kim | Band Alignment in WSe2-Graphene Heterostructures[END_REF][START_REF] Jin | Tuning the Schottky Barrier at the Graphene/MoS2 Interface by Electron Doping: Density Functional Theory and Many-Body Calculations[END_REF][START_REF] Pierucci | Band Alignment and Minigaps in Monolayer MoS2-Graphene van der Waals Heterostructures[END_REF][START_REF] Quang | Scanning Tunneling Spectroscopy of van Der Waal Graphene/Semiconductor Interfaces: Absence of Fermi Level Pinning 2[END_REF][START_REF] Ulstrup | Nanoscale mapping of quasiparticle band Alignment[END_REF]. To address these questions, TMDs deposited on epitaxial graphene, itself grown on SiC(0001), constitute an excellent template. In the case of MoS2, it has notably been shown that the graphene-TMDC interlayer interaction goes beyond a simple van der Waals interaction, in part due to charge transfer effects [START_REF] Dau | Beyond van Der Waals Interaction: The Case of MoSe2 Epitaxially Grown on Few-Layer Graphene[END_REF]. By using SiC(0001), one can easily obtain a surface featuring terraces with single-layer graphene (SLG) or bilayer graphene (BLG). This allows investigating at once the influence of the electronic properties of the substrate on the band offset and charge transfer mechanisms. In the case of WSe2/Gr and WS2/Gr heterojunctions (HJs), scanning tunneling spectroscopy [START_REF] Quang | Scanning Tunneling Spectroscopy of van Der Waal Graphene/Semiconductor Interfaces: Absence of Fermi Level Pinning 2[END_REF] and angle resolved photoemission [START_REF] Ulstrup | Nanoscale mapping of quasiparticle band Alignment[END_REF] measurements have in particular shown that the band offsets track the work function of the underlying single layer and bilayer graphene.

These observations point towards the absence of Fermilevel pinning in these HJs [START_REF] Liu | Van der Waals metal-semiconductor junction: Weak Fermi level pinning enables effective tuning of Schottky barrier[END_REF], which band alignment shall therefore follow the prediction of the Schottky-Mott model. A key ingredient of the Schottky-Mott model is that there should be no electrostatic discontinuity at the metal-semi-conductor (MS) interface (or equivalently that the vacuum level alignment is kept on both interface sides), since there are no interface states.

In the present work, we probe surface and interface electrostatics of tungsten diselenide (WSe2) layers deposited by molecular beam epitaxy (MBE) on graphene mono-layers (SLG) and bi-layers (BLG), themselves grown on SiC(0001) substrates. Figure 1 provides an atomic representation of these interfaces. To carry out these investigations, we combine experimental Kelvin probe force microscopy (KPFM) measurements and Density Functional Theory (DFT) calculations, both techniques giving complementary access to the structural and electrical properties of the two-dimensional (2D) van der Waals heterojunctions (vHJ). KPFM is well suited to probe the band offsets at 2D HJs, charge transfer and/or interfacial dipoles. This technique yields a measurement of the vacuum-level variations at a sample surface, which are directly related to the energy band offsets [START_REF] Cahen | Electron Energetics at Surfaces and Interfaces: Concepts and Experiments[END_REF] taking place at junctions between two different materials.

The KPFM surface potential images and the DFT calculations reveal the existence of an interface dipole resulting from an electron transfer from the graphene to the TMD. We show why this electrostatic discontinuity at the WSe2-Gr interface is not in contradiction with the prediction of the Schottky-Mott model. In addition, we also observe a lateral band bending effect attributed to negatively charged defects at the TMD flakes edges, which amplitude is found to be consistent with the results of former STS reports [START_REF] Quang | Scanning Tunneling Spectroscopy of van Der Waal Graphene/Semiconductor Interfaces: Absence of Fermi Level Pinning 2[END_REF].

Methods

Samples synthesis and processing

1-2 layer-thick graphene substrates were directly grown by chemical vapor deposition (CVD) on top of silicon-terminated face (0001) of n-doped (Nitrogen-doped) 6H-SiC. The CVD growth was performed by propane decomposition (C3H8) at 1550°C, under a mixed atmosphere of argon (Ar) and dihydrogen (H2) at 800 mbar (1 H2 for 10 Ar ratio). The CVD growth under high partial pressure of Ar demonstrated high quality graphene formation [START_REF] Portail | Growth Mode and Electric Properties of Graphene and Graphitic Phase Grown by Argon-Propane Assisted CVD on 3C-SiC/Si and 6H-SiC[END_REF]. It has been shown that a low H2/Ar ratio allows forming the so-called buffer layer in between the SiC substrate and the graphene layer [START_REF] Michon | X-ray Diffraction and Raman Spectroscopy Study of Strain in Graphene Films Grown on 6H-SiC[END_REF]. Increasing the H2/Ar ratio results in in-plane rotational disorder, and to the formation of wrinkles [START_REF] Michon | X-ray Diffraction and Raman Spectroscopy Study of Strain in Graphene Films Grown on 6H-SiC[END_REF].

All TMD monolayers (1L-TMD), namely tungsten diselenide (WSe2), were grown in van der Waals epitaxy [START_REF] Walsh | van der Waals Epitaxy: 2D Materials and Topological Insulators[END_REF] on top of the 1-2 layer-thick graphene surface in a Molecular Beam Epitaxy (MBE) chamber, allowing the formation of vertical 2D TMD/Gr van der Waals heterojunctions (vHJ). The MBE growths were performed with a base pressure of about 5.10 -10 mbar which increased up to about 2.10 -8 mbar during the co-deposition of W and Se. The W and Se were respectively evaporated by an e-gun evaporator and a homebuilt evaporation cell. The Se pressure measured at the sample position with a retractable Bayard-Alpert gauge is about 1.10 -6 mbar. The W deposition rate measured with a quartz balance monitor was 1.5 Å.min -1 . The post-annealing was performed at 720°C for 15 minutes under Se flux. A Se capping of about 100 nm is deposited after TMDs synthesis in order to prevent sample degradation during the air exposure from the MBE chamber to the AFM UHV-set-up. A medium-temperature annealing at 450°C is then performed for 1 hour under UHV in the preparation chamber of the VT-AFM to remove the Secapping before nc-AFM/KPFM measurements.

nc-AFM/KPFM measurements nc-AFM/KPFM

experiments on graphene and WSe2/graphene were achieved under UHV conditions within a beam deflection Omicron VT-AFM setup. Topographic imaging was performed in frequency modulation mode with frequency shifts (negative set-points for nc-AFM) and vibration amplitude of a few hertz and a few tens of nanometers, respectively. KPFM measurements were simultaneously acquired in frequency modulation mode (FM-KPFM) in single-pass mode. All data were obtained with super-sharp silicon cantilever (SuperSharpSilicon, Nanoworld, resonance frequency range: 250-320 kHz) annealed and ion sputtered under UHV in the preparation chamber of the VT-AFM to remove contaminants. The lateral resolution is estimated at 3-5 nm. The modulation bias Vac (0.5 V peak-to-peak) and the KPFM compensation potential Vdc were applied to the cantilever (tip bias Vtip = Vdc), and the Gr substrate was grounded. As Vdc is applied to the tip, the contact potential difference (CPD) becomes equal to -Vdc. In this work, the tip compensation bias (Vtip = -VCPD) was used to represent the potentiometric data and was called for simplicity KPFM potential or surface potential.

Theoretical Calculations

Our calculations have been performed using the very efficient DFT localized orbital Fireball code [START_REF] Harris | Simplified Method for Calculating the Energy of Weakly Interacting Fragments[END_REF]. Since the main features have been widely described elsewhere [START_REF] Jelínek | Multicenter Approach to the Exchange-Correlation Interactions in Ab Initio Tight-Binding Methods[END_REF], we just mention here that this code uses a self-consistent version of the Harris-Foulke LDA functional [START_REF] Harris | Simplified Method for Calculating the Energy of Weakly Interacting Fragments[END_REF][START_REF] Foulkes | Tight-Binding Models and Density-Functional Theory[END_REF] and that the selfconsistency is achieved over the occupation numbers. Optimized numerical basis sets have been used for tungsten, selenium, carbon and silicon with respective cutoff radii in atomic units of s = 4.7, p = 5.2, d = 4.5 for W, s = 4.3, p = 4.8, d = 5.3 for Se, s = 4.5, p = 4.5 for C [START_REF] Basanta | Optimized Atomic-Like Orbitals for First-Principles Tight-Binding Molecular Dynamics[END_REF] and s = 4.8, p = 5.4 for Si. The Fireball code has been used to optimize the different structures considered in this work, namely, a 3 × 3 unit cell of graphene buffer layer (or Zero Layer Graphene, ZLG) on top of a five layer SiC slab, with either one or two graphene layer on top (as reference systems), and one or two WSe2 layer on top of all. These different systems, represented in Figure 1, have been optimized until the forces were below 0.1 eV/Å. Since vdW interaction is not properly described in DFT, we have considered an extra perturbation theory method, the LCAO-S2 + vdW approach [START_REF] Dappe | Intermolecular Interaction in Density Functional Theory: Application to Carbon Nanotubes and Fullerenes[END_REF], to accurately determine the interlayer equilibrium configuration. This approach is based on the dipolar approximation for vdW interaction and has already been used successfully in graphitic and 2D materials, in good agreement with experimental determinations [START_REF] Dau | Beyond van Der Waals Interaction: The Case of MoSe2 Epitaxially Grown on Few-Layer Graphene[END_REF][START_REF] Pierucci | Self-Organized Metal-Semiconductor Epitaxial Graphene Layer on Off-Axis 4H-SiC[END_REF][START_REF] Dappe | Effect of van der Waals Forces on the Stacking of Coronenes Encapsulated in a Single-Wall Carbon Nanotube and Many-Body Excitation[END_REF].

Once the equilibrium configuration is found, the electronic structure of each system has been analyzed by means of band structure and charge transfer calculations. A set of 300 specific k-points along the Γ-K-M path has been used for band structure calculations.

Results and discussion

At first, it is mandatory to characterize finely the CVD grown graphene substrate on SiC, since it forms template for the TMD growth. Figure 2 displays topographic and potentiometric nc-AFM/KPFM images recorded on the graphene substrate after in situ (i.e. under UHV) annealing at ca. 720°C (same temperature that the one used for the WSe2 growth process). The topographic image shown in Figure 2(a) features several hundred nanometers wide parallel terraces. The step height deduced from height profiles (not shown) is mainly of 7.5 ± 0.3 Å (three interatomic Si-C distance each of 2.5 Å), as expected for the surface morphology of few-layer graphene grown on 6H-SiC (0001) substrates. In addition, no height differences can be measured within terraces, indicating the growth of a continuous graphene layer with a homogeneous thickness over entire terraces. In particular, the absence of wrinkles, generally characteristic of in-plane rotational disorder [START_REF] Michon | Direct Growth of Few-Layer Graphene on 6H-SiC and 3C-SiC/Si via Propane Chemical Vapor Deposition[END_REF], suggests a monocrystalline character over several microns (scale of the nc-AFM images) or at least on a given terrace. topographic image. This surface potential map permits to identify regions with two distinct contrast levels: a homogeneous surface potential level observed over the majority of the sample and darker patches, displaying a lower surface potential compared to the main level. Consistently with the results of the literature [START_REF] Filleter | Local Work Function Measurements of Epitaxial Graphene[END_REF] and with complementary characterization performed by scanning electron microscopy (SEM) co-localized with Raman spectroscopy mapping (see Figure S1 in the supplementary information), the main level can be attributed to single-layer graphene (SLG), while dark patches corresponds to bilayer graphene (BLG) regions.

This relies on the fact that the graphene doping by the substrate decreases with the graphene thickness. In the case of graphene grown on SiC(0001), it is well known that the first graphene layer is n-doped by by the buffer layer with a (6√3 6√3)R30° periodicity between the graphene and the SiC [START_REF] Ristein | Origin of Doping in Quasi-Free-Standing Graphene on Silicon Carbide[END_REF][START_REF] Miwa | Waals Epitaxy of Two-Dimensional MoS2-Graphene Heterostructures in Ultrahigh Vacuum[END_REF]. Due to the electrostatic screening by the first layer, the second layer is less n-doped, resulting in a larger work function in case of BLG regarding to the SLG work function. The resulting band alignment is depicted in Figure 2(d), which shows the electron-doped Dirac cones for SLG and BLG (less n-doped and with a gap opening, as discussed hereafter) before and after contact at thermal equilibrium. Due to the BLG higher work function, the local vacuum level is shifted upward over the BLG. This shift is compensated by applying a more negative compensation bias on the AFM tip, accounting for the observed contrasts.

This situation is in overall well reproduced by the DFT calculations for mono-and bi-layer graphene on ZLG/SiC (here ZLG stands for Zero Layer Graphene, a 3 × 3 unit cell of graphene buffer layer, see Methods). In Figure 2 (e) and (f) we represent the band-structure calculations on these two systems note that all the energies are referred to the Fermi level). In Fig. 2(e), we can observe the graphene n-doping induced by the buffer/SiC substrate with a Dirac point located at 450 meV below the Fermi level. Fig. 2(f) shows a gap opening of ~160 meV in the graphene bilayer, due to the symmetry breaking induced by the charge transfer with the substrate. Hence, the lower graphene layer in the bilayer is doped by the ZLG/SiC substrate and presents a different charge density with respect to the upper graphene layer. As a consequence, an electric field is generated between the layers, which breaks the symmetry and opens a gap in the structure [START_REF] Dau | Beyond van Der Waals Interaction: The Case of MoSe2 Epitaxially Grown on Few-Layer Graphene[END_REF][START_REF] Pierucci | Self-Organized Metal-Semiconductor Epitaxial Graphene Layer on Off-Axis 4H-SiC[END_REF]. From the DFT calculations, we can extract the electronic charge transfers of each layer in the unit cell. The lower graphene layer gains 0.18 electrons per unit cell, coming from the ZLG/SiC substrate, whereas the upper graphene layer gains only 0.04 electrons coming from the lower graphene layer (see Figure 1(b)).

The KPFM surface potential variation between SLG and BLG (see the histogram in Figure2 (c)), permits to estimate a work function difference of 135 ±2 mV, in excellent qualitative agreement with previous studies performed on few-layers graphene grown on 6H-SiC (0001) [START_REF] Filleter | Local Work Function Measurements of Epitaxial Graphene[END_REF][START_REF] Lee | The Interlayer Screening Effect of Graphene Sheets Investigated by Kelvin Probe Force Microscopy[END_REF][START_REF] Hibino | Dependence of Electronic Properties of Epitaxial Few-Layer Graphene on the Number of Layers Investigated by Photoelectron Emission Microscopy[END_REF][START_REF] Ziegler | Variations in the Work Function of Doped Single-and Few-Layer Graphene Assessed by Kelvin Probe Force Microscopy and Density Functional Theory[END_REF][START_REF] Panchal | Standardization of Surface Potential Measurements of Graphene[END_REF][START_REF] Melios | Surface and Interface Structure of Quasi-Free Standing Graphene on SiC 2[END_REF][START_REF] Leenaerts | The Work Function of Few-Layer Graphene[END_REF]. As a matter of comparison, the Fermi level calculated for the SLG/ZLG/SiC system is shifted above the one of the BLG/ZLG/SiC by +102mV. This slight underestimation may originate from the intrinsic defects of DFT which does not properly take into account the excitations in a system, leading to an underestimation of the gap. In the present case, the gap underestimation is also related to a slight underestimation of the Fermi level position.

We turn now to the structural and electronic properties of the WSe2/Graphene heterojunction (HJ). The growth conditions were chosen in order to get a sub-monolayer TMD coverage and keep uncovered graphene regions, so as to reference the KPFM surface potential and to be able to discriminate between SLG and BLG. Figures 3(a,b) shows the nc-AFM topography and the associated surface potential (SP) of a representative area of the WSe2/Gr sample. From the topographic image, substrate terraces delimited by step edges can still be observed. The images clearly reveal the growth of WSe2 flakes covering approximately 30% of the surface (this surface coverage value has been obtained by analyzing the distribution of the SP levels, see Figure S2 in the supplementary information). Thanks to the reference measurements carried out on the bare substrate (Figure 2), distinct areas can be unambiguously associated with SLG or BLG (the location of a BLG patch is highlighted by an arrow in Figure 3(b)). Both KPFM measurements (i.e. on the bare substrate and on the WSe2/Graphene sample) yield the same work function difference between SLG and BLG within the error bar (compare the SP histograms in Figure 2(c) and in Figure 3(c)).

Whether we are considering SLG or BLG areas, the surface potential measured over the WSe2/Gr heterojunctions is in average lower by several hundreds of mV than the one measured over the bare graphene layer. This effect can be understood in two ways: (i) the lower SP measured over HJ indicates a larger WF (ii) the SP variation can be interpreted by the presence of an effective out-of-plane surface dipole pointing downward, resulting from an interfacial charge transfer between the WSe2 flakes and the graphene substrate. The downward orientation of the surface dipole indicates that electrons have been transferred from the underlying graphene to the WSe2 flakes.

Further insight on the nanoscale morphology and electronic properties can be gained by recording images at smaller scales. Figure 4a presents a 3D plot of the surface topography of two WSe2 islands, grown at the vicinity and over a graphene step. In Figure 4b, the same topographic levels are combined with the KPFM potential information (used for the color code) in a unique 3D plot. This mixed representation facilitates the identification of correlations between the topographic and potentiometric contrasts.

These images reveal that the WSe2 domains consist in mono-layer (1L) islands partially covered by a second layer (2L), with stacking heights of ca. 0.68nm (see the topographic profile in Figure 4c). These values are close to the ones reported by others for MBE-grown WSe2 on graphite [START_REF] Park | Scanning Tunneling Microscopy and Spectroscopy of Air Exposure Effects on Molecular Beam Epitaxy Grown WSe2 Monolayers and Bilayers[END_REF]. In addition, the topographic image reveals the existence of defects, which appear in the form local sur-elevations (typ. 0.3nm) with regard to the expected stacking heights for the first and second layers. In Figure 4a, two arrows labelled D-1L and D-2L pinpoint the location of some of these defects in the first (D-1L) and second (D-2L) WSe2 layers (additional topographic profiles crossing these corrugations are provided in the supplementary information, see Figure S3).

A closer examination reveals the existence of other kinds of defects. First, one can see that the domain edge is slightly sur-elevated in most of its parts with respect to the flake interior. It is thus likely that the edges are in part decorated with another material such as Se clusters or residual molecules (such as hydrocarbons, O2, H2O), which has already been observed for MBE-grown samples [START_REF] Park | Scanning Tunneling Microscopy and Spectroscopy of Air Exposure Effects on Molecular Beam Epitaxy Grown WSe2 Monolayers and Bilayers[END_REF][START_REF] Zhang | Direct Imaging of Band Profile in Single Layer MoS2 on Graphite: Quasiparticle Energy Gap, Metallic Edge States[END_REF][START_REF] Quang | Band-Bending Induced by Charged Defects and Edges of Atomically Thin Transition Metal Dichalcogenide Films 2[END_REF]. Similar surelevations forming inversion domain boundaries (DB) are also present inside the WSe2 flake (one DB is highlighted by an arrow in Figure 4). The existence of such domain (or twin) boundaries is also not surprising; others [START_REF] Dau | Beyond van Der Waals Interaction: The Case of MoSe2 Epitaxially Grown on Few-Layer Graphene[END_REF] have already mentioned their existence for TMD layers grown by molecular beam epitaxy [START_REF] Alvarez | Impact of a van der Waals interface on intrinsic and extrinsic defects in an MoSe2 monolayer[END_REF].

Obviously, the flakes are not homogeneously charged, and the nanostructure has a deep impact on the electrostatic landscape. First, one observes an important edge effect. The WSe2 islands are indeed clearly surrounded by dark halos in the surface potential images (Figure 4(b)). In other words, the surface potential falls down at the edge of the 1L-WSe2 flake, i.e. goes to more negative values with respect to the ones measured above the monolayer interior. A simple argument rules out the possibility of an artefact induced by tip convolution/averaging at the flake's edge: in such a scenario, we would expect the surface potential to take a value comprised between the ones measured above adjacent materials (i.e. between the interior of the 1L-WSe2 flake and SLG surface), which is not the case here. We also note that the surface potential reaches its lowest (i.e. more negative) level over the local structural defects that correspond to the highest deviations with respect to the expected height values for the first and second TMD layers (D-1L and D-2L); in other words these defects are negatively charged, too. However, the exact nature of these defects remains uncertain; we will thus restrict the discussion to the case of flake edges (and domain boundaries) in the following.

In first analysis, the observation of "dark SP halos" at the flake edges is in good agreement with the results of scanning tunneling spectroscopy measurements (STS) reported by LeQuang et al. for similar samples [START_REF] Quang | Band-Bending Induced by Charged Defects and Edges of Atomically Thin Transition Metal Dichalcogenide Films 2[END_REF]. Their STS experiments pointed out a significant lateral (in the 1L-TMD plan) band bending occurring in the vicinity of WSe2/vacuum interfaces. This effect was explained by assuming the existence of n-type-like charged defects at the flake edges, which electrostatic field is screened out by the free carriers in the graphene layer (as depicted in Figure 4d).

The presence of such negatively charged sites is consistent with the lower surface potential observed at the edge of the flakes, as well as over the domain boundaries. A band-bending of ca. 300mV (over 5nm) was extracted from the STS data [START_REF] Quang | Band-Bending Induced by Charged Defects and Edges of Atomically Thin Transition Metal Dichalcogenide Films 2[END_REF] on 1L-WSe2, which is close to the values deduced from the KPFM surface potential profiles (see the potentiometric profiles in Figure4c and in Figure S3). In our case, however, the lateral extension of the band bending region (of the order of a few nanometers) could not be accurately measured, due to convolution effects with the AFM tip (the impact of these effects on the topographic and potentiometric profiles is illustrated in Figure S3).

The band-bending value is obtained by comparing the surface potential levels measured at the edge and at the interior (over the first WSe2 layer) of the flake. Here, it is crucial to note that the surface potential at the interior of the flake (out of the band bending region) is lower by ca. 300mV than the one measured over the SLG. This confirms that the WSe2 islands are intrinsically negatively charged, or (equivalently) that an interface dipole oriented downward exist at the WSe2/graphene interface. It means that the negative charging of WSe2 domains does not solely originate from charged defects at the edges. The vertical interface dipole and the lateral band bending contribute simultaneously to the global electrostatic landscape probed by KPFM, as depicted in Figure 4d.

On a theoretical point of view, DFT calculations confirm the existence of an important electron transfer from the monolayer graphene to the WSe2 monolayer. In Figure 5 (left and middle pannels), we compare the DFT calculated bandstructures of SLG (SiC/ZLG/graphene) and of 1L-WSe2 on SLG (SiC/ZLG/graphene/ WSe2). In addition, the charges per unit cell resulting from the interfacial charge transfers are given in Figure 1. The WSe2 layer gains 0.95 electrons per unit cell from the monolayer graphene, which loses 0.75 electrons (per unit cell). Thus, the charge transfer cancels the n-doping from the ZLG/SiC substrate, as it can be seen from Fig. 5 with the relocation of the Dirac point closer to the Fermi level. Consequently, we can state that the TMD monolayer neutralizes the extra charge given to graphene by the ZLG/SiC substrate. We can interpret this charge transfer as an electrostatic reaction of the system to the initial nonequilibrium state of the graphene doped by the substrate. In that respect, the TMD is less perturbed by the gained electronic charges, due to its semiconducting nature, which allows absorbing a small electronic density in the gap, which is not the case of graphene of semi-metallic nature.

These joint experimental and theoretical evidences of a negative charge transfer from the graphene substrate to the TMD invite us to reconsider the band alignment model for WSe2/Gr-SiC interface. Several theoretical and experimental works have already shown that graphene-TMD interfaces display very weak Fermi level pinning thanks to the absence of metal-induced gap states (MIGS) [START_REF] Quang | Scanning Tunneling Spectroscopy of van Der Waal Graphene/Semiconductor Interfaces: Absence of Fermi Level Pinning 2[END_REF][START_REF] Liu | Van der Waals metal-semiconductor junction: Weak Fermi level pinning enables effective tuning of Schottky barrier[END_REF]. The conclusions of these studies are all in line with the Schottky-Mott model, which predicts that in the absence of interface states, the energy difference between the semiconductor's conduction band minimum and the metal work function at the interface (i.e. the Schottky barrier, SBH) shall be equal to the difference between the metal work function (M) and the semiconductor electron affinity (). As abovementioned, the Schottky-Mott model implies that there is no electrostatic discontinuity at the metal-semi-conductor (MS) interface, since there are no interface states. However, we note that in its classical form (i.e. in the case of an interface between three-dimensional materials), this model does not exclude the existence of an effective interface dipole resulting from the charge transfer (necessary to Fermi level alignment) across the MS interface. This charge transfer results in a band bending on the semiconductor side that goes with an equivalent shift of the vacuum level (this situation is illustrated in the literature [START_REF] Vilan | Chemical Modification of Semiconductor Surfaces for Molecular Electronics[END_REF]). In our case, there is no more space charge or band-bending region in the vertical direction due to the two dimensional nature of the stacked materials. The DFT calculations and KPFM data demonstrate nonetheless the existence of an interface dipole originating from electron transfer on both sides of the WSe2/Gr interface. Representing at once this dipole (with the vacuum level shift linked to it) and the energy bands in a Schottky-Mott configuration (fulfilling a vacuum alignment condition) is conceptually difficult (see Figure S4 in the supplementary information); this brings once again to light the exceptional nature of two-dimensional TMD/graphene interfaces.

Comparing the interface dipoles for WSe2/SGL and WSe2/BLG provides another way to approach the problem. As pointed out by LeQuang et al. [START_REF] Quang | Scanning Tunneling Spectroscopy of van Der Waal Graphene/Semiconductor Interfaces: Absence of Fermi Level Pinning 2[END_REF], the validity of the Schottky-Mott rule (i.e. a linear dependence of the Schottky barrier as a function of the metal work function) is preserved only if the interface dipole is identical for both systems. Such a comparison is presented in Figure 6, which displays nc-AFM/KPFM data recorded over two WSe2 domains grown on both sides of a step edge delimiting BLG -SLG areas. The overall contrast of the surface potential map (Figure 6b) reveals clearly that the vacuum level variations over the WSe2 flake follow the substrate ones. A quantitative analysis can be performed by drawing a cross section profile that matches approximately the flake interior on both sides of the SLG-BLG border; this analysis (Figure 6c) confirms that the vacuum level shift follows rigidly the one of the substrate within the measurement error bar. DFT calculations confirm that the charge transfer is almost identical across WSe2/SLG and WSe2/BLG interfaces, with an excess of 0.95 electrons >and 0.93 electrons per unit cell for the TMD layer (Figure 1), respectively. Overall, these results confirm that the band alignment and charge transfer for the WSe2/Gr system can be effectively accounted in the frame of the Schottky-Mott model.

To complete this study on the charge transfer, we finally consider the case of the second TMD layer. The DFT calculations (Figure 5) indicate that the first TMD layer absorbs almost all the excess charges contained in the graphene. Figure 1(e) presents the atomic configuration for two WSe2 layers on SLG/ZLG/SiC, rotated by 180 degrees one with respect to the other (this configuration was chosen in order to avoid any symmetry consideration, since no specific angle between the layers can be deduced from the experimental data). One observes no significant difference between the different band-structures calculated with one or two WSe2 monolayer (Figure 5). The charge transfer analysis (as represented in Figure 1(e)) reveals that the lower and upper WSe2 layers gain +0.89 and +0.07 electron per unit cell, respectively. This last result indicates that the second WSe2 layer (2L) is unlikely to play a significant role in the electrostatics of the system. Nevertheless, compared to the interior of the first WSe2 layer (1L), the surface potential is in average lower by 300mV over the second layer islands, as shown by the crosssection profile in Figure 4(c) (note that the structural defects displaying an improper stacking height have been excluded from this analysis). This apparent contradiction can be resolved by assuming that similarly to the case of 1L-WSe2, n-type-like charged defects exist at 2L-WSe2 edges. This hypothesis is also supported by the fact that the SP decrease over 2L domains is almost equal to the one observed over 1L edges. The absence of well-marked halos in the SP image at the periphery of 2L domains can be for its part understood by taking into account the reduced dimension of the 2L domains and convolution effects, which prevent an accurate visualization of the band bending region extension. Further investigations on samples displaying larger bilayer domains would be needed to definitely confirm that the slight negative charge excess calculated by DFT has a negligible impact on the electrostatic potential probed by KPFM.

Conclusion

In summary, we have investigated the electronic properties of MBE-grown WSe2 layers on graphene (single-layer and bilayer graphene on SiC(0001)) by Kelvin Probe Force Microscopy and Density Functional Theory (DFT) calculations. Overall, our results show that the global electrostatic landscape probed by KPFM originates from the joint contributions of charged defects and charge transfers at the TMD-substrate interface. Charged defects at the flake edges affect the electrostatic landscape probed by KPFM, by inducing a lateral band bending which value is fully consistent with the results of former investigations by scanning tunneling spectroscopy. The initial n-doping of the graphene monolayer on SiC is neutralized by the adsorption of the TMD, resulting in the existence of an interface dipole at the TMD/graphene interface. The vacuum level shift induced by this charge transfer is identical within the experimental error over singlelayer and bilayer graphene, which confirms that these interfaces can be described in a 2D version of the Schottky-Mott model where the spatially extended band bending area in the vertical direction on the semiconductor side is replaced by an interface dipole.

Figure 2 (

 2 Figure 2 (a) Topographic nc-AFM image of the graphene substrate on 6H-SiC(0001). (b) KPFM image of the tip compensation bias (Vtip =VKPFM, surface potential). (c) Histogram of the surface potential values recorded in the area delimited by dotted contours in (b). (d) Schematic illustration of the band alignment of Single-layer Graphene (SLG) and Bi-Layer Graphene (BLG) before contact (thought experiment) and after contact. The lower surface potential pointed out over BLG areas can be interpreted as an upshift of the local vacuum energy level (Evac), regarding to the one of SLG, as depicted by the after contact band alignment. The KPFM loop restores the vacuum level alignment by applying a compensation bias shifted in the opposite direction. At thermal equilibrium, the Fermi energy being equivalent for SLG and BLG, an upshift of Evac over BLG implies a larger work function (ΔΦ ≈ 135±2 meV). (e,f) DFT calculated bandstructures of graphene on ZLG/SiC (e) and bilayer graphene on ZLG/SiC (f). The dashed red, blue and green lines show respectively the Fermi level, the position of the Dirac cone of the monolayer graphene and the mid-gap of the bilayer graphene. All the energies are referred to the Fermi level

Figure 2 (

 2 b) displays the surface potential image (KPFM tip compensation bias) acquired simultaneously with the

Figure 3 (

 3 Figure 3 (a) Topographic nc-AFM image of the system after growth of WSe2 flakes on top of the graphene substrate. (b) KPFM image of the tip compensation bias (Vtip =VKPFM, surface potential). The arrow pinpoints a BLG area. (c) Histogram of the surface potential values recorded in the area delimited by dotted contours in (b).

Figure 4 (

 4 Figure 4 (a,b) Topographic nc-AFM image (a) and surface potential image (b) of two WSe2 islands located over SLG (255255nm). D-1L, D-2L and DB arrows pinpoint structural (and charged) defects in the first (D-1L) and second (D-2L) WSe2 layers, and a domain (or twin) boundary (DB). (c) Topographic and potentiometric cross-section profiles corresponding to the paths in (a) and (b). The vertical arrows highlight local topographic sur-elevations at the step edge and over the domain boundary. Half-transparent red bars (guideline for the eye only) indicate approximately the position of band-bending regions. (d) Schematic representation of the contributions to the surface potential of charged defects at the step edges and of the effective interface dipole (D) resulting from electron transfer from the graphene to the WSe2. The black and white dots represent a negatively charged defect at the WSe2 edge and an effective counter charge of opposite sign due to electrostatic screening by the graphene substrate.

Figure 5

 5 Figure 5 From left to right, DFT calculated band-structures of SLG (SiC/ZLG/graphene), 1L-WSe2 on SLG (SiC/ZLG/graphene/ WSe2) and of 2L-WSe2 on SLG (SiC/ZLG/graphene/WSe2-WSe2 180). The dashed red and blue lines show respectively the Fermi level and the position of the Dirac cone of the monolayer graphene.

Figure 6 (

 6 Figure 6 (a,b) Topographic nc-AFM image (a) and surface potential image (b) of two WSe2 islands grown on both sides of an SLG-BLG border (351351nm). (c) Potentiometric cross-section profile corresponding to the paths in (a) and (b). Red bars indicate the lateral positions in between which average values of the SP have been calculated. The SP shift between WSe2/SLG and WSe2/BLG has been calculated by considering the areas displaying the highest (in average) SP values. The error bars have been estimated from standard deviation values.

Acknowledgments

Financial support by the Agence Nationale de la Recherche (France) with the 2DJ project (ANR-15-CE24-0017) is gratefully acknowledged.

Supporting Information

Graphene characterization by co-localized SEM imaging and µ-Raman mapping

Raman measurements were performed using a confocal Renishaw InVia Raman microspectrometer using a 532 nm laser excitation source with an average power of 0.7 µW focused on a ~ 0.5 µm beam diameter. The top view SEM images were carried out with a scanning electron microscope (Zeiss FE-SEM Ultra) operating at 1.5 kV. Raman and SEM images are colocalized. Two marks were drawn on the sample with a diamond tip enabling the repositioning when using different characterization equipments.

Figure S1 presents the results of the (graphene) substrate characterization by co-localized scanning electron microscopy (SEM) and micro-Raman (µ-Raman) spectroscopy imaging.

Raman spectra and co-localized images obtained by superposing SEM images and Raman µmaps are shown in the upper panels and lower panels, respectively. SEM images (Figure S1(c,d)) feature dark patches covering approximately 5% of the surface (Area 2), while the rest of the sample's surface is mainly dominated by a brighter uniform contrast (Area 1). In addition bright strips, separated by several hundred of nanometer, can be observed and associated to SiC step edges. According to its sensitivity to work function differences [1][2][3], secondary electron contrasts reveal the presence of regions with different work functions, due to various numbers of graphene sheets, as demonstrated hereafter.

In order to ascertain the graphene thickness and get further information on its properties [4], co-localized Raman spectroscopy and mapping has been performed. In Figure S1(c),a hyperspectral mapping of the Full Width at Half Maximum (FWHM) of the graphene 2D band, named Γ2D in the following, is superimposed to the SEM image. The main information is that brighter regions, which represent higher values of Γ2D, show a one-by-one correspondence with dark SEM regions (Areas 2).

Figure S1(a) presents two finely located Raman spectra, respectively extracted from area 1 (black dot) and area 2 (red diamond). Note that the Raman contributions of the raw SiC substrate have been subtracted. Both spectra display several characteristic Raman features of graphene: (i) the G-band around 1600cm -1 , (ii) the 2D-band (also known as G'-band) close to 2700 cm -1 and (iii) a small defects-induced D-band around 1350 cm -1 , demonstrating the growth of graphene-like material over the whole sample. Moreover, the weak intensity of the D-band regarding to the G-band (ID/IG≈0.1), indicates a good crystallographic quality [4], while the other contributions close to the D-band (from 1200 cm -1 to 1600 cm -1 ) point out a 6√3 surface reconstruction [5,6], associated with the formation of an interfacial buffer layer at the graphene/SiC interface. In this case, it is worthy to note that, according to the literature [7,8], the few-layer graphene is expected to be n-doped by the SiC substrate through the buffer layer.

The number of graphene layers can also be inferred by analyzing the 2D-band shape which is a sensitive probe of the graphene thickness [4,9] (an analysis of the integrated area of the G-band [10] brings the same conclusions, see Figure S1(d)). Indeed, Γ2D increases as the number of layers increases, i.e. the thicker the graphene the higher the Γ2D value. In area 1 (black bottom spectra in Fig S1(a,b)), the 2D-band is well fitted by a single Lorentzian peak with Γ2D ≈ 35 cm -1 , highlighting the presence of single-layer graphene (SLG) [9] in these areas.

In area 2 (red upper spectra in Fig S1(a,b)), Γ2D is estimated at Γ2D ≈ 55 cm -1 and a good experimental data fitting is only performed thanks to 4 Lorentzian peaks, indicating Bi-Layer Graphene (BLG) stacked in AB-Bernal configuration [4]. as the one displayed in Figure 2). (b) Image obtained by keeping the areas for which the SP is lower than 75mV (the areas for which the SP is higher than this threshold are colored in red). This threshold has been fixed accordingly to the SP histogram values (see Figure 2 in the main text). The WSe2 flakes cover approximately 30% of the surface. there should be no electrostatic discontinuity (i.e. no shift of the vacuum level) at the graphene/TMD interfaces. In turn, the charge transfer across the interface (and the resulting dipole, D) induces a shift of the vacuum level. In the case of bulk materials, the equalization of the Fermi levels can be accounted by a band bending region in the vertical direction (an illustration can be found in [11]). In the case of 2D van der Waals interfaces, there is no band bending region in the vertical direction due to the two dimensional nature of the stacked materials. In other words, it is impossible to provide a correct representation of the interface that includes the Schottky Mott configuration at the interface and the vacuum level shift associated to the dipole. This diagram is thus provided only to stress that conceptual difficulty. SLG: single-layer graphene. BKG: bilayer graphene. SBH: Schottky barrier height. : work function. : electron affinity. D: interface dipole. Evac: vacuum level.

Supporting Information

Additional cross section profiles

Supporting Information