Surface Structure of Bi(111) from Helium Atom Scattering Measurements. Inelastic Close-Coupling Formalism

Elastic and inelastic close-coupling (CC) calculations have been used to extract information about the corrugation amplitude and the surface vibrational atomic displacement by fitting to several experimental diffraction patterns. To model the three-dimensional interaction between the He atom and the Bi(111) surface under investigation, a corrugated Morse potential has been assumed. Two different types of calculations are used to obtain theoretical diffraction intensities at three surface temperatures along the two symmetry directions. Type one consists of solving the elastic CC (eCC) and attenuating the corresponding diffraction intensities by a global Debye–Waller (DW) factor. The second one, within a unitary theory, is derived from merely solving the inelastic CC (iCC) equations, where no DW factor is necessary to include. While both methods arrive at similar predictions for the peak-to-peak corrugation value, the variance of the value obtained by the iCC method is much better. Furthermore, the more extensive calculation is better suited to model the temperature induced signal asymmetries and renders the inclusion for a second Debye temperature for the diffraction peaks futile.


■ INTRODUCTION
The electronic density structure of a surface determines its chemical behavior. While on surfaces like platinum, which is widely used as a catalyst, the effects of crystal face, surface steps, and kinks are well-known, more complicated electronic surface structures still lack a detailed treatment. Recently, the (111) surfaces of the semimetals bismuth (Bi) and antimony (Sb) have raised a lot of interest. Not only do they represent the two main ingredients of topological insulators, 1,2 but they also both present a fairly strong electronic surface density corrugation despite exhibiting conducting surface states. 3,4 The temperature dependence of these peculiar electronic structures may change the binding character of adsorbed species remarkably; 5 thus, it is essential to determine a complete picture of an electronic surface structure before conducting adsorbate experiments on them. Helium atom scattering (HAS) experiments provide a low-energetic, completely nondestructive means of investigation to measure the pure surface properties of materials. The inert neutral helium atoms are already repelled from the electronic density corrugation above the surface, probing only surface effects. Close-coupling (CC) calculations 6 provide a significant improvement compared to oversimple approximate methods. While the essential accurate knowledge of the interaction potential requires numerous measurements and a careful analysis, the effort may be worthwhile because the quantum mechanical treatment of the scattering procedure provides by far better insight into the scattering processes.
Earlier CC investigations included the Debye−Waller (DW) factor to account for the thermal attenuation of scattering intensities. Heavy materials like bismuth, however, exhibit very low surface Debye temperatures, indicating an extremely fast decay of scattering intensities with the surface temperature. The inclusion of inelastic channels into the standard CC equations 7 provides a natural extension within this theoretical framework, and as a result, the attenuation of the diffraction intensities arises automatically. When the number of diffraction channels (elastic and inelastic) is increased, the number of exiting channels increases and the initial flux of He atoms has to be redistributed among them, the elastic peaks being reduced in intensity. This fact implies that we do not have to add any ad hoc global factor, as the DW factor, to account for such an attenuation. Even more, within the inelastic CC (iCC) formalism, the theory is unitary. The unitarity in the elastic CC (eCC) calculations is lost when the DW factor is used to attenuate the diffraction intensities. In this work, we are going to use the iCC equations to extract information about the interaction potential. In general, the number of fitting parameters for this type of scattering is quite high and settles around five in this specific case. For surface temperatures around 300 K, the number of total channels playing a role in this scattering is also quite high, typically more than 200. On the other hand, when working on a fixed geometry experimental setup, the whole experimental diffraction pattern is not fitted at once by solving a complete set of iCC equations (up to numerical convergence). Each experimental diffraction peak of the corresponding pattern has to be fitted by solving a complete set of those equations. As a result, the fitting procedure becomes rather cumbersome. To circumvent the expensive high-dimensional fitting procedure, previously fitted eCC potential parameters have been assumed as given and fixed.
■ BISMUTH SURFACE STRUCTURE AND EXPERIMENTAL SETUP Bulk bismuth, like all of the heavier pnictogens, crystallizes in the rhombohedral A7 structure with two atoms per unit cell (space group R3m). A typical structural property of this crystal structure is the existence of puckered bilayers of atoms perpendicular to the [111] direction, as illustrated in Figure 1a.
The bonding of the atoms within these bilayers is of a strong covalent type, while the interbilayer bonding is closer to a van der Waals type. This is reflected in the relative distances between the layers, as labeled in Figure 1a. Because of this strong contrast in binding energies, bismuth crystals preferably cleave perpendicular to the [111] direction. The topmost layer of this prepared crystal reveals a 6-fold symmetry, as illustrated in Figure 1b. Despite the crystal's 3-fold symmetry, the Bi(111) surface can be treated as being 6-fold symmetric in low energy HAS experiments. 9,10 The lattice constant of this hexagonal-like surface structure has been determined to be a = 4.538 Å by LEED and HAS measurements. 3,8 This hexagonal surface structure leads to two distinguishable high symmetry directions in reciprocal space that are commonly denoted as ΓM and ΓK, as illustrated in Figure 1c. The reciprocal directions in Figure 1c correspond with the real-space directions in Figure 1b.
All measurements mentioned above and used in this work have been carried out on a helium atom scattering apparatus with a fixed source-target-detector angle of 91.5°that has been described in a previous publication. 11 The helium-atom beam is produced via supersonic expansion of He-gas through a cooled nozzle at 50 bar which is followed by a skimmer creating a spatially and energetically narrow beam (ΔE/E ≈ 2%). The bismuth sample is positioned on a 7-axis manipulator in the main chamber at a base pressure of 10 −11 mbar. The sample can be cooled using LN 2 or heated using a button heater while the temperature is measured by a type K thermocouple. The scattered He atoms are then detected by a quadrupole mass spectrometer followed by a multichannel analyzer. Angular elastic scans can be carried out by rotating the manipulator leading to diffraction spectra (0.1°resolution). Time-of-flight measurements allow us to record inelastic scattering spectra and are realized using a pseudo-random chopper disk with subsequent deconvolution of the measured signal. The diskshaped Bi(111) single crystal sample with a diameter of 15 mm and a thickness of 2 mm has been cleaned using several cycles of Ar + sputtering followed by annealing at 150°C. Surface cleanliness and contamination were checked via Auger electron spectroscopy (AES) and the intensity of the diffuse elastic peak, its orientation has been aligned using a low energy electron diffraction (LEED) system. Experiments can be carried out within a beam energy range of 15−25 meV with the sample cooled (113 K) or at room temperature in the two main symmetry directions of the crystal surface.

