Article

PDF
Access to the PDF text
Advertising


Free Article !

Comptes Rendus Physique
Volume 17, n° 3-4
pages 302-321 (mars 2016)
Doi : 10.1016/j.crhy.2015.11.006
Friedel oscillations: Decoding the hidden physics
Les oscillations de Friedel : decoder la physique cachée
 

Cristina Bena a, b,
a Laboratoire de physique des solides, UMR 8502, bât. 510, 91405 Orsay cedex, France 
b Institut de physique théorique, CEA/Saclay, L'Orme des Merisiers, 91190 Gif-sur-Yvette cedex, France 

Correspondence to: IPhT CEA Saclay, 91190 Gif sur Yvette cedex, France.IPhT CEA SaclayGif sur Yvette cedex91190France
Abstract

We show that the impurity-induced Friedel oscillations allow one to probe in an unexpected and quite remarkable manner the electronic properties of two-dimensional systems such as graphene or high-temperature superconductors. In particular, we show that by studying these oscillations, one can get access not only to the constant-energy maps, but also to more hidden information such as the chiral properties of Dirac electrons in graphene, which cannot be observed directly by other methods. For graphene, this hidden information is revealed by comparing the theoretical predictions with scanning tunneling microscopy experimental measurements of the local density of states.

The full text of this article is available in PDF format.
Résumé

Nous montrons comment les oscillations de Friedel induites par la présence d'une impureté permettent de sonder de manière inattendue et tout à fait remarquable les propriétés électroniques des systèmes bidimensionnels tels que le graphène ou les supraconducteurs à haute température. En particulier, nous montrons que par l'étude de ces oscillations, on peut obtenir accès pas seulement à la structure de bande, mais aussi à des informations cachées comme les propriétés chirales des électrons de Dirac dans le graphène, qui ne peuvent pas être observés directement par d'autres méthodes. Pour le graphène, cette information est révélée par la comparaison entre les prédictions théoriques et les mesures expérimentales de la densité locale d'états par la microscopie à effet tunnel.

The full text of this article is available in PDF format.

Keywords : Friedel oscillations, Impurity scattering, Graphene

Mots-clés : Oscillations de Friedel, Desordre, Graphene


Introduction

One of the most remarkable feats achieved with an STM, beside the possibility to visualize material surfaces with atomic resolution, was the possibility to image the standing waves associated with the interference of quasi-free electron wavefunctions. The “standing waves” arising in the presence of surface inhomogeneities are also known as Friedel oscillations [[1]], where the term of Friedel oscillations has been first introduced to describe the asymptotic dependence of the perturbed density of states of a free electron gas in the presence of disorder. Such disorder can consist of vacancies (naturally occurring, or induced by irradiation), ad-atoms (magnetic or non-magnetic), substitutional atoms, lattice deformations such as curvatures, local deformations in the lattice structure and the number of bonds, etc. The observation of the Friedel oscillations induced in the vicinity of these impurities allows the illustration of some very important concepts of condensed matter physics. Thus, their analysis provides a direct observation of screening and of electron–electron interaction. Moreover, these oscillations lie at the foundation of the description of the indirect coupling between magnetic moments via the conduction electrons in a metal with the famous Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction potential [[2], [3], [4]], as well as of the long-range adsorbate interaction mediated by a two-dimensional electron gas [[5]].

Following their first description by Friedel, the possibility to use the Friedel oscillations to probe the electronic structure of materials was considered by many others, notably in the case of transitional metals [[6]]. They were mentioned in relation to magnetic impurities in bulk (3D) materials, for which an important damping factor of the amplitude of the oscillations is observed: the oscillations fall off with the distance r as   where α is the dimensionality of the considered electron gas. Furthermore, the observation of the Friedel oscillations was for the first time done indirectly by the observation of the coupling between two magnetic layers with a non-magnetic spacer (host metal) [[7]] for which the theoretical description [[8]] has revealed the Fermi surface of the spacer. The development of the scanning tunneling microscopy (STM) has offered the possibility to study the standing waves in the local density of states (LDOS) which are in fact energy resolved Friedel oscillations for 2D or 1D electron gases which provide longer coherence lengths than those observed in 3D. The local density of surface states has been first obtained by Hasegawa and Avouris [[9]] on confined Shockley states of Au(111) surfaces at room temperature, and by Crommie et al. [[10]] on Cu(111) at 4 K.

Subsequently, the dispersion relation   of surface-state electrons of Ag(111) and Cu(111) permitted the estimation of the surface state inelastic lifetime [[11]]. As large energies are accessible both below and above the Fermi level, it was also possible to study the deviation of the free-electron-like parabolic dispersion, moreover Bürgi et al. [[12]] have defined a method to directly image the potential landscape on Au(111) by STM.

For these measurements, the standing waves have been imaged not only at the Fermi level but also at different energies. This requires the spatial mapping of   at different applied voltages V , which allows one to focus on individual values of the energy specified by V , and avoids an integration over all the wavelengths corresponding to the energies between the considered energy level and the Fermi level. This was clearly demonstrated by Pivetta et al. [[13]] on Ag(111). In this study, standing waves are generated by step edges, and one focuses on a particular reciprocal lattice direction when analyzing the Fermi surface.

In 1997, Sprunger et al. [[14]] demonstrated the possibility to use STM to directly image the Fermi surface contour of a metal surface, by performing the power spectrum (Fourier transform (FT)) of a topographic image of a complex “electron sea” pattern generated by a random distribution of point-like surface defects. These studies were performed around the Fermi level (low bias voltage) and on simple isotropic Fermi surface of Au(111) and Cu(111) surfaces [[15]]. In this seminal paper, Friedel oscillations appear in the Fourier transform as a circular feature centered around the Brillouin zone center with a radius corresponding to  . This corresponds to the simple case of an isotropic Fermi contour centered around the center of the Brillouin zone. As we will see in the next section, this is a direct evidence of the singularity in the Lindhard susceptibility function of the two dimensional electron at scattering momentum vector  .

The dependence of the Lindhard susceptibility on   was subsequently elegantly demonstrated by the observation of standing waves on Be ( ) surfaces [[16]]. In this system the Fermi surface is no longer isotropic and standing waves can not be generated by all step edges (i.e. in all directions). As we will see later, this corresponds to the fact that the wavevectors for the quasiparticles with a given energy can isotropically point in any direction for a circular Fermi contour centered around the Γ, but acquire a dependence on the direction and position in momentum space if the Fermi surface is no longer isotropic, or if it becomes split.

More recently the FT-STM technique (obtaining and analyzing the Fourier transform of STM images) has begun to be applied also to LDOS images (  maps) which combine STM imaging and a spectroscopic measurement of the LDOS as a function of energy and position. This technique has been denoted FT-STS: “FT-Scanning Tunneling Spectroscopy”, though in the literature the terms FT-STM and FT-STS are sometimes both used to describe FT-STS measurements. The FT-STS technique has been applied on high-  superconductors, for a small range of energy around the Fermi level (see, e.g., [[17], [18], [19]]). Another important observation made using the FT-STS was about the spin of the quasiparticles. As two waves with opposite spin directions cannot generate quasiparticle interferences (QPI), Pascual et al. [[20]] demonstrated the ability to probe indirectly the orientation of the spin associated with a Fermi surface of a half magnetic material. The strength of this technique has also been tested on a semi-metal ErSi2 , for which some constant energy contours (CECs) are split into several bands, and where it has been demonstrated as well the possibility to determine the whole 2D band structure in a wide range of energy [[21], [22]].

Here we discuss how this technique has been used to study theoretically and experimentally various systems such as graphene, high-temperature superconductors and interacting one dimensional systems, and how in certain instances this has allowed one to extract important hidden information about the electronic structure of those systems that would not be available experimentally via other types of measurements. We first describe the theoretical technique most widely used to approach this problem, the T-matrix approximation (section 2). We then present applications to various systems such as monolayer and bilayer graphene (section 3), graphene in the quantum Hall regime (section 4), high-temperature superconductors modeled using the d-density wave model (section 5), and one-dimensional Luttinger liquids (section 6). We conclude in section 7.

T-matrix approximation

The T-matrix approximation has been successfully used to calculate the effects of disorder on the spectral properties of an electronic system. The basic theory is described for example in [[23], [24]], and for a few example of how this is applied we mention Refs. [[25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43]]. The basic principle of this technique lies in an infinite perturbative summation of the diagrams resulting from expanding perturbatively the Green function of a system to all orders in the impurity scattering. In order to be able to apply the T-matrix approximation one needs to have the exact form of the tight-binding Hamiltonian of a system. The T-matrix approximation is valid as long as the impurity potential considered is localized, since it is this condition that allows the re-summation of all orders in perturbation theory. For extended impurities, the T-matrix approximation is in general replaced by the Born approximation, for which only the first order term in the impurity potential series is considered; this is equivalent to a perturbative expansion in the impurity potential which is valid only when the impurity scattering is weak.

T-matrix approximation in momentum space

