Atomic relaxations at the ( 0001 ) surface of Bi 2 Se 3 single crystals and ultrathin films

Sumalay Roy,1 H. L. Meyerheim,1,* K. Mohseni,1 A. Ernst,1,2 M. M. Otrokov,3,4 M. G. Vergniory,1,3 G. Mussler,5 J. Kampmeier,5 D. Grützmacher,5 C. Tusche,1 J. Schneider,6 E. V. Chulkov,3,4,7 and J. Kirschner1,8 1Max-Planck-Institut für Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany 2Wilhelm-Ostwald-Institut für Physikalische und Theoretische Chemie, Universität Leipzig, Linnéstraße 2, 04103 Leipzig, Germany 3Donostia International Physics Center (DIPC), 20018 San Sebastián/Donostia, Spain 4Tomsk State University, 634050 Tomsk, Russia 5Peter Grünberg Institute 9, Forschungzentrum Jülich, D-52425 Jülich, Germany and JARA, Fundamentals of Future Information Technologies 6Department für Geowissenschaften, Ludwig-Maximilians Universität München, D-80333 München, Germany 7Departamento de Fı́sica de Materiales UPV/EHU, Centro de Fı́sica de Materiales CFM–MPC and Centro Mixto CSIC-UPV/EHU, 20080 San Sebastián/Donostia, Spain 8Institut für Physik, Martin-Luther-Universität Halle-Wittenberg, D-06099 Halle, Germany (Received 30 July 2014; revised manuscript received 21 September 2014; published 30 October 2014)