■ THEORETICAL BACKGROUND
Inelastic Close-Coupling Equations. The CC formalism provides a method for calculating the diffraction intensities of scattering experiments exactly (up to numerical convergence) in the elastic as well as in the inelastic regime. 6 The helium atom is considered to be a structureless and nonpenetrating particle, the one-phonon approximation is assumed, and the surface corrugation is described by a static as well as a dynamic time-dependent contribution. The time-dependent Schrodinger equation for a structureless particle is written as where squared wave vector quantities are given in energy units with ℏ 2 /2m = 1, with m being the mass of the incident particle. The standard notation is also used here where capital letters are for vectors parallel to the surface (2D) and small letters are for vectors in 3D. The gas−surface interaction potential, V, turns out to be dependent on time through the instantaneous position of the surface atoms, R + u (R,t), with u(R,t) being the deviation or displacement from the equilibrium position. If this displacement is considered to behave as a Gaussian function within the unit cell, it can be written along each surface degree of freedom as 7 with u z,0 the initial amplitude, ω the frequency of the active phonon mode, and σ c the parameter describing the width of the Gaussian function. As long as the relative displacements u are small compared to the lattice constant, the interaction potential can be Taylor expanded up to first order (within the so-called single-phonon approximation) as 12 From the layer description of lattice dynamics, it is well-known that the u displacement can be, in general, written as The Journal of Physical Chemistry C Article where the amplitude A includes the phonon polarization vector and the dependence on the surface temperature and ω ν (Q) is the frequency of the surface mode with quantum numbers (Q,ν). For most practical purposes, only displacements of atoms on the first layer significantly contribute to the interaction potential.
In the Taylor expansion given by eq 3, the zero order or static part of the interaction, V(r), is evaluated at zero displacements. Considering the periodicity of the lattice surface, this function can then be expanded into a Fourier series as with G being the 2D reciprocal lattice vector.
On the other hand, the wave function Ψ(r,t) has to take into account the double periodicity given by the Hamiltonian, in space and time. Thus, according to the Bloch theorem, Ψ(r,t) can be expanded as where n Q,ν stands for the number of phonons of the mode (Q,ν). In this work, we exclusively consider the inelastic effects on the elastic intensities. As commonly known, the Bragg law (for Q = 0) is written as Moreover, we also assume that only one mode is active in the scattering process and the coupling among phonons is neglected within the harmonic approximation. Thus, we can drop the subindex (Q,ν) in n Q,ν for the number of phonons, writing only n. Similarly, for the frequency of the active mode, we can simply write ω. After substituting eqs 3, 5, and 6 into eq 1, multiplying the resulting expression by exp[−i(K i + G)·R] and exp[−inωt], and then integrating over both time and the area of a single unit cell, one obtains 7 the following set of coupled differential equations for the diffracted waves ∑ ∑ and ∑ ∑ where is the z component of the kinetic energy for the (n,G)diffracted wave and is the contribution of the gradient of the interaction potential or vector force field (V′ represents the first derivative with respect to z); the first term represents the (x,y) components of the force, and the second one is its z component. The iCC equations are solved numerically by imposing the standard boundary conditions given elsewhere. 6 The theory is unitary; that is, the sum of diffraction probabilities (forming the diffraction pattern), for a given incident energy and angle, is equal to 1.
Inelastic and Elastic Channels: Floquet Blocks. Within this theoretical scheme, each diffraction channel is then represented by an effective potential formed by V 0 (z) plus the asymptotic energy, given by and is called an inelastic diffraction channel. Thus, any inelastic event, annihilation or creation, is represented by the transition from the entrance (or specular) channel to one of the channels with (n − 1) or (n + 1). Similarly, the wave functions associated with the discrete spectrum (bound states, labeled by v) of each channel are denoted by |K i + G, n, v⟩ and those associated with the continuum one (diffracted beams) by |K i + G, n, k G,n,z 2 ⟩. In the literature, it is also said that the inelastic channels are dressed by the phonon field. The number of channels dressed by a given number of phonons form a block, called a Floquet block. Thus, if only single-phonon scattering is considered, at least three Floquet blocks must be included in the calculation: the blocks dressed by minus and plus one phonon of the active mode and the block dressed by zero phonons or pure elastic channels (those used for an eCC calculation). The number of diffraction channels within a given Floquet block is formed at least by those used to obtain numerical convergence in an elastic CC calculation. Multiphonon contributions of the same active mode are taken into account by including more Floquet blocksthose dressed by two, three, or more phonons by following the staircase structure of eqs 8 and 9 through n ± 1.
Intrablock and Interblock Couplings. Furthermore, two coupling terms of very different nature are now present: V G−G′ (z) is responsible for the intrablock coupling, and the scalar function for the interblock one. The latter is responsible for the thermal attenuation of the diffraction intensities (see the second term on the right-hand part of eqs 8 and 9) described many times from a phenomenological viewpoint by a DW factor. In previous publications on the phonon dispersion of the Bi(111) surface, 9,10 the lowest lying, isolated Rayleigh mode was identified as the shear-vertical mode corresponding to the sole phonon dispersion line in the Debye model. According to the shear-vertical polarization of the suggested mode, the horizontal displacement of the lattice atoms can be neglected, simplifying the force term eq 11 to the vertical term and consequently the inelastic coupling term to The average thermal displacement A z (T) is related to the effective mean square displacement and has been estimated 7 to be  (16) with T the actual surface temperature, Θ D the surface Debye temperature, M the mass of a surface particle, a the lattice constant, k B the Boltzmann constant, and Q c is a fitting parameter for the Gaussian cutoff factor given by Q c = 2/σ c with σ c being the width parameter introduced in eq 2.
Averaging over Phonon Frequencies. On the other hand, when solving the iCC equations, frequency-dependent diffraction intensities are obtained, and these have to be averaged by assuming a density of phonons in order to compare with the experimental ones. The corresponding integration over phonon frequencies can be weighted by the Debye spectral density given by with ω D the Debye frequency. Thus, the final diffraction intensities are due to virtual phonon events only; no real phonon events are taken into account since the corresponding momenta are not involved in Bragg's law. The term "virtual events" denotes that when the phonon is created in the dynamics it has to be annihilated in order to have a net energy balance equal to zero. The origin of the attenuation in the iCC formalism is precisely due to these virtual phonon events since they are responsible for the appearance of the new inelastic channels. Note that the quadratic dependence on phonon frequency is strictly speaking only valid for the bulk; the surface would be better represented by a linear dependence. However, the quadratic term was chosen since the Debye model that the simulation is compared to is built upon the bulk description of the material. A simulation comparing the intensities using a linear and a quadratic term for the Debye spectral density yielded no mentionable difference for the relative elastic diffraction intensities. Debye−Waller Attenuation Factor. As previously mentioned, an alternative way to obtain diffraction intensities from the eCC equations that can be compared with the temperature-dependent experimental results and iCC calculations is by including a global attenuating factor, the DW factor. 13 By doing this, the unitarity of the attenuated eCC intensities is lost. This is an important theoretical inconsistency of this procedure. As observed in earlier measurements, 3 diffraction peak intensities are surface-temperature dependent. As known, the DW factor relates the intensity I(T S ) of diffraction peaks at temperature T S to the eCC intensity I 0 at zero surface temperature by means of where exp(−2W(T S )) represents the DW factor. Although the theoretical basis for the DW factor has been developed for neutron and X-ray diffraction, 14 a reasonable approximation for surfaces can be given by assuming zero parallel momentum transfer to the surface and final angles near the specular angle. In this expression, ⟨u z 2 ⟩ describes the average squared displacement of the atom perpendicular to the surface and Δk z is the momentum transfer perpendicular to the surface during the scattering event.
Assuming a harmonic oscillator within the Debye model, W(T S ) becomes 14,15 The applicability of the conventional DW factor which was introduced for X-ray diffraction in atom-surface scattering model, 16,17 has been discussed extensively in theory as well as in experiments. Different models have been discussed by Levi,18 who, for example, predicts an increase of diffraction intensities for soft potentials. Deviations from the predicted temperature dependencies of DW factors especially at high surface temperatures have been experimentally observed on the He− Cu(001) system. These deviations have been analyzed with special focus on the role of the interaction potentials and scattering from surface defects. 19 Multiphonon and resonance effects concerning the dependency of the initial particle energy on the DW-factor for a coupled channel approach which cannot be described by the Born approximation are described in a comparative work by Brenig. 20 Multiphonon effects in the evolution of the DW factor have already been observed on standard scattering targets like LiF 21 and still lack a proper treatment in the standard DW model. The surface Debye temperature of Bi(111) has been determined to be Θ D = 71(+7/−5) K using LEED and Θ D = (84 ± 8) K using the specular beam in HAS experiments. 3,8 The surface Debye temperature determined from the attenuation of the first-order diffraction peaks was Θ D ′ = (75 ± 8) K. In the literature, the appearance of two different values for the surface Debye temperature for specular and scattered contributions has been justified because the DW factor relies on scattering processes without momentum transfer parallel to the surface. In particular, the specular HAS value reproduces the theoretical approximation of van Delft 22 that estimates the surface Debye temperature to be lower than the bulk value, which is 120 K, 23,24 by a multiplicative factor of 1/√2. Static Surface Corrugation. For antimony, the elastic coupling parameters for a lattice with the same surface periodicity have been calculated in a previous publication 25 for a corrugation represented as a sum of cosine functions from a Fourier expansion up to the second term, such as ξ ξ π π ξ π with ξ 0 as the corrugation amplitude. The same corrugation function is, with a different lattice constant, assumed for the Bi(111) surface. By assuming a corrugated Morse potential written as the intrablock coupling is given by ν where α = 2κξ 0 . Equation 24 is only exact for a corrugated Morse interaction potential. While modified versions of the interaction potential describe the long-range interaction more accurately, the sole change of V 0 in eqs 8 and 9 to a hybrid potential poses a certain inconsistency. As the inclusion of inelastic effects seems to eliminate the limitations of the Morse potential to model second-order diffraction intensities, 4,25 the usage of the plain Morse potential avoids these inconsistencies while producing excellent results. A more in-depth analysis of a more complicated potential structure calculated numerically from ab initio simulations might further improve the results.