We will briefly review here the principle of the T-matrix approximation [[23], [24]]. The impurity scattering problem can be solved both in the real space and in momentum space. We will first focus on the momentum space calculation. For a given system one can define a finite temperature (imaginary time) generalized Green's function,
(1)G(k1,k2,τ)=−Tre−β(K−Ω)Tτψk1(τ)ψk2+(0) where  ,  , and   is the imaginary time ordering operator. For a translationally-invariant, disorder-free system, the generalized Green's function defined above is non-zero only if   (momentum is conserved), and the generalized Green's function reduces to the standard Green's function which depends only on one momentum. However, the generalized Green's function acquires a non-zero component for   if the system is inhomogeneous, such as in the presence of impurity scattering. This component can be calculated using the T -matrix formulation [[31], [32], [33], [34]], where in what follows we will focus solely on the zero-temperature limit:(2)G(k1,k2,iωn)=G0(k1,iωn)T(k1,k2,iωn)G0(k2,iωn) where(3)G0(k,iωn)=(iωn−Hk)−1 is the unperturbed Green's function of the homogeneous system,   is the Hamiltonian, and(4)T(k1,k2,iωn)=V(k1,k2)+∑k′V(k1,k′)G0(k′,iωn)T(k′,k2,iωn)

As detailed in Fig. 1, this expression stems from a perturbative expansion to all orders in the impurity scattering. If only the first term of the expansion is preserved, one recovers the first order perturbation theory in the impurity scattering potential, also denoted the Born approximation. Oftentimes, one assumes that the impurity scattering potential is very close to a delta function so that V is independent of k and  
(5)Himp=∫dxVc†(x)c(x)δ(x)=∫k→,k→′Vc†(k→)c(k→′) For this case we can solve Eq. ((4)), and obtain(6)T(iωn)=[1−V∫d2kSBZG0(k,iωn)]−1V where   is the area of the first Brillouin zone (BZ) of the system, and the integral over k is performed over the entire Brillouin zone.



Fig. 1


Fig. 1. 

Diagrammatic description of the T-matrix approximation: the Green function for a system in the presence of an impurity can be obtained from adding to the unperturbed Green's function the contributions corresponding to all-order impurity scattering processes. We note that the sum of the contributions appearing to the right of G 0 is identical to T , which allows one to write the self-consistency equation on the third line.

Zoom

In the neighborhood of the impurity, spatial oscillations of the local density of states are induced. The Fourier transform of these oscillations can be related to the generalized Green's function,
(7)ρ(q,ω)=i∑kg(k,q,ω) where  , and   is the generalized retarded Green's function obtained by analytical continuation   from the Fourier transform of the imaginary-time Green's function  .

T-matrix approximation in real space

If one is interested in the real space spectral properties of a system, such as the space dependence of the LDOS, one can use the real space T-matrix formalism. The relations described above can be Fourier transformed to real space, yielding for the retarded space-dependent generalized Green's function of the system
(8)G(R→1,R→2,E)=∫R→G0(R→−R→1,E)T(R→,E)G0(R→2−R→,E) where   is the space-dependent T -matrix. In the same way as for the momentum-dependent formalism, the generalized Green function depends on two (spatial) variables. In the absence of disorder, for a homogeneous system, it is only a function of the difference  , so that only one spatial index is preserved in the notation:  . However, in the presence of disorder, the generalized Green's function will depend on both variables independently.

The spatial dependent Green's function characterizes the propagation between two points in space   and  . If one is not interested in the propagation of a particle between two spatial points, but would rather want to determine the spectral properties, for example the number of allowed states having a given energy at a given position, one needs to focus only on the limit  . Consequently, the LDOS can be obtained from the retarded Green's function by using the conventional relation:
(9)ρ(r→,E)=−ImG(r→,r→,E)

In general one focuses on a delta-function impurity localized at  , which yields
(10)G(r→,r→,E)=G0(−r→,E)T(E)G0(r→,E) where  , and the integral over   is performed on the first BZ, whose area is denoted by  .

Friedel oscillations in graphene: chirality measurements

Graphene has been studied extensively in recent years. Its most fascinating aspect is the existence of linearly-dispersing gapless excitations in the vicinity of the Dirac points. This gives rise to very interesting electronic properties such as Friedel oscillations in the local density of states at low energy which decay as   [[44], [45]], instead the usual   characteristic to two-dimensional systems [[33], [46]]. Comparison with experiments can provide information about graphene's fundamental physics, and about the nature of the impurities.

Here we analyze the Fourier transform of the density of states measurable by Fourier transform scanning tunneling spectroscopy (FTSTS). Such measurements have recently been developed for graphene [[47]], as well as for other two-dimensional materials such as ErSi2 [[21], [22], [48]] and the cuprates [[16], [17], [18], [19], [49]].

Our first observation is that the FTSTS spectra can be used to distinguish between monolayer and bilayer graphene. In particular, for monolayer graphene with a localized (delta-function) impurity potential, the Friedel oscillations generated by intranodal scattering decay as   at low STM bias, consistent with previous analysis [[44], [45]]. In the FTSTS spectra this is manifested by a filled circle of high intensity in the center of the Brillouin zone (BZ), with a radius proportional to the STM bias. On the other hand, the Friedel oscillations generated by the scattering of quasiparticles between different Dirac points decay as  . In the FTSTS spectra, these oscillations are translated into circular contours of high intensity centered around the corners of the BZ and around sites of the reciprocal lattice. Due to the form of the underlying Hamiltonian, the distribution of intensity on some of these circles is not rotationally invariant.

For the bilayer system, at low energy the oscillations due to both intranodal and internodal scattering have a   dependence, corresponding to circular lines of high intensity close to the center and the corners of the BZ, and around the sites of the reciprocal lattice. At higher energies, the splitting of the bands for the bilayer sample is also observable in the FTSTS spectra.

We also note that the FTSTS spectra can allow us to distinguish between different types of impurities. For example, for the case of a screened-charge impurity, the effect of internodal scattering is greatly reduced compared to the effect of intranodal scattering. This gives a clear signature in the FTSTS spectra which can be observed in an experiment.

Our last observation is that the FTSTS spectra, besides providing information about the band structure of graphene, can also give insight into the underlying Hamiltonian. In particular, the shift of the decay of the Friedel oscillations from   to  , and the rotational asymmetry of some of the high-intensity spots, are strongly dependent on the peculiar form of the tight-binding Hamiltonian, and cannot be deduced solely from band-structure arguments.

The tight-binding Hamiltonian for monolayer graphene is:
(11)H=∫d2k→[ak→†bk→f(k→)+h.c.] where the operators  ,   correspond to creating electrons on the sublattice A and B respectively, and  . Here  ,  ,  , t is the nearest-neighbor hopping amplitude, and a is the spacing between two adjacent carbon atoms, which we are setting to 1.

We will use this form of the Hamiltonian when performing our numerical analysis of the FTSTS spectra. However, it is useful to expand the Hamiltonian close to the corners of the BZ, which we also denote as nodes or “Dirac points”, and use the linearized form to solve the problem analytically at low energies. The momenta of the six corners of the Brillouin zone are given by  ,  ,  . Close to each corner, m , of the BZ we can write  , where   denotes the distance from the respective corner. Also  ,  ,   and  ,  .

The corresponding Green function,  , derived from the tight-binding Hamiltonian in Eq. ((11)) can be expanded at low energy around the six nodes (denoted m ), and in the   (A, B) sublattice basis can be written as:
(12)G(k→,ω)≈Gm(k→,ω)=1ω2−|ϕ˜m(k→)|2(ω+iδϕ˜m(k→)ϕ˜m⁎(k→)ω+iδ) where δ is the quasiparticle inverse lifetime. The Fourier transform of the linearized Green function is given by:(13)Gm(r→,ω)∝ω(H0(1)(z)iϕm(r→)H1(1)(z)iϕm⁎(r→)H1(1)(z)H0(1)(z)) where  ,   are Hankel functions,  , and  .

We first focus on a delta-function impurity localized on an atom belonging to sublattice A . In the   basis, the impurity potential matrix V has only one non-zero component  . We start with a T -matrix study of the full Hamiltonian presented in Eq. ((11)), and we analyze our results numerically for various energies. The resulting FTSTS spectra (corresponding to the real part of the Fourier transform of the LDOS) are plotted in Fig. 2. There are several interesting features that should be noted. First, there are regions of high intensity in the FTSTS spectra corresponding to intranodal quasiparticle scattering (central region) and internodal scattering (outer regions). Similar features have also been observed experimentally [[47]]. The high-intensity regions that we find are point-like at zero energy, and acquire distinct features as one increases the energy (STM bias). Thus, at low energy the central high-intensity region is a filled circle, while the outer regions are empty. Also, the rotational symmetry of the high-intensity regions located at the corners of the BZ is broken, while it is preserved for the high-intensity regions centered on sites of the reciprocal-lattice. With increasing the energy even higher, other effects such as the changing of the shape of the equal-energy contours from circular to triangular (trigonal warping) start playing an important role. At very high energy the FTSTS intensity map becomes quite intricate. Note that the negative values in the FTSTS are not surprising: on one hand we plot the FT of the LDOS due to the impurity which can be positive or negative as long as the total LDOS is positive. Moreover, for a system with hexagonal symmetry taking the FT may yield both positive and negative values, as well an imaginary component that we do not depict here. In experiments oftentimes what is considered is the absolute value of the FT of the LDOS to avoid these problems.



Fig. 2


Fig. 2. 

FTSTS spectra for a monolayer graphene sample with a single delta-function impurity. Fig. 2a), Fig. 2c) and Fig. 2d) correspond to energies 0.15 t , 0.6t , and 1.2t respectively, at δ =0.07t . The BZ is indicated by dashed lines. The actual lowest (0) and highest (1) normalized values of the FTSTS intensity correspond to different interval of values for each energy ((−1.3,2.6) for 0.15t , (−0.8,5.9) for 0.6t , and (−6.2,7.2) for 1.2t , in arbitrary units). Fig. 2b) shows a cross section of the FTSTS intensity as a function of k x for k y =0, and for energy 0.15t .