I. INTRODUCTION
Bismuth selenide (Bi 2 Se 3 ) is being extensively investigated as an archetype three-dimensional (3D) topological insulator (TI) because of its large intrinsic band gap (0.3 eV), which makes it a primary candidate for potential spintronic device applications at room temperature [1][2][3].In this context it has been shown that the topological surface state is robust under the presence of surface structural disorder and upon deposition of foreign species as long as time reversal symmetry is conserved [4][5][6][7][8].While the investigation of the surface electronic structure of TIs has been the focus of experimental and theoretical studies, the atomic structure is much less investigated.Scanning tunneling microscopy (STM)-as useful as it is for analysis of the morphology and growth mode and, in special cases, for adsorption site determination [4,5,[9][10][11][12][13][14][15]-it does not allow either the precise determination of interatomic distances or the detection of subsurface atomic species and relaxations.
The first report of the tetradymite (Bi 2 Te 2 S)-type structure, to which Bi 2 Se 3 belongs, was made by Lange [16].Later, the bulk crystal structure of Bi 2 Se 3 was investigated by several authors [17][18][19].Bulk Bi 2 Se 3 has a rhombohedral crystal structure [space group D 5 3d (R32/m)], whose lattice parameters in the hexagonal setting are equal to a 0 = 4.14 Å and c 0 = 28.64Å.We note that the proper labeling of the basal surface is "(0001)" rather than (111) as often used; the latter is valid in the rhombohedral setting.Along the hexagonal c axis the crystal * hmeyerhm@mpi-halle.destructure is characterized by a sequence of Se-Bi-Se-Bi-Se quintuple layers (QLs).The hexagonal unit cell contains three QLs and hence 15 atomic layers.The QLs are connected by weak van der Waals (vdW) bonds, which make these crystals easily cleavable [20].It is generally assumed that cleaving leaves the crystal with a bulk truncated surface, i.e., with a topmost selenium layer, but a recent low-energy ion scattering study proposed a bismuth bilayer termination [21].This model has been challenged by a low-energy electron diffraction (LEED) analysis of the very similar phase Bi 2 Te 3 [22] and by a recent LEED and surface x-ray diffraction (SXRD) study of Bi 2 Se 3 [23,24] which found evidence of a single tellurium or selenium termination, respectively.
Considerable efforts have also been made to prepare ultrathin films of Bi 2 Se 3 on suitable substrates.One major motivation comes from the fact that, owing to the presence of selenium defects, single crystals are n-type semiconductors but not insulating in the bulk.In consequence, experiments to study the transport properties of the topologically protected metallic surface state are hampered [25][26][27].One possible solution comes from the possibility of preparing thin Bi 2 Se 3 films which are almost free of selenium vacancies.Epitaxially grown thin films are important for potential quantum device fabrications also.The preparation of (0001)-oriented Bi 2 Se 3 films on Si(111) substrates was investigated earlier [28][29][30] and studies mostly used transmission electron microscopy and electron diffraction experiments as a structure analysis tool.Recently, also LEED has been employed to determine the strain developed in thin films [31], which, in the case of Bi 2 Se 3 , was found to tune the Dirac state [32].
In order to elucidate the near-surface atomic structure we have carried out an in situ surface SXRD study of singlecrystal and thin-film samples.In addition, we have carried out an x-ray powder diffraction analysis to independently derive the structure parameters of bulk Bi 2 Se 3 , since the available data [17,19,33] show considerable deviations in detail.
The SXRD analysis provides clear evidence that the surface structure is characterized by an expansion of the first interlayer spacing ( d 12 /d 12 ), i.e., the vertical distance between the surface selenium layer and the second bismuth layer.This is unexpected since it is generally accepted that d 12 is contracted rather than expanded.We show that the expansion is related to trace amounts of impurities like carbon and oxygen near the surface, which may escape detection by conventional Auger electron spectroscopy (AES).An expansion in the range between 2% and 4% is observed.Contaminations which exceed the "trace" level are found to drastically increase d 12 /d 12 , to values up to 17%.

II. PREPARATION
Bulk Bi 2 Se 3 was grown by the vertical Bridgman method.Preparation of thin-film samples was carried out by molecular beam epitaxy (MBE) on Si(111) wafers.Prior to deposition, the Si substrates were chemically cleaned by the HF-last RCA procedure to remove the native oxide and to passivate the surface with hydrogen.The substrate was subsequently heated in situ to 600 • C for 20 min to desorb hydrogen.Deposition of bismuth and selenium was carried out using effusion cells operated at 470 • C (Bi) and 110 • C (Se).The shutter of the selenium cell was opened 2 s before that of the bismuth cell in order to saturate the silicon substrate surface with selenium.Throughout the growth, the substrate temperature was kept at 300 • C. Films of 13.5 nm thickness at a growth rate of 4.5 nm/h were prepared in order to obtain a smooth and uniform sample surface.
Subsequently, single crystals and MBE-grown samples were transferred into an ultrahigh-vacuum diffractometer equipped with standard surface analytical tools like LEED and AES.Surface cleaning was carried out by mild Ar + ion sputtering (E = 0.5 keV) followed by annealing.While single crystals can be sputtered for several hours and heated up to T = 500 • C, in the case of ultrathin films care has to be taken not to completely evaporate them.Using a type K thermocouple attached closely to the sample in combination with a calibrated infrared pyrometer (ε = 0.07), we found a limit of T ≈ 430 • C above which rapid evaporation of the Bi 2 Se 3 film sets in.Furthermore, sputtering with an approximate current density of j = 0.7 μA/cm 2 leads to the complete removal of the film within about 80-130 min of sputtering time.In summary, these limitations impose considerable restrictions on the extent of the surface cleaning process.
Preliminary characterization was carried out by LEED and AES. Figure 1 shows LEED patterns of the bulk [Fig.1(a)] and the MBE-grown sample [Fig.1(b)] recorded at an electron energy of approximately 40 eV.In both cases first-order diffraction spots are visible, but other significant differences exist.First, the LEED pattern of the bulk sample exhibits the 3m point-group symmetry in correspondence with the p3m1 plane-group symmetry of the surface structure.Sharp spots on a low background are observed.By contrast, the LEED pattern of the MBE sample has a 6mm symmetry; the spots are considerably broader and the background is enhanced, indicating a less ordered structure in terms of terrace size and defect density.The 6mm symmetry of the diffraction pattern is a consequence of a stacking fault(s) within the films equivalent to a 180 • rotation of the structure about the c axis.This leads to an overlap of the crystallographically inequivalent (10L) and (01L) rods.
Some representative AES spectra are shown to scale in Fig. 2 for MBE (M) and bulk (B) samples.The AES transitions are labeled.The most important ones are the Bi-NOO transitions at E = 96 and 101 eV and the Bi-NNO transitions at E = 249 and 268 eV.The latter appear as a characteristic double-peak with a "w"-like profile.The correct notation for the two groups of bismuth lines is N 6,7 O 4,5 O 4,5 and N 4 N 7 O 4,5 respectively.However, for simplicity they are denoted by NOO and NNO transitions, respectively.Details can be found in Ref. [34].Note that the C-KLL transition (E ≈ 272 eV) cannot be resolved from the Bi-NNO transition at E = 268 eV, but carbon can be detected by an increase in the second peak relative to the first one (E = 249 eV; see Fig. 2).Sample M4 can be considered as (quasi-) "clean" on the basis of AES.The Bi-NNO peak at 268 eV has only about 75% the amplitude of the first one at 249 eV, in correspondence with the reference spectrum in Ref. [34].By contrast, samples M2 and B1 are increasingly contaminated by carbon.In addition, M2 also shows a significant amount of oxygen contamination and can be seen as an example of a "dirty" sample.For interpretation of the AES spectra the relative sensitivity of the AES transitions (for 3 keV primary electron energy) is considered: the Bi-NNO transitions at E = 249 and 268 eV are the least sensitive transitions in the spectrum.They are about a factor of 8 less sensitive than the Bi-NOO transitions at E ≈ 101 eV and roughly a factor of 3 less sensitive than the C-KLL transition.Finally, the O-KLL transition is more sensitive than the C-KLL transition by a factor of 2.5.We summarize that-although carbon belongs to the group of elements which are not well detected by AES-the estimate on the basis of the relative peak amplitudes between Bi-NNO and the C-KLL peak allows some conclusions regarding the degree of surface contamination.We estimate that the "clean" sample M4 is contaminated by less than 0.2 monolayer (ML) of carbon.Here and in the following we refer to 1 ML as 1 ad-atom per surface unit cell; i.e., 1 ML = 6.8 × 10 14 atoms/cm 2 .Support for the estimate is provided by a comparison between experimental and calculated top layer expansion (see below).Subsequent characterization was carried out by x-ray diffraction.In the first step the bulk structure was analyzed by powder diffraction.

III. STRUCTURE OF BULK Bi 2 Se 3 BY X-RAY POWDER DIFFRACTION
In the first step a powder diffraction analysis of bulk Bi 2 Se 3 was carried out using Mo-K α radiation (λ = 0.7092 Å) in the Bragg-Brentano geometry.Data were fitted using the Rietveld analysis.A very good fit could be achieved; the results are listed in Table I.In space group R32/m (No. 166) bismuth and selenium atoms occupy Wyckoff positions 3a at (0, 0, 0) and 6c at (0, 0, z).In consequence, only the z position of the atoms in (6c) are free positional parameters.In addition, for each atom an isotropic Debye parameter (B = 8π 2 u 2 , with u 2 being the mean squared displacement amplitude) was refined.Based on the refined coordinates the bulk interlayer spacings were calculated.The structure model in Fig. 4    shorter than the inner ones (d 23 = d 34 = 1.940Å).We estimate the uncertainty to be equal to 0.01 Å.
Table II compares the results with those of different previous studies.A reasonably good agreement exists with respect to the most recent study by Vicente et al. [19], where the differences are mostly below 1%.However, there is a larger disagreement (mostly concerning d 12 and d 56 ) with respect to the older studies by Nakajima [17] and-more markedly-the electron diffraction study by Semiletov [33].Due to the relatively small differences as compared to the work of Vicente et al. [19], in the following we refer to their bulk values when relaxations of d i,i+1 are normalized to the respective bulk values.

IV. STRUCTURE OF Bi 2 Se 3 (0001)
SXRD experiments were carried out in situ using laboratory equipment at the Max-Planck-Institute in Halle and at beamline 25B of the European Synchrotron Radiation Facility (ESRF) in Grenoble, France.For the laboratory experiments a microfocus x-ray source (Cu-K α ) was used.Reflected intensities were collected by a two-dimensional (2D) pixel detector (Pilatus 100k) [35].Similarly, the experiments at the ESRF employed a Maxipix two-dimensional detector using an x-ray wavelength of 0.82 Å.In the following we focus on the laboratory experiments and discuss several particularities related to data collection.
Integrated intensities were collected along the integer-order crystal truncation rods (CTRs) under total reflection (incidence angle = 0.3 • ). Figure 5(a) qualitatively shows the intensity distribution along the (10L) and the (01L) CTRs, where bright (yellow) and dark (green) shading represent high and low intensity, respectively.The CTRs arise due to the truncation of the bulk crystal.This makes the third reflection index, L = q z /c * , a continuous parameter.Here, q z and c * represent the momentum transfer normal to the surface and reciprocal lattice unit, respectively [36][37][38][39].
In the hexagonal setting of the rhombohedral structure bulk Bragg reflections appear for the condition (−h + k + L) = 3n, with n ∈ Z. Data collection is carried out by rotating the sample about its normal.This is carried out step by step while the detector position is adjusted simultaneously.Figures 5(b) and 5(c) show snapshots of diffraction patterns collected for the single crystal and the film where the samples were aligned so as to scatter close to the bulk (015) Bragg reflection.The detector is aligned in such a way that the surface normal approximately corresponds to the horizontal axis (L = q z /c * ) as indicated by the labels.Note that the reciprocal coordinate is not linear along L. In both images, the exposure time is about 30 s per image.For bulk Bragg reflections up to several 10 7 counts are collected by integrating over the spot, decreasing to only a few tens in the case of weak reflections along the CTRs.The image of the bulk sample shows several peaks here originating from different grains.With increasing distance from the position L = 5 these reflections rapidly disappear due to the large misalignment of the corresponding rods relative to the surface normal.
Also, for the thin-film sample two reflections are observed.Here the reason is that the film has a twinned structure and the bulk-like Bragg reflections (104) and (015) of the two twin domains are detected.Although the crystal is aligned close to the (015) reflection of one domain, also the (104) reflection of the other is observed because of its finite width and high intensity.The diffraction pattern of the MBE-grown sample has a remarkably low background and is free of spurious spots of other grains.In this respect, the diffraction pattern of the film is clearly different from that of the single crystal, which exhibits a high background and several powder-diffraction-like rings.This is attributed to the presence of polycrystalline grains and structural defects.In summary, on an atomic scale the surface of the single crystal is well ordered and of high quality in terms of long-range order (see the LEED pattern in Fig. 1), while on a mesoscopic/macroscopic scale the film is clearly of superior quality.

A. Data analysis
Integrated intensities (I obs ) were derived by setting an appropriate window around the reflections as discussed in Ref. [35].Structure factor intensities |F obs | 2 ∝ I obs were derived from the integrated intensities (I obs ) by correcting for instrumental factors (see, e.g., Refs.[35,40,41]).Symbols in Fig. 6 represent the |F obs | 2 along four symmetrically independent CTRs.The total (tot) structure factor intensity (|F tot | 2 ) scattered by the sample is given by the interference sum of the bulk and surface contribution: Here, F B represents the structure factor of the bulk CTR at integer-order positions (h,k) in reciprocal space, while F S corresponds to the surface contribution.The interference of F B and F S allows the determination of the position of adlayer atoms relative to the (1 × 1) surface unit cell.This can formally be written as |F tot | 2 = |F B + θ j f j exp[i2π (hx j + ky j + Lz j )]| 2 , where the summation runs over all atoms (j) within the (1 × 1) surface cell.The parameters θ j and f j are the relative occupancy of a position by the atom (j) and its atomic scattering factor, respectively.
Data analysis was carried out by calculating the bulk CTR and placing three QLs as a surface contribution onto the bulk truncated Bi 2 Se 3 crystal.In principle, this allows the assumption of far-reaching relaxations, but it turned out that only relaxations within the first QL are significant.
The structure refinement was carried out by least-squares fitting of the calculated structure factor amplitude |F calc | to the |F obs | using the program Prometheus [42].The best fit is shown by solid lines, which follow the experimental data very well over four orders of magnitude.We note that reflections in the vicinity of the bulk Bragg peaks (±0.15 reciprocal lattice unit) were omitted for the analysis.This is motivated by the fact that, close to the bulk reflections, spots from differently oriented grains pass through the detector window, thereby leading to an artificial enhancement of the collected intensity.
Selenium and bismuth atoms occupy only high-symmetry sites in the plane group p3m1.For this reason, only the z parameters are free parameters (i.e., one per layer), in addition to the displacement factor T (q) representing thermal and static disorder.It is given by T (q) = e −B•q 2 /4 , where q is the magnitude of the scattering vector and B represents the Debye parameter as noted above.For details we refer to Refs.[43,44].In general, for each QL one overall Debye parameter was refined.
The fit quality is quantified by the unweighted residual (R U ) [45].For the single-crystal samples R U values in the range between 0.12 and 0.16 were obtained when based on structure factor amplitudes.These values can be considered as reasonably good, although for SXRD better fits are achieved in many cases.The worse-fit quality achieved for the singlecrystal samples here is attributed to systematic errors related to the mediocre surface quality, which is characterized by macroscopic roughness and surface defects, the latter giving rise to spurious spots and powder rings.On the other hand, the fit quality achieved for the thin-film samples, with their flat and homogeneous surfaces, is quite good as outlined in the following.
Figure 7 shows one data set including the same rods as shown for the bulk crystals.Here, the structure factor intensity (|F (hkL)| 2 ) is plotted versus q z .It is given by the incoherent summation over the contributions originating from the two twin domains, (1) and ( 2 , where F 1 (hkL) and F 2 (hkL) represent the corresponding structure factors and the relative abundance of twin domain (1), respectively.The analysis yields ≈ 0.50; i.e., both twin domains are present with equal probability.
Sample twinning directly manifests itself by the equivalence of the (10L) and the (01L) rod as well as by the simultaneous appearance of bulk-like Bragg peaks at L values (4, 5, 7, 8, 10, 11, . ..) characteristic for both rods [see Fig. 5(a)].For the thin-film samples the intensity dynamics along L is less than that observed for the single crystals, which is related to the finite film thickness probed by the x rays.This is due to the dramatic increase in (static) disorder with depth, making the fourth QL almost invisible to the x rays (see below).
In general, fits to the MBE sample data are quite good, and unweighted residuals in the range of R u ≈ 0.16 to 0.20 are obtained based on intensities.We emphasize that the comparison of the fit quality achieved for the single crystals and the thin-film samples on the basis of the residuals requires a caveat, since R U values expressed in terms of intensities are-as a rule of thumb-larger by a factor of approximately 2 than those based on amplitudes.In consequence, the fits to the thin-film intensities are about ∼20%-50% better than those for the single crystals.In this context one might argue that the goodness-of-fit (GOF) parameter is the more suitable parameter for quantifying the fit quality.The GOF is given by GOF = 1/(N − P ) • [(I obs − I calc ) 2 /σ 2 ] [46], where the difference between observed and calculated intensities is normalized to the uncertainties expressed by the standard deviation (σ ) and to (N − P ), i.e., the difference between the number of independent data points (N) and the number of parameters (P) which are varied.However, it should be noted that there are some issues, which are important and must critically be treated when the GOF is used in SXRD.
First, the standard deviations (σ ) need to be properly estimated, which can be done from the quadrature sum of the statistical and the systematic uncertainty (see, e.g., Ref. [47]).While the statistical uncertainty is usually derived from the integration of the reflection intensity, the systematic uncertainty requires the collection of symmetryequivalent reflections, which-owing to the inhomogeneity and roughness of the surface of the single crystals-were not always accessible (e.g., beam blocking) or did not lead to acceptable results (disagreement of more than ∼30%).Second, the straightforward application of the formula including the denominator (N − P ) is problematic when applied to the CTRs, since the reflections along a rod cannot be viewed as independent in general.Consequently, the parameter N is not simply equivalent to the total number of reflections collected.There is an ongoing discussion about how to treat this issue.One suggestion is to derive N along a rod by referencing the total q z range to the width of the bulk Bragg peaks or by performing a Fourier analysis of the intensity distribution where the number of necessary Fourier components is equal to N [48].Finally, CTRs are often dominated by weak reflections (with larger σ values) in the presence of only a few strong ones (with smaller σ 's) close to the bulk Bragg reflections, which might bias the fit.For these reasons the reference to R U is preferable in the context of the present CTR analysis.

B. Structure model
The refinement of the z parameters of the first five layers, i.e., the first QL of the differently contaminated samples B1-B3 and M1-M4, led to the results which are shown schematically in Fig. 8.The modification of the interlayer relaxations relative to the bulk ( d i,i+1 /d i,i+1 , i = 1, . . .,5) is plotted as bars, where we refer to the bulk values (d i,i+1 ) published by Vicente et al. [19].The most important result is that for all samples an expansion of the top interlayer distance, i.e., of the vertical distance between the surface selenium layer and the second bismuth layer, is observed.This is unexpected, since for bulk truncated crystals a contraction of d 12 is generally expected.This is a consequence of the reduced coordination number the surface atoms experience.Ab inito calculations have indicated that the expansion is related to carbon "contamination" in the near-surface region.The energetically most favorable carbon site was found to be the interstitial between the first and the second layer as shown by the smallest, dark sphere in Fig. 4. The comparatively small coverage of 0.083 ML induces an expansion of d 12 /d 12 = 6% [24], as qualitatively indicated by the arrow in Fig. 4.This suggests that even trace amounts of carbon, which likely escape detection by AES, significantly affect the surface geometric structure and, with it, the electronic structure.This conclusion is supported by our study of a number of differently contaminated samples.For instance, sample M4 can be considered clean based on AES (see Fig. 2), but we find an expansion of d 12 /d 12 = 4%, suggesting that some nondetectable carbon impurities are present.This applies to other samples as well, such as M3, B2, and B3, which do not exhibit any traces of carbon and oxygen in the AES spectra (not shown) and whose top interlayer spacings are expanded by 3%.
On the other hand, AES indicates significant amounts of carbon for the bulk sample B1 and even additional amounts of oxygen for the MBE sample M2.We estimate the concentration to lie in the range of 0.2 to 0.5 ML.In parallel with the increasing amount of carbon and oxygen, d 12 /d 12 increases to 11% for B1 and up to 17% for M2.We conclude that the degree of top interlayer expansion is quite sensitive to the concentration of (near-) surface adsorbed carbon and oxygen.One can assume a linear relation between the carbon contamination and the degree of expansion.Based on the theoretical prediction in Ref. [24], which relates d 12 /d 12 = +6% to a concentration of 0.083 ML, the experimentally derived expansion d 12 /d 12 = +17% corresponds to a carbon concentration of 0.23 ML.This amount of contamination is readily observable by AES (see the spectrum for sample M2 in Fig. 2).
In general, relaxations of the deeper layers are not significant within the experimental uncertainty, which we estimate to be ±3% points.Inspection of Fig. 8 shows that the majority of all d i,i+1 /d i,i+1 with i > 1 show a scatter of ±3% around the corresponding bulk value.We emphasize one exception, namely, d 45 /d 45 for the MBE samples (M1 to M4), for all of which an expansion in the range between +2% and +6% is found.We may speculate that this spacing, which corresponds to the Bi-Se spacing at the lower end of the first QL at the vdW gap, is expanded due to some contamination located within or near the vdW gap.Finally, the vdW gap itself shows no sign of appreciable expansion or contraction.In order to investigate the effect of foreign gaseous species on the near-surface structure of Bi 2 Se 3 , also several experiments were carried out where the ultrahigh-vacuum chamber was flooded with molecular hydrogen and xenon up to a pressure of 1 mbar, which, however, had no detectable effect on the intensity distribution along the CTRs.In summary, with regard to the surface atomic relaxations and its response to carbon and oxygen adsorption we find a remarkable similarity between bulk and MBE grown Bi 2 Se 3 .
By contrast to the atomic relaxations some differences exist between bulk and MBE samples with regard to the structural disorder within deeper-lying layers.Figure 9 shows the Debye parameter (B) versus depth on a log scale for samples M1 and B3 as being representative for an MBE and a bulk sample, respectively.The refinement was carried out by allowing only one overall B factor to vary within one QL.The variation of this overall B factor versus the position of the QL within the sample is shown in Fig. 9. Whereas for the bulk sample B decreases from values in the range of 6 Å2 for the topmost QL (QL = 1) to B ≈ 2.5 Å2 for the fifth QL, it increases drastically, from B ≈ 10 Å2 for the surface QL to giant values, in the range of 40 Å2 , for the 4th QL in the case of the MBE-grown sample.
First, these numbers should be compared to the (thermal) disorder in bulk Bi 2 Se 3 , for which a maximum Debye parameter of B = 1.4 Å2 at room temperature is derived from the x-ray powder diffraction analysis.We recall that B = 1.4 Å2 is equivalent to the root mean square displacement amplitude of u 2 = (1/π ) √ B/8 = 0.13 Å, which corresponds to the half-width of the (Gaussian) probability density distribution function representing the displacement function of the atom in direct space.In general, the probability density distribution function is the Fourier transform of T (q) [43].Concerning the relatively small value of u 2 , we tentatively assume that it primarily represents thermal disorder, although thermal and static disorder cannot be separated without temperaturedependent measurements.
The disorder at the surface of the bulk single-crystal samples is characterized by B ≈ 6 Å2 , corresponding to u 2 = 0.27 Å, which is roughly a factor of 2 larger than the (maximum) bulk value.This is not unexpected and previous studies have found a similar behavior for high-index crystal surfaces such as Cu(001) [49].However, while for the latter the decay to the bulk value is rather fast (≈2 atomic layers), in the case of Bi 2 Se 3 single crystals the decay is very slow and B only decreases to ≈3 Å2 for the fifth QL, still significantly larger than the bulk value.This behavior is unusual and no clear-cut explanation can be given at present.We only note that bulk Bi 2 Se 3 is a relatively "soft" material (Debye temperature θ D ≈ 200 K [50]), and to some extent, static disorder may play some role.
The situation is dramatically different for the MBE-grown thin-film samples.In those cases, values of B = 20 and 40 Å2 correspond to u 2 = 0.50 Å and u 2 = 0.71 Å, respectively, which is attributed to static rather than thermal disorder.Its increase with depth suggests that its origin lies in the epitaxial growth of the MBE layers on the Si(111) substrate crystal.

V. THE ELECTRONIC STRUCTURE OF CARBON-DOPED Bi 2 Se 3 (0001)
In the following the effect of the top interlayer expansion and the surface carbon doping on the Bi 2 Se 3 (0001) electronic structure is theoretically investigated.First-principles calculations were used, which are known to yield a surface electronic structure in good agreement with experiment [51,52].The calculations were performed using the projector augmentedwave method [53] in the VASP implementation [54,55] and the generalized gradient approximation to the exchangecorrelation potential [56].The Hamiltonian contained the scalar relativistic corrections and the spin-orbit coupling was taken into account by the second variation method [57].We set the energy cutoff for the plane-wave expansion of wave functions to be equal to 400 eV and chose a -centered special k-point grid of 5 × 5 × 1 to sample the two-dimensional Brillouin zone.The Bi 2 Se 3 (0001) surface was simulated by a 6-QL-thick slab and (2 × 2) hexagonal supercell with 4 atoms per single layer.The latter enables us to simulate an adsorbate coverage of 0.25 ML.In the case of the carbon-contaminated surface, the carbon atom was placed in the interstitial site between the first (selenium) and the second (bismuth) layer, which was found to be energetically preferable in our previous work (see Ref. [24]).All calculations were performed using a model of repeating slabs separated by a vacuum gap of a minimum of 10 Å.
Figure 10 compares the calculated band structures of the pristine and perturbed (the latter with and without consideration of carbon doping) Bi 2 Se 3 (0001) along the K--M direction.As a perturbation, only the top interlayer expansion of d 12 /d 12 = 15% was considered since the relaxations of the deeper layers were found to be insignificant according to the SXRD analysis.The left panel in Fig. 10 shows the band structure of the pristine (unperturbed) Bi 2 Se 3 (0001).It is characterized by a linearly dispersing surface state with the Dirac point (DP) located at the Fermi level, in good agreement with previous studies [1,58,59].A 15% top interlayer expansion shifts the DP towards the conduction band by 138 meV without leading to additional states in the band gap (see middle panel in Fig. 10).Accounting for carbon doping (right panel in Fig. 10) results in electron doping, whereupon the Fermi level shifts into the conduction band.Simultaneously, the DP gets shifted towards the valence band (cf.middle panel in Fig. 10).Thus, as a result of the competition between these two factors, i.e., the upward DP shift due to the top interlayer expansion and the impurity-induced downward shift, the DP turns out to be moved upwards by 6 meV compared to the pristine surface.This is in qualitative agreement with a recent photoemission experiment [24].Note that in Ref. [24] the experimentally observed DP shift was somewhat larger, ∼40 ± 20 meV, for the top interlayer expansion d 12 /d 12 10%.We believe that this difference might be related to the higher concentration of carbon (0.25 ML) in the present calculations, which requires a larger expansion of the first interlayer spacing to induce the stronger upward shift of the DP.Finally, we note that the carbon p states are located above the Bi 2 Se 3 bulk conduction band minimum and therefore do not significantly affect the dispersion of the topological surface state inside the fundamental band gap.

VI. SUMMARY
In summary, we have investigated the near-surface structure of Bi 2 Se 3 (0001) using x-ray diffraction.Differently prepared samples were investigated: bulk (Bridgman method) and 13-nm-thick MBE-grown films deposited on Si(111).For both type of samples, a significant expansion of the top interlayer spacing was determined, which amounts to a maximum value of up to 17%.Comparison with Auger electron spectra collected from the same samples indicates a direct correlation of the expansion with the amount of near-surface contamination by carbon and oxygen.Even in the case where no contaminations could be detected by AES, a slight expansion of about 3% exists.Deeper layers as well as the vdW gap are not relaxed.Although bulk and MBE samples are similar with regard to the surface layer relaxations, they differ significantly with regard to structural disorder as measured by the Debye parameter B. According to first-principles calculations, the surface contamination by carbon leads to an upward shift of the DP inside the fundamental band gap but does not significantly modify the dispersion of the topological surface state.
FIG. 2. (Color online) DifferentialAuger electron spectra (primary electron energy, 3 keV) for a bulk (B) and two MBE (M)-grown Bi 2 Se 3 samples.Several AES transitions are indicated.Note that one Bi-NNO transition at E ≈ 268 eV almost coincides with the C-KLL line (E ≈ 272 eV).While sample M4 can be considered "clean," samples M2 and B1 can be viewed as being heavily contaminated with carbon and oxygen.

Figure 3
shows observed (I obs ) and calculated (I calc ) intensities.Vertical (green) bars in the lower part of the figure indicate the peak positions of the silicon standard (upper panel) and of Bi 2 Se 3 (lower panel), respectively.The difference I obs − I calc is shown by the horizontal (blue) line at the bottom.

FIG. 3 .
FIG. 3. (Color online) Powder diffraction pattern for bulk Bi 2 Se 3 .Circles (red) and black line correspond to the measured and fitted intensity, respectively.The difference is shown by the bottom horizontal (blue) line on the same scale.Reflection positions are indicated by the vertical (green) bars for the silicon standard (upper panel) and Bi 2 Se 3 (lower panel).

FIG. 4
FIG. 4. (Color online) Ball-and-stick model of Bi 2 Se 3 .Small (red) and large (blue) balls represent selenium and bismuth atoms, respectively.The interstitial carbon atom is shown by the smallest, dark ball.Interlayer distances are labeled d i,i+1 (see TableII).The arrow emphasizes the expansion of the top layer spacing.
FIG. 4. (Color online) Ball-and-stick model of Bi 2 Se 3 .Small (red) and large (blue) balls represent selenium and bismuth atoms, respectively.The interstitial carbon atom is shown by the smallest, dark ball.Interlayer distances are labeled d i,i+1 (see TableII).The arrow emphasizes the expansion of the top layer spacing.

FIG. 5 .
FIG. 5. (Color online) (a) Schematic of the intensity distribution along the (10L) and the (01L) rod of Bi 2 Se 3 (0001).Bright (yellow) and dark (green) color represents high and low intensity, respectively.(b), (c) Snapshots of detector images for a single crystal (b) and a thin-film sample (c).

FIG. 7 .
FIG.7.Experimental (symbols) and calculated (lines) structure factor intensities for the MBE-grown film.For details, see text.

FIG. 9 .
FIG. 9. (Color online) Debye parameter versus depth expressed by the number of the quintuple layer (QL) for samples M1 and B3.QL = 1 indicates the QL at the sample surface.Uncertainties are <10%.Note that B is shown on a log scale.

TABLE II .
Comparison of interlayer spacings (in Å) for bulk Bi 2 Se 3 derived from different studies.Parameters d i,i+1 refer to Fig. 4.