■ RESULTS AND DISCUSSION
Inelastic TOF and Interaction Potential. Previous investigations of the He−Bi(111) interaction potential 27 revealed three well-defined bound state energies when a 9−3 interaction potential model was used. However, preceding close-coupling studies using various potential shapes on Sb(111) 25 suggest that Morse-or Morse-like potential functions are much more suitable for representing the bound state energies of semimetal surfaces.
To propose a more accurate Morse-like interaction potential, features resulting from inelastic resonance effects in time-offlight (TOF) spectra were analyzed to identify an additional bound state level at smaller bound state energies. Figure 2 illustrates one of the spectra with an isolated feature originating from the fourth identified bound state as listed in Table 1.
The last line of Table 1 lists the obtained bound state energy levels of the fitted first Fourier coefficient of the corrugated Morse potential with a potential depth D of (7.898 ± 0.126) meV and a potential stiffness κ of (0.884 ± 0.024) Å −1 . With the highest identified bound state level much closer to the threshold, the attractive part of the fitted potential may be considered to describe the real interaction more accurately. eCC and iCC Analysis of Bi(111). Previous investigations 3 treated the electron density corrugation of the Bi(111) surface from the helium atom scattering (HAS) data using the GR method and the Eikonal approximation, including the Beeby correction. Thermal attenuation effects were included using the DW factor, with two different surface Debye temperatures to account for the two different attenuation features obtained from the measurements. Because the surface Debye temperature is an intrinsic surface property, given by the maximum energy of the phonons in the Debye model, it seems unsatisfactory to include a second temperature in order to account for the different attenuation of the scattering channels. Thus, all of the diffraction intensities issued from the eCC calculations plus the DW factor (eCC+DW) were achieved using only one surface Debye temperature Θ D = 85 K.
As mentioned previously, in order to reduce the number of fitting parameters in the CC calculations, only the surface corrugation amplitude was considered in the eCC calculations, while in the iCC calculations the parameter space was extended to include also the Gaussian cutoff value Q c . All six angular diffraction spectra (three temperatures at the two distinguishable high-symmetry directions ΓM and ΓK) were fitted separately using both methods. In all cases, the overall deviation of the measured diffraction intensities from the calculated intensities  (26) with N being the total number of scattering intensities considered per fit was minimized. The optimization algorithm is terminated after a relative accuracy of 0.1% in all of the considered parameters and conditions. Figure 3 displays the experimental diffraction intensities as well as the best-fitting eCC + DW (red stars) and iCC results (blue downward triangles). To obtain comparable values, the experimental diffraction peak areas were normalized to their respective specular peak area for each spectrum separately. As can be seen in both directions, but especially at the hightemperature measurement in the ΓK direction, the iCC method can almost perfectly account for the emerging asymmetry at varying temperatures, a feature that is vastly impossible for DW-attenuated features. In Table 2, for comparison, the corrugation amplitudes averaged over the three measured surface temperatures for each high symmetry direction given in percentage of the lattice length are listed for both applied calculation methods.
The corrugation amplitude values calculated by both methods are significantly lower than the ones obtained by approximate methods, 3 with the GR method assuming around  The Journal of Physical Chemistry C Article 10% and the Eikonal approximation settling around 11% of the lattice constant a. While both methods when averaged over both directions assume a peak-to-peak corrugation height of around 4.5−5.0%, the variance of the simple elastic model is somewhat higher. In addition, the direction-specific corrugation heights differ significantly from each other in the case of the eCC + DW calculations, while both calculations performed with the iCC method settle around the same corrugation value. The cutoff value implies a width of 3.1 Å, around two-thirds the value of the lattice constant, which justifies the use of the expansion in eq 4 within the single phonon approximation. The inelastic calculations were carried out using 150 scattering channels in three Floquet blocks in a grid going from −7 to +16 Å. The shortest wavelength of the open channels is always described by 100 points. The phonon frequency average is carried out by integration via a simple weighted Legendre quadrature from zero frequency to the Debye frequency with 10 evaluated roots. All parameters were tested for convergence for each of the scattering spectra involved. In particular for the highest temperature involved, the extension to five Floquet blocks was evaluated and found to present no advantage. While the small Debye temperature of Bi(111) influences both the DW factor and the inelastic coupling constant, the iCC approach is the only one that can account for deviation from the strict exponential characteristic of the elastic diffraction peaks with the surface temperature assumed up to now. Figure 4 shows the temperature-dependent attenuation as calculated by the iCC method in comparison with a simple DW factor for both specular and first-order diffraction intensities. The parallel but shifted behavior confirms a DW-like attenuation of the diffraction peaks with the same Debye temperature as for the specular. This also poses as an internal test of our iCC calculation. The left panel of Figure 4 confirms that the attenuation of the specular contribution in both directions follow a DW behavior with a Debye temperature of 85.9 K even though the set Debye temperature of the system introduced in the inelastic coupling constants in eq 16 was 85 K.
The experimental condition of a constant source-detector angle causes the angular spectrum to be recorded while changing the angle of incidence. This so-called "moving threshold" situation causes the beam on one side of the specular contribution to encounter a different scattering scheme Figure 3. Measured and calculated diffraction peak intensities in both distinguishable lattice directions at three surface temperatures and a beam energy of 17 meV. Black dots signify measured peak areas, red stars signify calculated peak intensities using elastic close-coupling with a DW attenuation, and blue downward triangles signify calculated peak intensities using the inelastic close-coupling approach. The "order" of the scattering peak refers to the number of reciprocal lattice vectors needed when fulfilling the Bragg condition (eq 7). Upper panel: Angular scans in ΓM direction at three different surface temperatures. Lower panel: Angular scans in ΓK direction at three different surface temperatures.  The Journal of Physical Chemistry C Article as on the other side, enabling the iCC calculation with its higher adaptability due to the inclusion of inelastic contributions to also model some of the experimentally encountered peak asymmetries as for example in the 400 K surface temperature measurement in ΓK direction depicted in Figure  3. Obviously, there are other experimental sources of asymmetry, as, for example, sample alignment or the changing visible surface area from the detector while rotating the sample. An extension of the coupling calculations, including a complete treatment of the overall force eq 11 and the correct geometry on the scattered helium atom in surface parallel directions, could account for the different polarization directions and improve the quality of the calculated scattering intensities even further. The overall ability of the iCC calculations to model the measured scattering features could be vastly improved if previously determined interaction potential parameters would also enter the fitting procedure directly, instead of being predetermined solely from bound state feature fittings. However, expanding the included parameter space to four dimensions (corrugation height, Gaussian cutoff, potential depth, potential stiffness) becomes prohibitive. Furthermore, including a more realistic interaction potential shape as well as mode-dependent lattice displacements, probably determined by ab initio approaches, would promote the iCC method into a remarkable tool for simulating the effects of inelastic scattering contributions in temperature dependent measurements. A further, essential advancement will be the inclusion of finite phonon momentum in the iCC calculations, extending the Bragg condition to K i − K f = G + Q. A so-enhanced inelastic scattering code could predictably be used to model the experimental time-of-flight spectra and extract information about the mode-specific electron−phonon interaction on conducting surfaces. Using the inelastic closecoupling approach to simulate the scattering from surfaces with finite temperatures clearly renders the inclusion of an additional surface Debye temperature futile. By not being bound to an exponential attenuation characteristic, the method is more adaptive and thus better suited to describe the temperaturedependent scattering behavior.