Zoom

We now turn to the analytical study of the dependence of the LDOS on the relative position with respect to the impurity ( ) at low energies. In this range, the physics is dominated by linearly dispersing quasiparticles close to the Dirac points. We find that spatial variations of the LDOS due to the impurity are given by:
(14)ρ(r→,E)∝−Im[G(−r→,E)T(E)G(r→,E)]≈−∑m,nIm[ei(K→m−K→n)⋅r→Gm(−r→,E)T(E)Gn(r→,E)] where   denote the corresponding Dirac points. Here   is the T -matrix,  , where I is the   identity matrix, and the integral over   is performed on the BZ, whose area is  .

Using Eq. ((13)) and expanding the Hankel functions to leading order in  , we find that far from the impurity ( ), the corrections to the local density of states due to scattering between the nodes m and n are given by:
(15)ρmn(r→,ω)∝ωrIm{t(ω)ei(K→m−K→n)⋅r→+2iωr/vi[1−ϕm⁎(r→)ϕn(r→)]} where   is the non-zero element of the T -matrix ( ), and we used the fact that  .

In the case of intranodal scattering ( ) the above expression vanishes and the LDOS is dominated by the next leading correction  . This is different from what usual wisdom would suggest for a two-dimensional system (  decay) [[33], [46]], and has also been described in Refs. [[44], [45]]. We should note that the two-dimensional FT of   is roughly  . This corresponds to a filled circle of high intensity in the FTSTS spectrum, which is consistent with the results of our numerical analysis for the central region of high intensity.

Nevertheless, for the decay of the Friedel oscillations generated by internodal scattering ( ), the leading order behavior is  . The FT of   is  , which translates into empty circles of high intensity in the FTSTS spectra, which is consistent with our numerical analysis. However, since the inverse quasiparticle lifetime δ is finite, there will be some broadening of the resonances and some weight inside the circular contours.

For some of the Friedel oscillations generated by internodal scattering, such as the ones between neighboring node pairs (e.g.  ), the rotational symmetry is broken:  . However, for next-to-nearest-neighbor node pairs (e.g.  ),  , and the oscillations are rotationally invariant. This is consistent with the results of our numerical analysis: the high-intensity regions centered on sites of the reciprocal-lattice are rotationally symmetric, while the high-intensity regions close to the corners of the BZ are not.

We now switch gears and consider the case of a single impurity in bilayer graphene. The bilayer graphene consists of two graphene layers stacked on top of each other such that the atoms in the sublattice A of the first layer occur naturally directly on top of the atoms in the sublattice   of the second layer [[50]], with a tunneling coupling of  . We consider the case of an impurity located on the sublattice A . The case of a single impurity located on a site of a different type, as well as the case of multiple impurities will be presented elsewhere. The resulting FTSTS spectra for the LDOS in the top layer are presented in Fig. 3.



Fig. 3


Fig. 3. 

FTSTS spectra for a bilayer sample. Fig. 3.a) and Fig. 3.c) depict the FTSTS intensity in arbitrary units at energies 0.1t , 0.4t , with t p =0.3t , and δ =0.05t . The actual lowest (0) and highest (1) values of the FTSTS normalized intensity correspond to an actual interval of values of (−6.9,3.6) (a) and (−6.0,7.4) (c). Fig. 3.b) and Fig. 3.d) are cross-sections at k y =0 of Fig. 3.a) and Fig. 3.c) depicting the FTSTS intensity as a function of k x .

Zoom

Note that there are similarities and discrepancies between the monolayer and bilayer cases. Like in the monolayer case, there are areas of high intensity centered on the corners of the BZ, as well as on the sites of the reciprocal lattice. The main difference at low energy is that the central region of high intensity is an empty circle, and not a full circle (as for the monolayer case). At high energy, we also note a doubling of the number of high intensity lines corresponding to the doubling of the number of bands.

An analytical study can be performed at low energies starting from the expansion of the Hamiltonian around the Dirac points m . In the sublattice basis   this yields [[51]]:
(16)Hmbilayer(k→)=(0[ϕ˜m(k→)]2[ϕ˜m⁎(k→)]20) where for simplicity we have set the effective mass of the quadratic spectrum to 1. The corresponding Green function in real space is given by:(17)Gm(r→,ω)∝(H0(1)(z)−[ϕm(r→)]2H2(1)(z)−[ϕm⁎(r→)]2H2(1)(z)H0(1)(z)) where we have denoted  . Starting from Eq. ((14)), we perform a similar analysis to the case of monolayer graphene. Thus we note that at large distances ( ), as opposed to the monolayer case, the leading ( ) contribution for intranodal scattering is non-vanishing:(18)ρm(r→,ω)∝1r|ω|cos⁡(r|ω|/v) This is consistent with the appearance of an empty circular contour at the center of the BZ, as opposed to the filled circle for the monolayer case. The leading contribution to the decay of the oscillations due to internodal scattering is also  . Note that in the monolayer case the amplitude of the Friedel oscillations corresponding to intranodal scattering is independent of energy, and the amplitude of the Friedel oscillations corresponding to internodal scattering increases as  , while for the bilayer case all oscillations decrease with energy as  .

These theoretical predictions have been confirmed experimentally [[47]]; these results are presented in detail in an accompanying article [[52]]. Indeed circular contours of high intensity centered around the corners of the BZ have been observed in monolayer graphene, as well as in bilayer graphene indicating internodal Friedel oscillations decaying as  . Moreover a filled high-intensity circular region at the center of the BZ has been observed in monolayer graphene indicating intranodal Friedel oscillations decaying as  , and an empty circular contour in bilayer graphene indicating intranodal Friedel oscillations decaying as  , as theoretically predicted [[34]]. Such measurements have also been performed in other groups [[53]], in Fig. 4 we present the results obtained by these measurement in bilayer graphene, again this shows circular contours of high intensity centered around the corners and the center of the BZ.



Fig. 4


Fig. 4. 

a) STM topographic image (10×10 nm2, −17 meV, 1 nA) showing an isolated defect. b) The FFT power spectrum of the 2D topographic image in a) as a function of  . The features indicated by the arrows correspond to inter-valley scattering.

Zoom

To conclude, we have computed the effect of single-impurity scattering on the LDOS, in particular on the Fourier transform of the LDOS measurable experimentally by FTSTS. We have found that the FTSTS spectra in the vicinity of an impurity are a very good tool to distinguish between monolayer and bilayer graphene. In particular, for monolayer graphene, the Friedel oscillations due to intranodal scattering decay as   and are rotationally invariant. In the FTSTS spectra they correspond to a filled high-intensity circular region at the center of the BZ. On the other hand, the Friedel oscillations due to internodal scattering decay as  . They lead to different FTSTS features, such as circular contours of high intensity centered around the corners of the BZ and on sites of the reciprocal lattice. Some of these contours display a breaking of rotational invariance. For the bilayer case both the internodal and the intranodal Friedel oscillations decay as  .

We have also showed that the FTSTS spectra can be used to distinguish between different types of impurities, for example between a delta-function impurity and a screened Coulomb scatterer.

Last, but not least, we have noticed that, while the FTSTS spectra in the presence of an impurity can give information on the band structure, they are not fully determined by it, but also contain very important information about the specific form of the Hamiltonian. We believe that this feature is very important, and could be also used in the case of cuprates to understand the physics of high temperature superconductivity. It would be interesting for example to compare ARPES measurements and FTSTS spectra in the presence of an impurity. For graphene, this would clearly distinguish the effects of the Hamiltonian from the effects of the band structure. For other systems such as the cuprates the differences between ARPES and FTSTS will be more complex, and may provide even deeper insights into the underlying physics. Note also that the resonances observed in the FTSTS can be related to the singularities of the Lindhard susceptibility which arise for momenta relating points in the Brillouin zone situated on equal-energy contours (e.g.,   at the Fermi energy), however the exact form of the FTSTS resonances can be obtained solely from the T-matrix calculations.

Friedel oscillations and Landau levels

The LDOS in the presence of impurity scattering in graphene can retrieve important information about the physics of its quasiparticles, specifically about their wavefunction. In what follows we use this to study the wavefunctions of the graphene quasiparticles in the quantum Hall effect (QHE) regime. This issue is of particular importance for understanding why the QHE arises in exfoliated, but not in epitaxial graphene, despite the presence of Landau levels (LLs) in both their spectra.

For a given value of the magnetic field we calculate the FTSTS spectra when the tip bias matches the energy of the Landau levels, focusing in particular on the zeroth and the second LL. For energies of the tip situated between two Landau levels, within the approximations we use (Landau levels constant throughout the sample, single-impurity scattering, energy conservation), the intensity of the FTSTS spectra is negligible. We observe that for the zeroth LL, only high-intensity regions corresponding to intra-nodal scattering and to scattering between equivalent nodes are present (at the center of the BZ and reciprocal lattice points respectively), while no features can be identified at the corners of the BZ corresponding to scattering between nonequivalent nodes. The scattering features are rotationally symmetric and decay in a Gaussian manner, consistent with the wavefunction of the quasiparticles in the zeroth LL of graphene.

For higher Landau levels, our calculations reveal both intra-nodal and inter-nodal scattering features. The intra-nodal ones, as well as the ones corresponding to scattering between equivalent nodes are rotationally symmetric. However, the patterns corresponding to scattering between nonequivalent nodes are asymmetric. Their asymmetry is a consequence of the chirality of the graphene quasiparticles. The radial dependence of these features stems from the electronic wavefunction of the Landau levels,1 and shows intensity minima and maxima. The scale of these fluctuations, as well as the scale associated with the Gaussian decay, are proportional to the inverse magnetic length.

In the absence of magnetic field, the eigenfunctions of the above Hamiltonian have been extensively studied (see e.g. [[54]] and references therein). In the presence of a large magnetic field (QHE regime), the eigenfunctions have also been determined in Refs. [[55], [56]] by noting that the Hamiltonian reduces to the Hamiltonian of the harmonic oscillator. The diagonalization of the Hamiltonian for  ,  , and   (around the point  ) can be done by building the eigenfunction [[55], [56]]:
(19)Ψ(r→)=∑keikxL(0ϕ0(y−klB2))ck,−1+∑k,n,αeikx2L(ϕn(y−klB2)αϕn+1(y−klB2))ck,n,α where   is the magnetic length, and   are the eigenfunctions of the one-dimensional harmonic oscillator (  are the usual Hermite polynomials). Also,  , the  's are the annihilation operators for quasiparticles in the  'st LL, with wavenumber k along the x direction and band α , and   is the annihilation operator for a quasiparticle in the zeroth LL. In the new “c ”-operator basis the Hamiltonian is diagonal, and the Green functions are:(20)Gn,k,α(ω)=〈cn,k,α†(ω)cn,k,α(ω)〉=1ω+iδ−En,k,α

We generalize this eigenfunction to take into account all the Dirac points (for the first BZ this reduces to the wavefunctions described in [[56]]) and we obtain:
(21)Ψ(r→)=∑ξ=±1,μ,νΨμνξ(r→)eiK→μνξ⋅r→ with(22)Ψμνξ(r→)=∑keikx2L((1−ξ)ϕ0(y−klB2)(1+ξ)ϕ0(y−klB2))ck,−1,ξμν+∑n,k,αeikx22L((1+ξ)ϕn(y−klB2)−α(1−ξ)ϕn+1(y−klB2)α(1+ξ)ϕn+1(y−klB2)+(1−ξ)ϕn(y−klB2))ck,n,α,ξμν The  -operators are annihilation operators that beside the wavenumber k , band α and LL index n have also the valley indices ξ , μ and ν . Their unperturbed correlation functions do not depend on the valley indices:(23)Gn,k,α,ξμν(ω)=〈cn,k,α,ξ†μν(ω)cn,k,α,ξμν(ω)〉=1ω+iδ−αEn and, in the absence of disorder, correlators of operators connecting two different valleys are zero.

We introduce a delta-function impurity localized on an atom belonging for example to the A sublattice with an impurity potential
(24)V=uΨA†(r→=0)ΨA(r→=0) where   is given by Eq. ((21)). Using Eqs. ((21)), ((23)) and the Born approximation, we find that in the presence of the impurity potential V ,   due to impurity scattering is given by:(25)δρ0(q→)∝∑ξ,ξ′,μ,ν,μ′,ν′∫−∞∞dk∫−∞∞dye−iqyyei(Kμ′ν′yξ′−Kμνyξ)y(1−ξ)2(1−ξ′)2×1(ω+iδ)2ϕ0(y−klB2)ϕ0(y−k′lB2)ϕ0(−klB2)ϕ0(−k′lB2)|k′=k−qx+Kμ′ν′xξ′−Kμνxξ for the zeroth LL ( ).

In this calculation, we have neglected the contributions coming from a quasiparticle being scattered between different LLs, as we work under the simplifying assumptions that the dominant scattering mechanism is elastic, and that the energy of each LL is constant throughout the sample. We have also focused on energies matching the energies of the Landau levels.

We can compute   analytically:
(26)δρ0(q→)=∑μ,ν,μ′,ν′δρ˜0(q→+K→μν−1−K→μ′ν′−1) where(27)δρ˜0(q→)∝∫−∞∞dye−iqyye−(y/lB)2/2e−(y/lB+qxlb)2/2∫−∞∞dkeiqyklB2e−(klB)2/2e−(klB+qxlb)2/2∝e−lB2(qx2+qy2)/2

The corresponding spectrum is plotted in Fig. 5. We note that there are regions of high intensity corresponding to quasiparticle scattering between two equivalent nodes (at the center of the first BZ and all equivalent points related via translation by a reciprocal lattice vector), but no high-intensity regions at the corners of the BZ corresponding to scattering between nonequivalent nodes. This can be seen directly from Eq. ((25)) as only the term proportional to   appears in   in the zeroth LL, and this term is nonzero only if  . This is related to the fact that in the zeroth LL the electronic wavefunctions have only one non-zero component (A or B), depending on the type of node ( ) on which the electron sits. In order for the electron to scatter between two nodes, it needs to be able to change the sublattice index during the scattering process. However, for the type of impurities we consider (potential disorder), this is not possible. The shape of the observed regions is rotationally symmetric, and the intensity decays with the distance from the center in a manner characteristic to the decay of the ground state of the harmonic oscillator ( ).



Fig. 5


Fig. 5. 

FTSTS spectrum as a function of   for a monolayer graphene sample with a single delta-function impurity, for an energy E =0 inside the zeroth LL, and l B /a =7. The BZ is indicated by dashed lines.

Zoom

For an energy corresponding to the  'st LL, the FT of the LDOS is given by:
(28)δρ˜αnμνμ′ν′ξξ′(q→)∝{(1+ξ)2(1+ξ′)2|Inn(q→)|2+(1−ξ)2(1−ξ′)2|In+1,n+1(q→)|2+(1+ξ)2(1−ξ′)2|In,n+1(q→)|2+(1−ξ)2(1+ξ′)2|In+1,n(q→)|2+e−2πi(μ+ν−μ′−ν′)/3[(1−ξ)2(1−ξ′)2Inn(q→)In+1,n+1⁎(q→)−(1−ξ)2(1+ξ′)2In,n+1(q→)In+1,n⁎(q→)−(1+ξ)2(1−ξ′)2In+1,n(q→)In,n+1⁎(q→)+(1+ξ)2(1+ξ′)2In+1,n+1(q→)In,n⁎(q→)]} where we can compute the   integrals analytically:(29)Imn(q→)=π2NM![(σnmqxlB−iqylB)/2]N−MLMN−M(q2lB2/2)e−(qx−iqy)2lB2/4 where   is a Laguerre polynomial,  ,   is the length of the   vector, and M and N denote the smaller and respectively the larger of m and n .

At this point, there are a few observations we can make. The first is that   contains terms proportional to Laguerre polynomials (which in general characterize the overlap between two LL eigenfunctions [[57]]). A Laguerre polynomial shows a number of zeroes (or “nodes”) given by the order of the polynomial (e.g.,   has 1 node,   has 2 nodes, and so on). As these Laguerre polynomials appear in the dependence of the LDOS   on q , we expect   to also have a “node” structure. Indeed, the scattering between nonequivalent Dirac points ( ) gives rise to terms that contain only combinations of   terms (all proportional to a Laguerre polynomial of order n ); these combinations will have n nodes. However, scattering between equivalent Dirac points   gives rise to terms that contain both   and  , and hence contain combinations of Laguerre polynomials of different orders. In general one cannot predict the exact number of nodes of such combinations.

The second observation is that the high-intensity features show, besides the Laguerre polynomial dependence, a Gaussian decay on a scale of a few times the inverse magnetic length  .

The third observation is that the expression in Eq. ((29)) contains rotationally asymmetric phase factors. The terms of the form   cancel in the final expression of   where only products of an integral of the type   and of a complex conjugate of such an integral appear. The other complex phase factors   only appear if  , and hence, as it can be seen from Eq. ((28)), only for processes involving scattering between nonequivalent nodes ( ). This type of scattering corresponds to coupling between the LL wavefunctions   and  . For the non-relativistic QHE this coupling can only arise for the transition of a quasiparticle between two distinct (energy separated) Landau levels. However, we see that for graphene, due to the spinorial structure of the wavefunction, such ( ) coupling occurs naturally inside the same LL, and appears to be a manifestation of the chirality of the graphene quasiparticles.

We have already discussed the particular case of the zeroth Landau level  . To illustrate the observations above we focus on the FTSTS spectra in the second Landau level ( ). This allows us to observe the node structure of the results in more detail than for   (first LL). The   integrals are given by:
(30)I11(q→)=π(2−q2lB2)e−(qx−iqy)2lB2/4I12(q→)=π(qx−iqy)(4−q2lB2)e−(qx−iqy)2lB2/4I21(q→)=−π(qx+iqy)(4−q2lB2)e−(qx−iqy)2lB2/4I22(q→)=π[8+q2lB2(q2lB2−8)]e−(qx−iqy)2lB2/4

In Fig. 5, Fig. 6, we plot the real part of the FT of the LDOS for a large magnetic field  ; these features should be observable for energies smaller or equal to 400 meV, for which the linear approximation is still reasonable [[58]].



Fig. 6


Fig. 6. 

The real part of the FTSTS spectra as a function of   for a monolayer graphene sample with a single delta-function impurity, for an energy equal to the energy of the second LL i.e. ω 400 meV, and l B /a =26.5 (B =50 T). The intensity of the features at the corners of the BZ was multiplied by a factor of two.

Zoom

For both values of the magnetic field there exist regions of high intensity corresponding to both intra-nodal and inter-nodal scattering which exhibit maxima and minima superposed over a Gaussian decay; the distance between the maxima and minima is proportional to the inverse magnetic length. The high-intensity regions corresponding to scattering between nonequivalent Dirac points ( ) have only one zero-intensity node ( ) in their radial dependence, as expected (see Fig. 6). The intensity also goes to zero at the center of these regions ( ) because of the   factors in Eq. ((30)). The high-intensity region corresponding to scattering inside the same Dirac point (located at the center of the BZ) has   nodes, while the high-intensity features corresponding to scattering between equivalent Dirac points ( ) show a structure of minima and maxima similar to that of the central feature, but no zero-intensity node.

Also, as noted above, the features corresponding to scattering between equivalent Dirac points are rotationally symmetric (depend only on the magnitude q ). However, for scattering between nonequivalent Dirac points,   is proportional to  , which breaks rotational symmetry.

To conclude, we have computed the effect of single-impurity scattering on the Fourier transform of the LDOS in the presence of a strong magnetic field. We have found that the FTSTS spectra contain high-intensity regions corresponding both to scattering processes in which quasiparticles remain at the same Dirac point (intra-nodal), and to scattering processes in which quasiparticles hop between different Dirac points (inter-nodal). Both types of processes give rise to features that contain information about the wavefunction of the electrons in the quantum Hall state (decay length, maxima, zeroes, minima). While scattering between equivalent nodes gives rise to rotationally symmetric features, the features coming from scattering between nonequivalent nodes break this symmetry, manifesting the chirality of the quasiparticles. A special situation arises for energies inside the zeroth LL when we observe no features corresponding to scattering between nonequivalent nodes. This should gives rise to a smoother spatial dependence of the LDOS than for the energies corresponding to higher Landau levels.

We focus on a localized impurity, but we expect our results to be quite similar for an extended Coulomb impurity. The most significant difference will be a reduction of the ratio between the intensity of the features away from the center and at the center. It would be interesting to see what happens if other types of disorder are considered which affect not only the electronic density, but also the hopping parameters in the neighborhood of the impurity. Unlike the potential disorder, the hopping disorder couples the nonequivalent valleys even in the zeroth LL, and one expects to observe the corresponding inter-nodal scattering features in the FTSTS spectra for all LLs. Therefore we propose to use the ratio between the intra-nodal and inter-nodal scattering features not only as a good indicator of the extension of the impurity potential, but also in the zeroth LL as an indicator of the form of the impurity potential (potential disorder versus hopping disorder).

The relevant physical regime for the magnetic field is between 20–50 T, and the dominant features we describe (for the second LL) are expected to arise for energies of the tip matching the energies of the LLs, i.e. 200–400 meV. It would be interesting to study what happens for energies of the tip situated between two Landau levels in the presence of multiple impurities. Because of disorder, the LLs are expected to broaden, which translates spatially into a LL energy that fluctuates spatially. Thus, for some energies, one may have overlapping contributions from various LL's, as well as more complicated (and not so neat) features arising in the FTSTS spectra.

We have found that the FTSTS spectra can give information about the electronic wavefunction in graphene at energies matching the LL energy (decay length, nodes, etc.), as well as about the nature of the disorder. We hope that a comparison between this theoretical study and experiments will shed light on the nature of the electronic states in epitaxial graphene under high magnetic field, where quantum-Hall features such as the LL's in the DOS are present, but no quantum Hall effect is observed. In particular it would be interesting to test whether the quasiparticles in epitaxial graphene are indeed described by typical quantum-Hall wavefunctions, and if so, whether it is possible to establish the nature of disorder in both epitaxial and in exfoliated graphene. Identifying the form of the wavefunction, as well as the nature of disorder will help understand why the quantum Hall effect is not observed in epitaxial graphene.

Friedel oscillations in high-temperature superconductors: testing possible models

One of the intriguing features of the cuprate high temperature superconductors is the existence of a pseudogap phase in the normal state [[59]]. Elucidating the nature of the pseudogap would be an important step in understanding the physics of high-temperature superconductivity. Among many scenarios, one concrete proposal is that the pseudogap is due to a hidden broken symmetry [[60]] of   type in the particle–hole channel [[61]]. Here we examine the consequences of the presence of such d-density wave (DDW) ordering on the Friedel oscillations in the pseudogap phase, as measurable in STM experiments of the local density of states and its Fourier transform (measured via Fourier transform scanning tunneling spectroscopy – FT-STS). Such STM measurements [[16], [17], [18], [19], [49]] on Bi-2122 have been performed in the superconducting state at low temperatures, and the results have been interpreted as interference patterns due to elastic scattering of the quasiparticles from impurities. [[17], [25], [31], [62], [63], [64]]. The underlying explanation for these observations is the fact that the scattering of quasiparticles between regions of the Brillouin zone with high densities of states yields peaks in FT-STS. In a d-wave superconductor (DSC), the regions of high local density of states (LDOS) are situated at the tips of the banana-shaped contours of the quasiparticle excitation spectrum [[17]]. The wavevectors of the observed peaks are consistent with the wavevectors that connect the tips of these bananas. Thus, the experimental results are consistent with the known pictures of the Fermi surface and the LDOS in the pure DSC state.

Moreover, it has even been possible to perform such experiments in the pseudogap regime above   [[19]] and also to study their doping dependence in the superconducting state [[18]]. These experiments give us a window on the gap (or depletion of low-energy states) which is present in these materials and on its variation from place to place in a given sample. One striking aspect of all of these measurements is the presence of peaks in the Fourier transform of the local electronic density of states (LDOS) at wavevectors  , and  , where the wavelength λ is between 4 and 7 lattice spacings [[18]]. This problem has been analyzed from a variety of theoretical perspectives [[29], [31], [32], [46], [64], [65], [66], [67], [68], [69], [70]].

Here, we investigate the possible form of STM spectra in both the superconducting and normal states from the perspective of DDW. The T-matrix formalism combined with a numerical analysis, are used to obtain the quasiparticle interference patterns of the various ordered states. In the coexisting state, at energies below the DSC gap, the patterns are almost identical to those in the pure DSC state with the same DSC gap. However, they are significantly different for energies greater than or equal to the DSC gap. This transition at an energy around the DSC gap can be used to test the nature of the superconducting state of the underdoped cuprates by scanning tunneling microscopy.

In the presence of DDW and DSC order, the general mean field Hamiltonian for a high   superconductor is given by
(31)H=∑kσ[ϵ(k)−μ]ckσ+ckσ+∑k[∑σiW(k)ckσ+ck+Q,σ+Δ(k)ck↑+c−k↓+]+h.c. where   and   are the DDW order (see also Refs. [[61]]) and DSC order parameters, respectively, and  . The sums over k include all the wave vectors in the first Brillouin zone (BZ),  ,  . We take the lattice constant  .(32)ϵk=t0+t1(cos⁡kx+cos⁡ky)/2+t2cos⁡kxcos⁡ky+t3(cos⁡2kx+cos⁡2ky)/2+t5cos⁡2kxcos⁡2ky+t4(cos⁡2kxcos⁡ky+cos⁡2kycos⁡kx)/2 and  ,  . We also take the chemical potential shift (from the above dispersion) to be  . The chemical potential is chosen such that – and this is consistent with the ARPES measurements [[71]] – no electron pockets open in the band structure. The equal energy contours for this ARPES band structure as well as for the simpler band structure used in Ref. [[32]] in which   vanish are given for comparison in Fig. 7. We focus on the case of non-magnetic impurity scattering.



Fig. 7


Fig. 7. 

Equal energy contours in the momentum space (k x ,k y ), for the pure DDW state with a) t t band structure b) ARPES consistent band structure. The magnitude of the DDW gap is taken to be 40 meV. The innermost contour corresponds to +40 meV energy, while the outermost contour corresponds to a −40 meV energy. Scattering between regions of high density of states is indicated by the arrows.

Zoom

To simplify the calculation, we introduce a four-component (spinor) Nambu field operator,  . In this new basis, the above Hamiltonian can be written as
(33)H=∑kψk+A(k)ψk where k is summed over half of the original Brillouin zone (reduced Brillouin zone – RBZ), namely,  .   is a four-by-four matrix given by(34)Ak=(ϵ1(k)+ϵ2(k)−μiW(k)Δ(k)0−iW(k)−ϵ1(k)+ϵ2(k)−μ0−Δ(k)Δ⁎(k)0−ϵ1(k)−ϵ2(k)+μiW(k)0−Δ(k)⁎−iW(k)ϵ1(k)−ϵ2(k)+μ) The eigenvalues of   are   and  , where(35)E1,2(k)={(ϵ1(k)2+W(k)2±[ϵ2(k)−μ])2+Δ(k)2}12 We consider impurity scattering of the form:(36)Himp=∑k,k′∑α,β=↑,↓Vkk′αβckα+ck′β Up to a constant, we can write it as(37)Himp=∑k,k′∈RBZψk+V(k,k′)ψk′ where   is a four by four matrix.

The local density of states  , is given by
(38)ρ(q,ω)∼i2π∑k∈RBZg(k,q,ω). where   is defined as follows. Let  . If   is in the RBZ,(39)g(k,q,ω)=∑i=14[Gii(k,k′,siω)−Gii⁎(k′,k,siω)] where   for   and   for  . If   is not in the RBZ, let  . For this case,(40)g(k,q,ω)=∑i=1,3[Gi,i+1(k,k″,siω)−Gi,i+1⁎(k″,k,siω)+Gi+1,i(k,k″,siω)−Gi+1,i⁎(k″,k,siω)] Here   is obtained by analytical continuation   of   from imaginary frequencies to real frequencies. The complexity of the above formula stems from translational symmetry breaking in the presence of DDW order.

We compute the local density of states   using the   impurity scattering matrices. For potential scattering given by a δ -function, they are
(41)V(k,k′)=VN(1100110000−1−100−1−1). for a non-magnetic impurity, and(42)V(k,k′)=VM(1100110000110011) for a magnetic impurity.

Our results are plotted in Fig. 8 for energies between −42 meV and 45 meV for a pure DDW state with a gap of  . The results are reported for the representative values:  . The imaginary part of the energy   is used for the entire numerical calculation. We have checked that the results are unchanged for smaller values of δ . Following Ref. [[31]], a   lattice is used in our analysis, and the results are displayed in the   interval on a   grid for any given frequency. The choice of these parameters is representative. We have repeated our calculations for a number of different set of parameters; the conclusions remain unchanged.



Fig. 8


Fig. 8. 

Quasiparticle interference spectra for a DDW state with DDW gap of 40 meV for non-magnetic impurity scattering V =0.1 eV. The results are displayed for energies ranging from −42 meV to 45 meV on a linear gray scale. Each plot represents the spectral intensity as a function of momentum in the 1st BZ, for |q x |,|q y |<π . The positions of the peaks are indicated by circles.

Zoom

We note the appearance of peaks corresponding to scattering between the tips of the ellipses, as indicated by arrows in Fig. 7. The position of the peaks is marked by circles in Fig. 8.

The dispersion of the wavevectors with energy is plotted in Fig. 9. As indicated, for positive energies, the dispersion is quite small, and the magnitude of the wavevectors ranges from   at 3 meV to   at 39 meV. The dispersion is larger for negative energies, and the peaks are not as well defined. We also note that if one does not shift the chemical potential, while electron pockets open, the dispersion of the peaks position with energy is much smaller and thus more in agreement with the experimental results. This indicates that both the position of the peaks and their energy dispersion are very sensitive to band structure parameters.



Fig. 9


Fig. 9. 

The energy dispersion of the wavevectors of the observed peaks. The magnitude of the wavevectors is plotted in units of 2π , and the energy is measured in meV. The magnitude of the DDW gap is 40 meV.

Zoom

Thus, the position of the peaks resulting from quasiparticle interference in a pure DDW with an ARPES-consistent band structure is similar to experimental observations, while the dispersion with energy is larger. However, this is a crude single impurity scattering T-matrix approximation calculation and other factors may need to be taken into account. Electron–electron interactions have been neglected, but are likely to be important for understanding ARPES experiments [[72]] in which they smear out antinodal quasiparticles. It is hard to imagine that this would not have an impact on LDOS measurements. Furthermore, the renormalization of the order parameter due to disorder in the presence of many impurities may need to be taken into account through a self-consistent Bogoliubov–de Gennes calculation [[73]].

In conclusion, we obtained the quasiparticle interference spectra in a DDW state, which would correspond to the pseudogap phase of the cuprates. We used an ARPES consistent band structure. We observed the emergence of peaks with wavevectors of magnitude   (at positive energies) along the ( ) direction. The magnitude of the wavevectors is similar to those seen in recent STM experiments in the pseudogap phase [[17], [18], [19]]. The energy dispersion of the peaks is larger than what was observed experimentally, though other theoretical and experimental factors may need to be taken into account when trying to connect our observations to the experimental data. However, the peaks seen in STM experiments will always disperse with energy – even in the case of charge order disrupted by impurities [[46], [65]]. The issue at hand is the quantitative one of how much.

Friedel oscillations in Luttinger liquids: effects of interactions

The interplay between interactions and disorder is a long-standing problem in condensed matter physics. In one-dimensional systems, in which the interactions can be treated exactly using the Luttinger liquid theory [[74], [75], [76]] and bosonization [[77], [78]], a lot of progress to understand the effects of impurities has been made over the last 20 years. It was shown that repulsive interactions such as the Coulomb interactions renormalize the impurity strength, such that at low energy even a weak impurity has a very strong effect and can cut the wire into two pieces [[79], [80], [81], [82], [83]]. This translates into a reduction of the local density of states (LDOS) at low energies, and the LDOS decays to zero as a power law [[84], [85], [86], [87], [88]]. At high energies, the effect of the impurity consists in a small power-law correction of the unperturbed LDOS. The two power laws are characterized by two different exponents which depend on the interaction strength.

In this work, we develop an approach that allows one to study the interplay between impurities and interactions and the Freidel oscillations in the vicinity of an impurity for arbitrary size impurities. Our approach is based on writing and solving the Dyson equations – a technique similar to the T-matrix approximation. We consider an infinite homogeneous interacting QW and the corresponding form for the fermionic Green's functions.

We calculate the form of the Friedel oscillations as well as the dependence of the LDOS with energy. For weak impurities we retrieve the expected Luttinger liquid power-law dependence of the impurity with energy at both low and high energies. For strong impurities at either high energies or large distances to impurity, we recover a power-law dependence with the same exponent   as for the weak-impurity regime, consistent with the Luttinger liquid predictions. However, at low energy and small distance, our approach fails to recover the transition to a different power-law exponent characteristic to breaking the wire into two independence pieces, i.e., with exponent ( ). This comes from a drawback in the approximation used in our approach, which consists to neglect terms mixing impurity potentials and Coulomb interactions in the Dyson equation. This nevertheless does not affect the behavior at large distances/energies, and the validity of the main results of this paper, i.e., the dependence of the LDOS in a wire with two impurities and the transition from the weak-impurity regime to the strong-impurity regime.

We consider a one-channel interacting QW with an impurities at positions   The impurity is described by a backward scattering potential   and a forward scattering potential  . The Hamiltonian can be written as  , where   describes the non-interacting QW without impurities:
(43)H0=−iħvF∑r=±r∫−∞∞ψr†(x)∂xψr(x)dx with   is the Fermi velocity,   and   are the creation and annihilation fermionic operators associated with the right movers ( ) and left movers ( ). The Hamiltonian   describes the Coulomb interaction in the wire:(44)Hint=12∫−∞∞∫−∞∞ρˆ(x)V(x,x′)ρˆ(x′)dxdx′ where   is the density operator, and V is the Coulomb potential, which is assumed to be short range due to screening effects by metallic gates or inter-wire coupling. The impurity Hamiltonian contains two types of contribution  , backward-scattering term(45)HB=∑r=±∫−∞∞λB(x)ψr†(x)ψ−r(x)dx and forward-scattering term(46)HF=∑r=±∫−∞∞λF(x)ψr†(x)ψr(x)dx In the following, we assume that the impurity is localized, i.e.  . Notice that in some works dealing with impurities in a Luttinger liquid, the forward-scattering terms were not included with the justification that such terms could be incorporated in the kinetic part [[89], [90], [91], [92]]. This is correct in the weak-impurity limit, but does not hold in the strong-impurity one. Indeed, the density of states is strongly affected by the forward-scattering terms at strong  ; in particular, these terms need to be taken into account explicitly in order to recover the Coulomb-blockade regime.

In order to calculate the form of the Green functions in the presence of impurities, we establish the Dyson equation associated with the Hamiltonian H . Assuming that Coulomb interactions are strongly attenuated with distance and weak in comparison to the energy and neglecting the contributions mixing the impurity potentials and the Coulomb potential, we find
(47)Gr,r′R(x,x′;ω)=grR(x,x′;ω)δr,r′+grR(x,x1;ω)[Γ1BG−r,r′R(x1,x′;ω)+Γ1FGr,r′R(x1,x′;ω)] Consequently, for the Green functions, we derive(48)Gr,rR(x1,x′;ω)=grR(x1,x′;ω)[1−Γ1Fg−rR(x1,x1;ω)][1−Γ1F[grR(x1,x1;ω)+g−rR(x1,x1;ω)]+grR(x1,x1;ω)[(Γ1F)2−(Γ1B)2]g−rR(x1,x1;ω)]−1 and(49)G−r,rR(x1,x′;ω)=g−rR(x1,x1;ω)Γ1BGr,rR(x1,x′;ω)1−Γ1Fg−rR(x1,x1;ω) where   are the Green functions of a clean interacting homogeneous wire, associated with  . They can be obtained in the framework of the Tomonaga–Luttinger theory [[74], [75]]. For an infinite QW with uniform interactions, their form has been derived explicitly in Ref. [[93]]. These Green functions depend on a single chiral index r since the chiral states are eigenstates of the interacting QW in the absence of impurities:(50)grR(x,x′;ω)=−eirkF(x−x′)K2ω+2ħvFωcπΓ(1+γ)(2i|x−x′|ωcaKω+)12−γ[Kγ−12(K|x−x′|ω+iaωc)−sgn(r(x−x′))Kγ+12(K|x−x′|ω+iaωc)] where  ,   is the Fermi momentum,  , a is the small-distance cutoff of the Tomonaga–Luttinger liquid theory [[74], [75]], Γ and K are, respectively, the gamma and modified gamma functions, and  . Here, K is the interaction parameter which is related to the interaction potential by the relation  . In the non-interacting limit,  , Eq. ((50)) recovers the chiral Green's functions of a clean non-interacting QW:(51)grR(x,x′;ω)=−ieirkF(x−x′)ħvF×eiωr(x−x′)/vFΘ(r(x−x′))

The above formulas contain all the information necessary to calculate the LDOS of any chiral wire in the presence of one or two impurities as a function of energy and position. It is in fact a generalization of the solutions obtained in Ref. [[94]] for non-interacting systems to take into account the effects of interactions. An important improvement with respect to Ref. [[94]] is the fact that our equations for the Green functions are chirality resolved in order to include appropriately Coulomb interactions.

In Fig. 10, Fig. 11, we show the profiles of the LDOS for increasing impurity potential in the presence and in the absence of Coulomb interactions. While the LDOS is asymmetric in energy at weak-impurity potential (Fig. 10), it becomes symmetrical for the strong-impurity potential (Fig. 11). The effect of the impurity is to introduce spatial oscillations whose amplitude increases with the impurity potential. In the presence of interactions, the period of these Friedel oscillations is modified, and the value of the LDOS is reduced. As shown in Fig. 12, the amplitude of oscillations and the density of states at the impurity position are both reduced for   (left graph) in comparison to   (right graph).



Fig. 10


Fig. 10. 

LDOS (left graph) in the presence of Coulomb interactions (K =0.7), and (right graph) for a non-interacting wire (K =1) as a function of position (horizontal axis) and energy (vertical axis) for a weak impurity  . We take k F =0.

Zoom



Fig. 11


Fig. 11. 

The same as Fig. 10 for a strong impurity  .

Zoom



Fig. 12


Fig. 12. 

LDOS in the presence of Coulomb interactions, K =0.7 (left graph) and for a non-interacting wire, K =1 (right graph) at ω =0.01ω c , and for   (red dashed lines),   (blue dashed–dotted lines), and   (black solid lines). We take k F =0. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Zoom

The reduction of the LDOS at the impurity position (here  ) when increasing the strength of the interactions is observed for all values of the impurity potential (see the left graph in Fig. 13). However, the LDOS is drastically reduced for the largest impurity potentials for all values of K (black full line); this is because a large impurity effectively cuts the wire into two disconnected pieces.



Fig. 13


Fig. 13. 

LDOS (left graph) at the impurity position (x =0) as a function of K , at ω =0.01ω c , and (right graph) close to the impurity position (x /a =1) as a function of energy at K =0.7. The right graph is plotted in a logarithm scale. On both graphs, we have   (red dashed lines),   (blue dashed–dotted lines), and   (black full lines). We take k F =0. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Zoom

In the right graph of Fig. 13 is shown the LDOS as a function of energy (on a logarithmic scale) for a position close to the impurity. For a weak impurity (red dashed line), the LDOS exhibits a power-law dependence with energy:
(52)ρ0(ω)=|ω|(K+K−1−2)/2πΓ(K+K−12) which is just the density of states of a clean interacting wire. Here, Γ is the gamma function. When the impurity potential increases, the LDOS deviates from this power-law at small energy but converges and oscillates around this power law behavior when the energy increases, as expected (see the dashed–dotted blue line and the black line). This behavior is in full agreement with the results obtained in Refs. [[92]] and [[95]]. However, within our approach (and its limitations), we are not able to recover the expected behavior of the LDOS at low energy/distance and strong impurity potential, i.e., a power law behavior of   characteristic to injecting an electron into the end of a semi-infinite wire. This is due to the fact that we neglect the terms mixing the Coulomb interactions and the impurity potential in the Dyson equation. This approximation is justified when   and   are both weak. For strong  , we need an addition assumption which is  . This is the reason why our approach fails at low energy and strong impurity potentials. At high energy/large distances, our results are valid.

We have developed an approach based on the Dyson equations which has allowed us to study the LDOS of an infinite interacting QW with two impurities of arbitrary strength. For an infinite homogeneous interacting wire with a single impurity, we have calculated the form of the Friedel oscillations as well as the dependence of the LDOS with energy. We have found that for weak impurities, as well as for strong impurities at high energies/large distances, our approach recovers the expected Luttinger liquid power-law dependence; however, it breaks down for strong impurities at low energy/small distance.

We have applied this approach to study the transition from the weak-impurity regime to the strong-impurity regime in a wire with two impurities, focusing in particular on the regime of large distances and energies. We have found that the main effect of interactions is to reduce the amplitude of the Fabry–Perot oscillations in the weak-impurity limit, as well as of the Coulomb-blockade peaks in the strong-impurity limit. In addition, the interactions affect the periodicity of the oscillations and the distance between the peaks. Moreover, we see that strong interactions also reduce to zero the LDOS on the impurity sites. Also, at non-zero   and for strong impurities, our results are consistent with those obtained in Ref. [[96]] for the LDOS of a Luttinger liquid in a box, in particular, we recover an extra-modulation of the LDOS oscillations in the presence of interactions, which is absent in the non-interacting system. This gives an extra confirmation that our approach is valid in that regime.

Conclusions

We have described how the impurity-induced Friedel oscillations allow one to probe specific electronic properties of two-dimensional systems such as graphene or high-temperature superconductors, as well as one-dimensional strongly interacting systems. We have shown that the study of the Friedel oscillations allows one to get information not only about the band structure of a system, such also possible using ARPES, but also more hidden information such as the chirality of the electrons in graphene, the wave-function structure in quantum Hall systems, or information about the underlying model describing high temperature superconductors. For some of the systems described here the theoretical predictions can be compared to experimental measurements, as such measurements are already available; it would be interesting to see more such measurements preformed in the future for a more diverse range of systems. While not described here, a similar analysis has been performed also for the two-dimensional surface states of three dimensional topological insulators; a more detailed analysis of topological systems using these tools is foreseeable in the future. The Friedel oscillations are a spectacular and extremely useful tool to analyze and decode the physics of low-dimensional physics, and the STM measurements of these oscillations remain in my opinion one of the fundamental and most important tools of characterizing such a system.

References

Friedel J. Nuovo Cimento 1958 ;  7 : 287 [cross-ref]
Ruderman M.A., Kittel C. Phys. Rev. 1954 ;  96 : 99
Yosida K. Phys. Rev. 1957 ;  106 : 893
Kasuya T. Prog. Theor. Phys. 1956 ;  16 : 45
Knorr N., Brune H., Epple M., Hirstein A., Schneider M.A., Kern K. Phys. Rev. B 2002 ;  65 : 115420
Gautier F., Lenglart P. Phys. Rev. 1965 ;  139 : A705
Petroff F., Barthélemy A., Mosca D.H., Lottis D.K., Fert A., Schroeder P.A., Pratt W.P., Loloee R., Lequien S. Phys. Rev. B 1991 ;  44 : 5355
Bruno P., Chappert C. Phys. Rev. Lett. 1991 ;  67 : 1602
Hasegawa Y., Avouris Ph. Phys. Rev. Lett. 1993 ;  71 : 1071
Crommie M.F., Lutz C.P., Eigler D.M. Nature 1993 ;  363 : 524
Jeandupeux O., Bürgi L., Hirstein A., Brune H., Kern K. Phys. Rev. B 1999 ;  59 : 15926
Bürgi L., Brune H., Kern K. Phys. Rev. Lett. 2002 ;  89 : 176801
Pivetta M., Silly F., Patthey F., Pelz J.P., Schneider W.-D. Phys. Rev. B 2003 ;  67 : 193402
Sprunger P.T., Petersen L., Plummer E.W., Lægsgaard E., Besenbacher F. Science 1997 ;  275 : 1764
Petersen L., Hofmann Ph., Plummer E.W., Besenbacher F. J. Electron Spectrosc. Relat. Phenom. 2000 ;  109 : 97
Hoffmann Ph., Brinner B.G., Doering M., Rust H.P., Plummer E.W., Bradshaw A.M. Phys. Rev. Lett. 1997 ;  79 : 265
Hoffman J.E., McElroy K., Lee D.H., Lang K.M., Eisaki H., Uchida S., Davis J.C. Science 2002 ;  297 : 1148 [cross-ref]
McElroy K., Simmonds R.W., Hoffman J.E., Lee D.H., Orenstein J., Elsaki H., Uchida S., Davis J.C. Nature 2003 ;  422 : 592
Vershinin M., Misra S., Ono S., Abe Y., Ando Y., Yazdani A. Science 2004 ;  303 : 1995
Pascual J.I., Bihlmayer G., Koroteev Y.M., Rust H.P., Ceballos G., Hansmann M., Horn K., Chulkov E.V., Blügel S., Echenique P.M., Hofmann Ph. Phys. Rev. Lett. 2004 ;  93 : 196802
Vonau F., Aubel D., Gewinner G., Pirri C., Peruchetti J.C., Bolmont D., Simon L. Phys. Rev. B 2004 ;  69 : 081305(R)
Vonau F., Aubel D., Gewinner G., Peruchetti J.C., Bolmont D., Simon L. Phys. Rev. Lett. 2005 ;  95 : 176803
Mahan G. Many-Particle Physics  New York: Kluwer Academics/Plenum Publishers (2000). ISBN 0-306-46338-5 chapt. 4.
Bruus H., Flensberg K. Many-Body Quantum Theory in Condensed Matter Physics – An Introduction  Oxford: Oxford University Press (2005). ISBN 0-19-856633-6 chapt. 5, etc.
Byers J.M., Flatté M.E., Scalapino D.J. Phys. Rev. Lett. 1993 ;  71 : 3363
Salkola M.I., Balatsky A.V., Scalapino D.J. Phys. Rev. Lett. 1996 ;  77 : 1841
Ziegler W., Poilblanc D., Preuss R., Hanke W., Scalapino D.J. Phys. Rev. B 1996 ;  53 : 8704
Hirschfeld P.J., Wölfle P., Einzel D. Phys. Rev. B 1988 ;  37 : 83
Polkovnikov A., Sachdev S., Vojta M. Phys. Rev. Lett. 2001 ;  86 : 296
Podolsky D., Demler E., Damle K., Halperin B.I. Phys. Rev. B 2003 ;  67 : 094514
Wang Q.H., Lee D.H. Phys. Rev. B 2003 ;  67 : 020511
Bena C., Chakravarty S., Hu J., Nayak C. Phys. Rev. B 2004 ;  63 : 134517
Bena C., Kivelson S. Phys. Rev. B 2005 ;  72 : 125432
Bena C. Phys. Rev. Lett. 2008 ;  100 : 076601
Wehling T.O., Balatsky A.V., Katsnelson M.I., Lichtenstein A.I., Scharnberg K., Wiesendanger R. Phys. Rev. B 2007 ;  75 : 125425
Peres N.M.R., Klironomos F.D., Tsai S.W., Santos J.R., Lopes dos Santos J.M.B., Castro Neto A.H. Europhys. Lett. 2007 ;  80 : 67007 [cross-ref]
Peres N.M., Guinea F., Castro Neto A.H. Phys. Rev. B 2006 ;  73 : 125411
Vozmediano M.A.H., López-Sancho M.P., Stauber T., Guinea F. Phys. Rev. B 2005 ;  72 : 155121
Ando T. J. Phys. Soc. Jpn. 2006 ;  75 : 074716
Pogorelov Y.G.arXiv:condmat/0603327
Skrypnyk Y.V., Loktev V.M. Phys. Rev. B 2006 ;  73 : 241402(R)
Skrypnyk Y.V., Loktev V.M. Phys. Rev. B 2007 ;  75 : 245401
Katsnelson M.I., Geim A.K. Philos. Trans. R. Soc. A 2008 ;  366 : 195
Cheianov V.V., Fal'ko V.I. Phys. Rev. Lett. 2006 ;  97 : 226801
Mariani E., and al. arXiv:cond-mat/0702019
Kivelson S.A., and al. Rev. Mod. Phys. 2003 ;  75 : 1201
Mallet P., and al. Phys. Rev. B 2007 ;  76 : 041403(R)

Rutter G.M., and al. Science 2007 ;  317 : 219
Dupont-Ferrier E., and al. Europhys. Lett. 2005 ;  72 : 430
Fang A., Howald C., Kaneko N., Greven M., Kapitulnik A. Phys. Rev. B 2004 ;  70 : 214514
Van Mieghem P. Rev. Mod. Phys. 1992 ;  64 : 755

McCann E. Phys. Rev. B 2006 ;  74 : 245426
McCann E., Fal'ko V.I. Phys. Rev. Lett. 2006 ;  96 : 086805

Nilsson J., and al. Phys. Rev. B 2006 ;  73 : 214418
Mallet P., and al. C. R. Physique 2016 ;  17 (3–4) : 294-301in this issue.  [cross-ref]
Simon L., Bena C., Vonau F., Aubel D., Nasrallah H., Habar M., Peruchetti J.C. Eur. Phys. J. B 2009 ;  69 : 355

Mallet P., and al. Phys. Rev. Lett. 2008 ;  101 : 206802

Zhang Y., Brar V., Girit C., Zettl A., Crommie M.F.arXiv:0902.4793L. Simon, C. Bena, F. Vonau, D. Aubel, H. Nasrallah, M. Habar, J.C. Perruchetti, European Phys. Journal B, http://dx.doi.org/10.1140/epjb/e2009-00142-3.
Castro Neto A.H., Guinea F., Peres N.M.R., Novoselov K.S., Geim A.K. Rev. Mod. Phys. 2009 ;  81 : 109 [cross-ref]
Peres N.M.R., Guinea F., Castro Neto A.H. Phys. Rev. B 2006 ;  73 : 125411
Brey L., Fertig H. Phys. Rev. B 2006 ;  73 : 195408
Roldan R., Fuchs J.-N., Goerbig M.O.arXiv:0809.2667
Plochocka P., and al. Phys. Rev. Lett. 2008 ;  100 : 087401
Timusk T., Statt B. Rep. Prog. Phys. 1999 ;  62 : 61 [cross-ref]
Chakravarty S., Laughlin R.B., Morr D.K., Nayak C. Phys. Rev. B 2001 ;  63 : 4503
Schulz H.J. Phys. Rev. B 1989 ;  64 : 2940

Affleck I., Marston J.B. Phys. Rev. B 1988 ;  37 : 3774

Kotliar G. Phys. Rev. B 1988 ;  37 : 3664

Ivanov D.A., Lee P.A., Wen X.-G. Phys. Rev. Lett. 2000 ;  84 : 3958

Nayak C. Phys. Rev. B 2000 ;  62 : 4880
Howald C., Eisaki H., Kaneko N., Greven M., Kapitulnik A. Phys. Rev. B 2003 ;  67 : 014533
Flatté M.E., Byers J.M. Phys. Rev. Lett. 1998 ;  80 : 4546
Capriotti L., Scalapino D.J., Sedgewick R.D. Phys. Rev. B 2003 ;  68 : 014508
Podolski D., Demler E., Damle K., Halperin B.I. Phys. Rev. B 2003 ;  67 : 094514
Chen C.-T., Yeh N.-C. Phys. Rev. B 2003 ;  68 : 220505
Chen H.-D., Vafek O., Yazdani A., Zhang S.-C.arXiv:cond-mat/0402323
Fu H.C., Davis J.C., Lee D.-H.arXiv:cond-mat/0403001
Pereg-Barnea T., Franz M. Phys. Rev. B 2003 ;  68 : 180506
Misra S., and al. arXiv:cond-mat/0405204This paper uses the lowest order Born approximation, as opposed to the full T -matrix calculation. In addition, it ignores the anomalous part of the Green function when applied to the DDW. The results are therefore in disagreement with ours.
Norman M.R., Randeria M., Ding H., Campuzano J.C. Phys. Rev. B 1994 ;  52 : 615
Chakravarty S., Nayak C., Tewari S. Phys. Rev. B 2003 ;  68 : 100504(R)
Ghosal A., Kee H.-Y. Phys. Rev. B 2004 ;  69 : 224513
Tomonaga S. Prog. Theor. Phys. 1950 ;  5 : 544
Luttinger J.M. J. Math. Phys. 1963 ;  4 : 1154 [cross-ref]
Mattis D.C., Lieb E.H. J. Math. Phys. 1965 ;  6 : 304 [cross-ref]
Haldane F.D.M. J. Phys. C 1981 ;  14 : 2585
Gogolin A.O., Nersesyan A.A., Tsvelik A.M. Bosonization and Strongly Correlated Systems Cambridge: Cambridge University Press (1998). 
Glazman L.I., Ruzin I.M., Shklovskii B.I. Phys. Rev. B 1992 ;  45 : 8454
Kane C.L., Fisher M.P.A. Phys. Rev. B 1992 ;  46 : 15233
Kane C.L., Fisher M.P.A. Phys. Rev. Lett. 1992 ;  68 : 1220
Furusaki A., Nagaosa N. Phys. Rev. B 1993 ;  47 : 4631
Fisher M.P.A., Glazman L.I. Mesoscopic Electron Transport Dordrecht: Kluwer (1997).  331
Maslov D.L. Phys. Rev. B 1995 ;  52 : R14368
Fabrizio M., Gogolin A.O. Phys. Rev. Lett. 1997 ;  78 : 4527
Egger R., Grabert H. Phys. Rev. B 1998 ;  58 : 10761
Dolcini F., Grabert H., Safi I., Trauzettel B. Phys. Rev. Lett. 2003 ;  91 : 266402
Dolcini F., Trauzettel B., Safi I., Grabert H. Phys. Rev. B 2005 ;  71 : 165309
Oreg Y., Finkel'stein A.M. Phys. Rev. Lett. 1996 ;  76 : 4230
Furusaki A. Phys. Rev. B 1997 ;  56 : 9352
von Delft J., Schoeller H. Ann. Phys. 1998 ;  7 : 225
Grishin A., Yurkevich I.V., Lerner I.V. Phys. Rev. B 2004 ;  69 : 165108
Braunecker B., Bena C., Simon P. Phys. Rev. B 2012 ;  85 : 035136
Song X.L., Zhao Z.Y., Wang Y., Shi Y.M. J. Shanghai Univ. 2003 ;  7 (4) : 361
Eggert S., Johannesson H., Mattsson A. Phys. Rev. Lett. 1996 ;  76 : 1505
Anfuso F., Eggert S. Phys. Rev. B 2003 ;  68 : 241301R

1  The intensity is proportional to an integral of a Gaussian and two Hermite polynomials which ends up proportional to a Laguerre polynomial.


© 2015  Académie des sciences@@#104156@@
EM-CONSULTE.COM is registrered at the CNIL, déclaration n° 1286925.
As per the Law relating to information storage and personal integrity, you have the right to oppose (art 26 of that law), access (art 34 of that law) and rectify (art 36 of that law) your personal data. You may thus request that your data, should it be inaccurate, incomplete, unclear, outdated, not be used or stored, be corrected, clarified, updated or deleted.
Personal information regarding our website's visitors, including their identity, is confidential.
The owners of this website hereby guarantee to respect the legal confidentiality conditions, applicable in France, and not to disclose this data to third parties.
Close
Article Outline