text:
# Estimation of matrix element effects and determination of the Fermi surface in BSCCO systems using angle-scanned photoemission spectroscopy

## I Introduction

Since the beginning of the field of the high $`T_c`$ superconductors (HTSC), ARPES has taken a special role in the experimental study of these systems.  Among the numerous remarkable ARPES experiments on the HTSC, a special place belongs to the investigation of the Fermi surface (FS) of these systems. The difficulties encountered upon applying ’traditional’ techniques for determining the FS to the HTSC (such as de Haas - van Alphen and positron annihilation) focussed the attention on the possibility offered by ARPES to obtain a direct image of the basal-plane projection of the FS. The vast majority of the earlier ARPES work  agreed in the view that the FS in the most HTSC is hole-like, and centred at the X,Y points of the 2D Brillouin zone. Recently, however, a controversy regarding the FS topology of these systems has flared up. Some groups have suggested the presence of an electron-like, $`\mathrm{\Gamma }`$-centred FS in the Bi2Sr2CaCu2O8+δ (BSCCO) \[Refs. \] and Pb-BSCCO \[Ref. \] systems, which would represent a complete revision of our thinking regarding the fundamentals of the electronic structure of the HTSC. On the other hand, other groups have also re-visited this question  and confirmed the ’old’ picture of a hole-like FS in the Bi-cuprates. At the same time in La2-xSrxCuO4 a cross-over from hole-like to electron-like FS has been suggested from ARPES data on going to the overdoped side of the phase diagramme. Consequently, there still exists no consensus as to the correct picture for the normal state Fermi surface topology in the HTSC, making this question an important one to clarify.

The current debate as regards the Fermi surface topology is based to a large extent upon ARPES intensity maps. This means that the issue of matrix elements, which could be strongly photon energy and k-dependent in the 2D cuprate-based materials  has to be treated seriously. Thus we should be able to identify situations in which the matrix elements dominate and, where possible, develop practical methods of exctracting the underlying true information from the raw photoemission intensity signal. A second issue is that of how to accurately determine the Fermi momentum vectors, kF from real photoemission data. The accuracy within which this kF determination can be tested has increased dramatically in the past few years as a result of a new generation of electron energy analysers which offer resolutions in k-space one order of magnitude superior to what was previously available.

In this paper, we address the question as to how one can best locate Fermi momentum vectors in the HTSC with the aid of the angular distribution of photoemission intensity. As a case study we take BSCCO, but in fact the conclusions arrived at are quite general to the high resolution ARPES investigation of quasi-2D systems. We experimentally demonstrate the strong distortion of the ’pure’ photoelectron angular distribution caused by matrix element effects, thus making an appropriate further analysis of the raw data compulsory. To give such a statement a firm foundation, a detailed discussion of the experimental conditions and their influence on the photocurrent together with a critical overview of the existing methods of kF determination are presented. We then demonstrate that, using what we call a ’self-normalization’ procedure, one can significantly reduce the dependence of photoemission spectra on the matrix elements and finally show that this approach can be successfully applied to the BSCCO compounds. Consequently on the basis of the self-normalization method we formulate a criterion of determining the Fermi surface of the HTSC from ARPES data.

## II Methodology

Because of rapid evolution of the modern ARPES experiment and rising number of possible techniques for the solution of a given problem, we include in this section not only traditional experimental subsection (A) but also a quite detailed description of our approach to the FS mapping. In subsection B we discuss the quantities which are in principle accessible by ARPES in our implementation. Then we discuss how to optimize the experimental parameters for the study of the HTSC cuprates (subsection C) before closing this section with a critical evaluation of the available methods of kF determination (subsection D).

### A Experimental

Two types of experimental set-up have been used. The majority of the data discussed here were recorded with an overall energy resolution of 19 meV (FWHM) using a SCIENTA SES200 analyser coupled to a high intensity He resonance source (GAMMADATA VUV5000) via a toroidal grating monochromator (giving a degree of linear polarization of ca. 40%). The SES200 analyser provides an angular resolution of down to 0.2. The single crystals were mounted on a purpose-built, high precision cryo-manipulator which allows the sample to be rotated with a precision of better than 0.2 about three perpendicular axes in a wide range of angles. The synchrotron based data were recorded as described in Ref..

### B ARPES with analysis of the energy and momentum distributions on an equal footing

The new generation of electron energy analysers mentioned above have enabled a jump in angular resolution performance as a result of electron optics possessing an angular dispersion capability in the direction parallel to the analyser entrance slit. One can visualise the new mode of data collection with the help of Fig. 1. This shows the information landing on the analyser’s 2D-detector within 7 minutes’ measuring time in a static experiment - i.e. without moving the sample or the analyser. One direction on the detector represents an angular interval and the other direction - an energy interval (in this case $`\pm 7^{}`$ and $``$0.6eV respectively).

The dataset shown in Fig. 1 is a ’snapshot’ taken along the $`\mathrm{\Gamma }`$-$`(\pi ,\pi )`$ high symmetry direction of the Brillouin zone of Pb-doped BSCCO at 120 K, with the photoemission intensity plotted as a function of both binding energy and momentum. The left panel shows cuts of the intensity distribution $`I(𝐤,\omega )`$ parallel to the energy axis, i.e. energy distribution curves (EDC), which are uniquely defined by the fixed value of momentum. One can also cut the same $`I(𝐤,\omega )`$ distribution parallel to the momentum axis. These cuts are shown in the right panel and are termed momentum distribution curves or MDCs. An MDC should reflect the vector nature of the momentum, and thus is uniquely defined by the chosen frequency (binding energy) and an arbitrary path in a two-dimensional k||-space. It is clear from Fig. 1 that the modern ARPES machinery enables the simultaneous measurement of the energy and angular distribution of the photoelectrons leaving the sample and therefore permits a treatment of the energy and momentum dependence of the photocurrent on an equal footing.

The use of a grey scale to represent the intensity enables an efficient presentation of the three dimensional data set shown in Fig. 1. The result, shown in the upper panel of Fig. 2 is called an energy distribution map (EDM). An EDM can be thought of as an array of EDCs (or MDCs) taken along particular path in the BZ within particular range of binding energies. If we now fix binding energy and grey-scale-code the intensity in a series of MDCs covering an area in (kx, ky)-space together, then we arrive at a momentum distribution map or MDM, which represents a constant energy surface. The lower panel of Fig. 2 shows such an MDM for $`E`$=0 eV binding energy covering a part of the first Brillouin zone of Pb-BSCCO. In our case, the sample rotation involved in the recording of an MDM is such that the individual MDCs representing the angular breadth of the 2D detector are radial in nature, with the origin at the $`\mathrm{\Gamma }`$ point.

The informative capacity provided by such a map is easy to estimate, although not only the $`E_F`$-MDM is important. The dataset still possesses the binding energy axis, and so a series of MDMs then represents the evolution of the momentum distribution of the electronic states when going from the Fermi energy towards higher binding energies.

The inset to Fig. 2 summarizes the completeness of the information available in our ARPES experiment. Here we illustrate the three dimensional (kx, ky, $`\omega `$)-space which can be probed with high E and k resolution. The fourth dimension here is symbolised by the grey scale and represents the photoemission intensity. Recording the ARPES intensity while moving along any vertical direction, i.e. parallel to the energy axis  will give an EDC, whereas an MDC is the intensity distribution along the arbitrary path which belongs to any of horizontal planes in this space. Horizontal planes themselves are MDMs and the vertical surface (plane) defined by a given path (line) in an MDM is the EDM. Thus, the portion of (k, $`\omega `$)-space shown in the inset to Fig. 2 is confined by three EDMs (one is not visible) and two MDMs (one is not visible) and consists typically of approximately 100000 data points.

The foregoing discussion has illustrated the potential offered by our experimental setup. Nevertheless, despite significant advances there remain factors such as the resolution or lifetime of the sample surface which make it necessary to optimize the other experimental conditions for the treatment of the physical problem at hand. Moreover, given the widespread use of intensity maps in the literature regarding the FS topology of the HTSC, it is evident that the chosen experimental conditions (such as the experimental geometry, the excitation energy, the photon polarization, the temperature etc.) could decisively change the final picture obtained. Therefore, in the next sub-section we briefly deal with the different experimental parameters which could strongly influence the photoemission intensity distribution in (k, $`\omega `$)-space.

### C Factors influencing the $`I(𝐤,\omega )`$ distribution of photoelectrons

In order to discuss different parameters affecting the measured photoemission intensity, we first write an expression for the photocurrent as a function of (k, $`\omega `$)-space. For quasi-2D systems and under the assumptions that the ’sudden approximation’ applies  and that only a single initial state is involved, the photocurrent can be written in the following form

$$I(𝐤,\omega )=G_𝐤\{M(𝐤)[A(𝐤,\omega )f(\omega )]R_{\omega ,𝐤}+B(\omega )\},$$
(1)
where $`G_𝐤`$ is a mainly geometrical prefactor which will be described below, $`M`$ represents the square of the matrix element linking the initial and final states, $`A`$ is the single particle spectral function, $`f`$ is the Fermi function and $`R_{\omega ,𝐤}`$ is the energy and momentum resolution function. $`B`$ is the background, which contains extrinsic effects such as inelastic scattering of the photoelectrons. As an approximation we assume a negligible k-dependence of the extrinsic background and likewise a negligible $`\omega `$-dependence of both the matrix elements and prefactor $`G_𝐤`$ within the energy interval of interest ($``$ 0.3 eV). Eq. 1 makes it clear in a formal manner, that the measured signal is not simply the spectral function, and thus that a number of parameters must be known before $`A`$ can be extracted.

The pre-factor $`G_𝐤`$ describes the combined effects of extrinsic parameters such as occur upon the rotation of the sample with respect to the analyser (i.e. changing effective photon density in the area of the sample ’seen’ by the analyser) or the inequal efficiency of the different channels of the parallel detection system. The raw-data MDM shown in Fig. 2 (lower panel) illustrates the effect of $`G_𝐤`$, as it can be seen that at the interjoins of the two separate arcs of individual (radially arranged) MDCs there is an intensity misfit, mainly due to inequal detector channel efficiencies. In order to minimize the effects connected with $`G_𝐤`$, calibration scans can be carried out by measuring an isotropic photoemitter such as an amorphous gold film. Another, simpler way of overcoming this problem, is the self-normlaization which will be described later.

The generic step-like background $`B`$ observed in the ARPES of the HTSC is still a puzzle. The authors of Ref.  demonstrated that the contribution from secondary electrons, which could be estimated upon the basis of electron energy-loss spectroscopic data , is not sufficient to explain the background which is in correspondence with earlier assumptions in this regard. However, our assumption that the background is approximately k-independent and of practically negligible intensity at $`E_F`$ in comparison with the main signal is supported by the similar $`B(\omega )`$ line shape for k-points from unoccupied part of the BZ or for those k for which spectral function peaks at higher binding energies (e.g. close to the $`\mathrm{\Gamma }`$-point). Therefore, where needed, the background can be safely subtracted. We have found that as a good representative for the $`B(\omega )`$, an EDC from the vicinity of ($`\pi `$/2, $`\pi `$/2) or ($`\pi `$, $`\pi `$) points could be taken.

The most important component of Eq. 1 other than the spectral function itself is the $`M`$(k) term describing the matrix elements, which depends upon both the photon energy and the photoelectron momentum via the operator which couples the final and initial state wave functions. The choice of the energy of exciting photons is far from being unimportant even in the case of quasi-2D electronic systems. Firstly, upon changing the photon energy one alters the momentum resolution. Secondly, the 2D-CuO2 plane materials exhibit extremely strong variations in the ARPES intensity of their lowest lying ionisation states as a function of the photon energy. Thus, in extreme cases, by an unfortunate choice of the photon energy the contribution to the angular distribution of the photoelectrons from a particular initial state can be significantly suppressed.

The same goes for the ’angular’ part of the $`M`$(k) when using highly polarized radiation. If the experimental geometry can be controlled so as to give a clearly defined symmetry condition - such as can be the case along high symmetry lines in k-space - the strong polarization dependence of the photoemission signal can be used as a probe of the symmetry of the initial states involved. If, however, the strongly polarised radiation is used to measure ARPES spectra away from the high symmetry directions in the Brillouin zone - then the observed photoemission intensity represents only a part of the whole picture. In this context we note that it has been clearly demonstrated that ARPES intensity maps recorded from BSCCO using the same photon energy differ very significantly when recorded with differing polarization geometries.

Two ways around this problem spring to mind. Firstly, one could use unpolarized light. The laboratory He source and monochromator used here generates VUV radiation with ca. 40% linear polarization. Thus, the majority of the MDM intensity is coming from excitation with unpolarised light, meaning that although the polarized component will favour emission from particular states, the global effect is quite small and we are consequently able to ’see’ all the states involved. The second possibility is to use the variable polarization offered by modern insertion devices at synchrotron radiation sources to record intensity maps in pairs with complementary polarization geometry.

Moving further through the factors separating a real ARPES experiment from $`A`$(k, $`\omega `$) we come to the energy and momentum resolutions, represented in Eq. 1 by the function $`R_{\omega ,𝐤}`$. When using an angle-multiplexing analyser, momentum resolution can be projected onto two mutually perpendicular directions - parallel to the entrance slit of the analyser and perpendicular to the slit: $`R_{\omega ,𝐤}=R_\omega R_{k_{slit}}R_{k_{slit}}`$. The resolution parallel to the slit $`R_{k_{slit}}`$ is defined by the electron optical characteristics of the spectrometer whereas $`R_{k_{slit}}`$ is futher controlled by the aperture and entrance slit size. In most cases we used 19 meV $`\times 0.2^{}\times 0.5^{}`$ FWHM resolution.

Apart from the angular resolution of the analyser, the flatness of the sample surface as well as the excitation energy and the absolute values of the momenta define the momentum component of the (k, $`\omega `$)-resolution. For example, the k-resolution for Fermi level emission in experiments using ’high’ photon energies (e.g. 55 eV) is up to a factor 3 worse than while using typical ’low’ photon energies (e.g. 21.2 eV as here) for k-vectors in the second Brillouin zone of BSCCO for the same angular resolution.

The practical effect of the energy and momentum resolutions on the measured data depends strongly on the dispersion of the feature in (k, $`\omega `$)-space. The influence of each component increases as the direction of the most rapid change of intensity $`I(𝐤,\omega )`$ approaches the corresponding axis onto which the total resolution is projected. The direction of the most rapid intensity variation roughly coincides with the normal to the bare band (which is a surface in (k, $`\omega `$)-space). Fig. 3 illustrates three exemplary EDCs recorded from pure BSCCO at 40K. In each case at least one of the resolution contributions is zero for our experimental geometry. Fig. 3(a) shows the $`𝐤_F`$-EDC taken from the $`\mathrm{\Gamma }`$X cut, for which $`R_{k_{slit}}`$ is negligible - the width of the feature is defined jointly by $`R_\omega `$ and $`R_{k_{slit}}`$. In fact, for the strongly dispersing states along the nodal line ($``$ 2 eVÅ), the momentum resolution along the slit is the dominating factor. For the data of Fig. 3(a) the energy resolution was 19 meV FWHM and $`R_{k_{slit}}`$ was 0.015 Å-1 (i.e. $`0.2^{}`$). The latter causes an energetic broadening of some 30 meV, which easily outweights the contribution from $`R_\omega `$ itself. In choosing an EDC from the MX cut as shown in Fig. 3(b), we switch off the influence of $`R_{k_{slit}}`$. As in the $`\mathrm{\Gamma }`$X case, the momentum resolution (this time $`R_{k_{slit}}`$) still plays the leading role as regards the instrumental broadening of the observed peak. Fig. 3(c) illustrates an EDC from the M-point, where it is well-known that the sharp peak observed below Tc varies little in binding energy as a function of k,  meaning that the instrumental contribution to the width of this structure is determined effectively by the energy resolution. The three examples shown in Fig. 3 show that the conditions we apply result in a good balance between all three resolution components, when considered in the light of the typical Fermi velocities in the BSCCO-based materials.

Having discussed the influence of all the parameters entering into Eq. 1, we conclude that the most unpredictable and therefore difficult factor to deal with which separates the photocurrent from the spectral function are the matrix elements, which, as was mentioned above, are strong in quasi-2D CuO2-plane materials.

We now turn our attention to a discussion of the best manner in which kF vectors can be derived from the ARPES data of the HTSC. This point is much more than merely a detail of the ARPES data evaluation, and lies at the heart of the current Fermi surface controversy. In particular, given the conclusion that it is only the matrix elements that severely hamper a clear view of the spectral function, our discussion of how to determine kF is thus centered on the question as to which level each method is immune, if at all, from matrix element effects.

### D What is the best way to determine when k=kF?

#### 1 Dispersion method (EDC maximum)

The relatively coarse k-mesh available in the majority of earlier ARPES investigations meant that FS crossings could only be located by the analysis of a series of EDCs (i.e. an EDM) containing a dispersive feature. The most simple and intuitively direct method is then to follow the energy position of the EDC maximum and to extrapolate obtained dispersion relation to E = $`E_F`$. This procedure, however, suffers from a number of drawbacks. Firstly, the influence of the Fermi cut-off distorts the picture within ca. 2kT of $`E_F`$. Secondly, if the self-energy is frequency dependent, following the EDC maxima will lead to the wrong result. This point can be simply visualised with the help of the three-dimensional plot, shown in the central panel of Fig. 1. As the intensity varies along the bare band, the trace formed by joining the maxima of cuts through this object taken parallel to the binding energy axis (EDCs) can never agree with that obtained by joining the maxima of the cuts taken parallel to the momentum axis (MDCs). If the self-energy only weakly depends on k, it is evident that the MDC dispersion is much closer to the ’true’ dispersion (i.e. the bare one plus the real part of the self-energy) than that from the EDC maxima. As a consequence, even though it does possess the advantages that it is insensitive to the normalization procedure, to the effects of finite energy resolution and is quite robust with respect to matrix elements effects, the EDC method gives only approximate values of kF.

#### 2 $`\mathrm{\Delta }`$T method

It has been proposed from ARPES measurements of TiTe2 and from simulations  that $`E_F`$-MDCs shift as a function of temperature in such a way that the difference of such MDCs turns out to be zero only for k=kF. The practical application of this method to the HTSC is blocked by two points. Firstly, everything has to be measured twice (for T1 and T2) with very high k-space location precision, with all other parameters being kept equal. This is often impossible due to the finite lifetime of the cleaved surfaces of the HTSC. Secondly, and more fundamentally, the $`\mathrm{\Delta }`$T method cannot function if the width of the EF-MDCs concerned is temperature dependent. As this is very clearly the case in the HTSC, the $`\mathrm{\Delta }`$T method is invalid in the context of the HTSC. Without wishing to persue this point further here, we refer the reader to the Appendix for the analytical evidence which forms the basis of these statements.

#### 3 Symmetrization

The symmetrization method is based upon an analysis of the lineshape of the result obtained by mirroring the photoemission EDC’s around E=$`E_F`$ and summing up the two spectra for each k-point, a procedure that can be described as $`I_S(𝐤,\omega )=I(𝐤,\omega )+I(𝐤,\omega )`$. Within this method kF is defined as the point at which the dip at EF in the symmetrized EDC’s (k$`<`$kF) turns to maximum. It is evident, however, that upon approaching kF the two peaks originated from the $`\omega `$ and -$`\omega `$ spectra will approach each other and become indistinguishable, giving a maximum in the symmetrized EDC over a range of k. One can estimate this k-range quantitatively by solving the equation $`(d^2I_S(𝐤,\omega )/d\omega ^2)_{\omega =0}=0`$. For the model spectral function considered in the Appendix, the solution gives k-kF $``$ 0.02 Å-1 at 300K. Consequently, the k point at which the dip in the symmetrized EDC’s transforms into a peak is shifted away from the true kF introducing considerably larger error than, for instance, the MDC maximum method gives (see below). Symmetrization does possess the advantage that it eliminates the Fermi cut-off from the kF-EDC and answers the question whether a FS crossing occured or not in a given EDM even with the presence of strong matrix element effects.

#### 4 Maximal gradient of the integrated intensity

It is well known that even for interacting Fermi systems, kF is characterised by a jump in the momentum distribution, $`n(𝐤)`$. For finite temperatures one could still, in principle, detect rapid variations in $`n(𝐤)`$ and estimate kF from the $`\mathrm{max}\left|_𝐤n(𝐤)\right|`$. It has been proposed  that the integrated intensity ($`I_{int}`$) of an EDC could give a measure of $`n(𝐤)`$ at one particular k-point and thus the analysis of a series of EDC’s could present an opportunity to estimate kF. This method has been applied to different systems, including BSCCO, and is still intensively used.

There exist, however, the following arguments against the $`\left|I_{int}(𝐤)\right|`$ method. Firstly, $`I_{int}(𝐤)`$ is not equal to $`n(𝐤)`$, for the same reasons that the raw photoemission intensity is not equal to the spectral function. Secondly, a single EDC does not represent the photoemission intensity for a single k-point, but rather for a range of k-points. For lower photon energies (e.g. 21.2 eV), the finite k-interval associated with the finite energy width of an EDC can even be comparable with the momentum resolution. Thirdly, the band structure of the system has to be amenable, in the sense that a single state needs to be isolated and well away in frequency from other features. Even if this is the case, the intensity integration should be carried out over all frequencies - in practice however, energy windows varying between 100-600 meV in width are taken for the integration. Taking a narrower window reduces the similarity with $`n(𝐤)`$, whereas a wider window results in enlargement of the k-interval and increases the contribution from deeper lying valence band states. A further difficulty arises from the data analysis in that numerical differentiation introduces additional errors. Such transformation of the raw data also produces a set of additional ”false” features on the map  which have to be identified as such and neglected at a later stage. In any case, the quantitative precision of the kF determination is in direct relation to the width of $`\left|I_{int}(𝐤)\right|`$. This width is much broader than, for example, a typical $`E_F`$-MDC for the 2 eVÅ-like dispersive features in the BSCCO compounds. The factors mentioned above make the $`\left|I_{int}(𝐤)\right|`$ method intrinsically inaccurate. In the Appendix, we use simulations to show further that this method can result in substantial systematic errors determining kF, which is in agreement with the results of other authors.

#### 5 Maximum intensity method (MDC maximum)

The maximum intensity method, as first introduced for BSCCO in Ref.  is based on measuring the $`E_F`$-MDM. In this case, the photocurrent is recorded only in a narrow energy window centred on the chemical potential, thus enabling the coverage of large areas in momentum space within a relatively short time.

The physical basis of this method of kF determination is straightforward. One starts with the reasonable assumption that at finite temperatures the spectral function (independent of the model used to describe it) has a peak at $`\omega `$=0 only for k=kF. It then immediately follows from Eq. 1 that every $`E_F`$-MDC corresponding to a path in an $`E_F`$-MDM which crosses the Fermi surface will show a maximum. This property of $`E_F`$-MDCs has been recognized and successfully applied for kF determination by a number of groups. For the MDC peak to lie exactly at kF, the influence of the $`G_k`$-factor, the matrix elements ($`M(𝐤)`$) and resolution ($`R_{\omega ,𝐤}`$) should not be strong enough to shift the peak position of the spectral function. In this context, the essentially symmetric resolution functions can certainly be regarded as harmless.

#### 6 Influence of matrix element effects on kF determination

It is much more tricky to evaluate how strong the dependence on the matrix elements is for a particular set of experimental conditions. Numerical calculations of the photoemission intensity including the matrix elements have been carried out, and predict that the matrix elements should have a dramatic effect on the angular distribution of photoelectrons that would be detected in an ARPES experiment, given an identical underlying Fermi surface topology.

In Fig. 4(a) we show $`E_F`$-MDCs and in Fig. 4(b) the integrated ARPES intensity $`I_{int}`$ both recorded for a path in k-space along the $`\mathrm{\Gamma }`$-M-Z direction in BSCCO compounds. In each case the four panels show data measured using different photon energies, whereby the 21 eV panel was recorded using radiation from a He resonance source.

Comparing the shape of the MDCs and $`I_{int}`$ traces at different excitation energies it is clear that a kF determination method which involves the photoemission intensity such as the maximal gradient of integrated intensity or the MDC maximum methods may function poorly in such a case. The vertical dashed lines marking the apparent location of Fermi vectors show different results between the two methods, and, more importantly, different results for the same method, depending on the experimental conditions. This holds even for the maximum-MDC method.

This strong dependence of the apparent location of kF on the experimental conditions used to measure the ARPES data is, in fact, at the root of the current controversy regarding the Fermi surface topology of the HTSC. At this stage one could even be led to doubt the value of ARPES as a method of determining the Fermi surface in the HTSC. Obviously, there is an urgent need to find a method which is able to accurately reflect the Fermi surface even in the presence of strong matrix element effects. In the following, final results section, we describe an approach which allows not only an estimation of the distorsions caused by matrix elements but also enables a robust and precise determination of kF vectors.

## III Fermi surface mapping

From the discussion of the existing methods of kF determination it follows that, in the absence of strong matrix element effects, the most precise, simple and physically transparent approach is the MDC maximum method. Therefore one natural way forward is to improve this method by minimizing its sensitivity to the matrix elements. The form of Eq. 1 suggests the possibility of being able to ’divide out’ both the $`G_𝐤`$ and $`M(𝐤)`$ prefactors, providing the denominator is proportional to their product. Within the energy range under consideration, one can take both $`M(𝐤)`$ and $`G_𝐤`$ to be frequency independent. Thus, a perfect candidate for such as division would be a signal from the same k-point, i.e. from the same EDC. At the same time, the reference signal should be a slowly varying function of momentum in the vicinity of the expected MDM (MDC) maxima. In principle, we should restrict ourselves to a rather small energy interval, so as give a narrow kF-window for each EDC.

Although several possibilities exist for such as matrix-element elimination, in the following we concentrate on the division by the integrated intensity, as this has already been successfully applied to the determination of the FS topology in BSCCO systems.

As discussed above in the context of the $`\left|I_{int}(𝐤)\right|`$ method, the integrated intensity versus k proves to be quite a slowly varying function in BSCCO in the vicinity of the expected Fermi surface crossings in comparison with the k-dependence of the Fermi energy intensity (i.e. the $`E_F`$-MDM). Within the framework of Eq. 1, the intensity after division, $`I_{norm}`$ is given by:

$`I_{norm}`$ $`=`$ $`{\displaystyle \frac{I(𝐤,0)}{_ϵI(𝐤,\omega )𝑑\omega }}=`$ (2)
$`=`$ $`{\displaystyle \frac{G_𝐤\{M(𝐤)[A(𝐤,\omega )f(\omega )]R_{\omega ,𝐤}|_0+B(0)\}}{G_𝐤\{_ϵM(𝐤)[A(𝐤,\omega )f(\omega )]R_{\omega ,𝐤}𝑑\omega +_ϵB(\omega )𝑑\omega \}}},`$ (3)
where $`ϵ`$ is the energy window of integration, normally chosen between 600 meV and -100 meV. It is easy to see that prefactor $`G_𝐤`$ cancels out immediately without any additional assumptions, thus automatically solving the problem of the detector efficiency calibration.

As this method involving the division by the integrated intensity has been criticised as being ’unphysical’, we now consider the behaviour of the function described in Eq. 2 above in different parts of the Brillouin zone in detail. We consider three regions in the zone: ’definitely occupied’ ($`E_𝐤>\omega _{max}`$, where $`E_𝐤`$ is the quasiparticle dispersion and $`\omega _{max}`$ is the higher binding energy of the integration window), ’definitely unoccupied’ ($`E_𝐤<0`$) and ’close-to-the-Fermi-surface’ ($`E_𝐤0`$) regions, and discuss the contribution to signal at $`E_F`$ from the spectral function and extrinsic background. The first region contains k-values for which the spectral function peaks more than 0.5 eV from the chemical potential. As regards the ’definitely unoccupied’ region, the Fermi cutoff means that the signal at $`E_F`$ from the spectral function disappears faster going into the unoccupied part of the Brillouin zone, and thus we consider the ’definitely unoccupied’ part as being located immediately after the FS crossing. In both of these regions, the $`I_{norm}`$ signal is determined only by the background lineshape. The weak k-dependence of the background means that $`I_{norm}`$ should be approximately constant in the ’definitely unoccupied’ and ’definitely occupied’ regions, although it is somewhat noisy due to the small values of $`B(0)`$. In the immediate vicinity of the expected Fermi surface, we can neglect both terms connected with the background, as $`B(0)`$ is typically an order of magnitude lower than I(kF,$`E_F`$) (as can be estimated from, for example, the data of Fig. 1) and the corresponding ratio of the integrated intensities is $``$ 0.15. In this case, the $`I_{norm}`$ function can be rewritten as:

$`I_{norm}\{\begin{array}{cc}\frac{B(0)}{_ϵB(\omega )𝑑\omega }const\hfill & \text{if }E_𝐤>\omega _{max}\text{ or }E_𝐤<0;\hfill \\ \frac{[A(𝐤,\omega )f(\omega )]R_{\omega ,𝐤}|_0}{_ϵ[A(𝐤,\omega )f(\omega )]R_{\omega ,𝐤}𝑑\omega }\hfill & \text{if }E_𝐤0.\hfill \end{array}`$ (6)
These relations show that, in the vicinity of kF, the function is independent of both the matrix elements and the $`G_𝐤`$-factor. It should be noted, however, that in the case of severe suppression of the intensity related to the spectral function, the background contribution becomes substantial and therefore $`I_{norm}`$ cannot be considered as matrix element independent - even in a first approximation - but still can be used for identifying such a situation.

Although not immediately obvious, the physical meaning of the $`I_{norm}`$ function is quite transparent. Consider its behaviour along a single cut through the 2D BZ (e.g. that shown in Fig. 1). The numerator is simply an $`E_F`$-MDC and the denominator is, in this case, the $`I_{int}`$ along $`\mathrm{\Gamma }`$X. We expect the integrated intensity to show a slow drop in the region of kF. At the same k-region the true, underlying MDC has a sharp maximum exactly at kF. Thus, the division of the narrow Lorenzian-like MDC function by the slowly falling $`I_{int}`$ function does not even result in a significant shift of the maximum, and makes itself felt only in the asymmetric shape of the renormalized MDC. In other words, by normalizing the MDC in such a manner, we do not change the position of its maximum, thus preserving its most important property as an indicator of kF. We note that the asymmetry of the normalised MDC peak can act as an additional guide in determining whether the FS has been crossed from the occupied part to the unoccupied or in the opposite direction.

Having now discussed the normalised MDC approach on paper, as it were, we now put it to the test. As was evident in Fig. 4, the determination of whether there is a kF vector along the $`\mathrm{\Gamma }`$-M-Z direction in pure Bi2212 represents a quite a challenge. In Fig. 5 we show the result of dividing each of the panels of Fig. 4(a) with the corresponding panel of Fig. 4(b) - i.e. Fig. 5 contains the renormalised $`E_F`$-MDCs obtained by division of the raw MDCs by the $`I_{int}`$ curves.

The major advance here over Fig. 4 is that now the normalized MDCs are all very similar, which, of course, is logical as there can be only one spectral function for this system. The nearly perfect coincidence between data recorded with a wide range of excitation energies is remarkable, and leaves no doubt that it was the matrix elements which had led to the differences in the raw MDC’s seen in Fig. 4. This example shows that in this way we can get an estimate of the influence of the most unpredictable part of the k-dependent matrix elements - i.e. that part which is not due to the symmetry selection rules.

Comparison of the raw \[Fig. 4(a)\] and normalized \[Fig. 5\] MDCs indicates that the strongest photon energy dependent variations in the raw intensity, $`I(𝐤)`$, are in the vicinity of the M-point. For $`h\nu =21.2eV`$, $`M(𝐤)`$ is asymmetric with respect to M-point, whereby emission from the states in the first BZ is favoured. The same asymmetry holds for $`h\nu =32eV`$, which is in agreement with other experiments. The most noticeable difference for 32 eV photons is the strong suppression of the spectral weight around the M-point itself - i.e. the raw MDC has a minimum at M. There is less pronounced M-point suppression for $`h\nu =40eV`$, and for $`h\nu =50eV`$ the suppression is no longer observed and the asymmetry has changed to favour emission from the 2nd BZ. This last point is a forerunner of the severe matrix elements effects at $`h\nu =55eV`$ where intensity from the second zone is completely dominant. Thus, as regards the spectral weight suppression, the k-interval corresponding to that part of the MDC which is ’eaten away’ by matrix element effects moves from the second BZ to the first when the excitation energy is varied from 21 to 55 eV, thus resulting in highest intensity suppresion at the M point for $`h\nu `$ ca. 32 eV, as had been predicted from photoemission calculations.

The remaining differences between the self-normalized MDCs presented in Fig. 5 could be due to a number of factors. For example, the M-point spectral weight suppression can be drastic enough (see, for example, Fig.1(a) in Ref.  for $`h\nu =32eV`$) to mean that one is no longer able to distinguish between the real signal and the background. Three of the curves presented in Fig. 5 are from pristine Bi2212, whose $`E_F`$-MDM is highly complex in the vicinity of the M point. Thus some of the observed fine structure could be due to the ’FS crossing’ of the two diffraction replicas.

Here we wish to stress that using the self-normalization method we are able to get an impression of the effects of the matrix elements. To determine the matrix element itself in a rigorous way, we would need to know the spectral function. In this sense, the self-normalized spectra shown in Fig. 5 are not equal to $`A(𝐤,\omega )`$ across the whole range of $`𝐤\omega `$-space but are very close to the spectral function for k near kF. It is this property of self-normalized MDC’s which makes them so suitable for FS mapping.

As a final comment to Fig. 5 we note that the universality of the lineshape of the normalized MDC’s along $`\mathrm{\Gamma }`$MZ in the BSCCO-based HTSC - in which no sharp peaks characteristic of FS crossings occur either side of the M-point - strongly suggests the absence of main band FS crossings in this direction in the BZ and thereby lends weight to the hole-like FS side in the Fermi surface controversy.

The foregoing discussion was limited to the self-normalization procedure using the integrated intensity as the normalising quantity. In Fig. 6 we show both raw and normalized MDMs recorded at room temperature within 3 hours after the cleavage of a Pb-BSCCO (OD 72K) single crystal.

Firstly, dealing with Fig. 6(b), the self-normalized MDM confirms the behaviour expected for the function given in Eq. 6. Both the ’definitely occupied’ and ’definitely unoccupied’ regions have approximately the same intensity (mid to dark grey \[red\] tone). On approaching the Fermi surface from the occupied side, the signal is reduced (giving dark areas) as here the contribution from the signal at the $`E_F`$ is still small but the integrated intensity is already quite large. There are two futher possibilities: if a band crosses the Fermi level (e.g. along $`\mathrm{\Gamma }`$-X) a sharp increase in intensity is observed giving a bright feature on the map. Alternatively, in the case of the $`\mathrm{\Gamma }`$-M-Z cut (see Fig. 5), the MDC has a plateau near the M-point, reflecting the behaviour of the flat band which does not cross the FS but approaches near enough to it to contribute to the signal at $`E_F`$. This leads to a fairly uniform intensity (mid grey \[red\] tone) around the M-point. Upon careful analysis of the locations of the MDC maxima in both raw and self-normalised maps, we find that there is no detectable shift between the two datasets, which confirms again that the $`I_{int}`$ function varies much more slowly than the MDC does in the vicinity of the FS. We now turn to Fig. 6(c), showing the self-normalization result using the intensity at high binding energies (in this case 250 meV below EF) as the denominator. The self-nomalised MDM in Fig. 6(c) displays all the characteristics of the map shown in 6(b), indicating firstly that the high binding energy signal is also sensitive enough to the matrix element effects to enable their elimination for k near to kF. Secondly, this points to the soundness of the assumption that the matrix element effects are insensitive to energy on the range of 0.5 eV. Thirdly, the good agreement between the self-normalization based upon the integrated intensity or high binding energy denominators proves the physically sound basis of the former procedure, in contrast to what is claimed in Ref. . Finally, we return to the BSCCO FS controversy, and point out that Fig. 6 shows without a doubt that the hole-like Fermi surface topology, which can be the only conclusion upon looking at our FS maps, is not a product of the integrated intensity normalization procedure, but is a robust result.

To summarise this section regarding the use of self-normalization to reduce the strong matrix element effects in ARPES of the HTSC we can say: (i) self-normalization conserves all the advantages of maximum intensity method of kF determination and comes close to the ideal of a method which delivers robust, precise results even in the presence of strong matrix element effects (ii) the denominator used in the self-normalization should be a signal which ’feels’ the matrix element effects, but which varies relatively slowly in k-space . We have shown both the integrated intensity and the intensity at higher binding energies to be two candidates which function well. (iii) the effectiveness of this method means that one can overcome the doubts raised earlier as regards the power of ARPES to determine the Fermi surface: the self-normalized MDM represents directly the Fermi surface map.

The self-normalization method also implies a formal criterion for the FS determination: a k-point in the 2D BZ belongs to the Fermi surface when and only when all possible normalized MDCs in its vicinity, except may be one (to account for inequal intensity distribution along the FS itself), have a local maximum at that point. To take into account the finite resolution of the experiment one would have need to introduce a tolerance angle within which the MDCs could be considered as belonging to the FS.

As regards the extension of the self-normalization method to quasi-2D systems other than BSCCO, we see no limitations in terms of its validity since the method is based on very general grounds. Even in the case of very low Fermi velocities the self-normalization effectively compensates already observable shift of the $`E_F`$-MDC due to the finite energy resolution, i.e. gives exact locations of kF-vectors.

## IV Summary

In this paper we have discussed the factors which separate the data of a real ARPES experiment from the spectral function, which is a highly topical subject in the light of the current controversy regarding the ARPES-derived Fermi surface topology in the HTSC. Based on high quality ARPES data of the BSCCO and Pb-BSCCO systems recorded under a variety of experimental conditions (photon energies, degrees of polarization), we have suggested a simple method which enables an estimation of the strength of the matrix elements effects in the ARPES of the HTSC, which at the same time allows the precise determination of the Fermi surface even when the matrix elements are strongly k-dependent. In this approach, a self-normalization effectively immunises the momentum distribution curves - whose maxima deliver precisely the Fermi wavevectors - against the matrix elements and extrinsic factors separating the photoemission signal from the spectral function for k near to kF. Consequently, the self-normalized momentum distribution map of the photoemission intensity at $`E_F`$ gives the most faithful reproduction of the underlying Fermi surface topology achievable from real ARPES data, and thus provides easy access to the quantitative analysis of the Fermi surface in these materials.

## V Acknowledgments

This work was funded in part by the BMBF (05 SB8BDA 6), the DFG (Graduiertenkolleg ’Struktur- und Korrelationseffekte in Festkörpern’ der TU-Dresden) and the SMWK (4-7531.50-040-823-99/6). H.B. is supported by the Swiss National Science Foundation and EPFL. We are grateful to G. Reichardt (BESSY GmbH), R. Müller and Ch. Janowitz (Humboldt Universität Berlin) for assistance and to Mike Norman for correcting the error in the calculations.

## VI Appendix

In this section we communicate in a little more detail certain points as regards the question of how best to determine kF. To enable a more quantitative analysis of the strengths and weaknesses of the various kF methods, we have generated a simulated dataset, based upon a fit to a real $`\mathrm{\Gamma }`$X EDM dataset from BSCCO. Starting from equation (1) given earlier and assuming that the problems of matrix elements, detector calibration and background have been adequately dealt with, we adopt a model in which the photocurrent can be calculated as

$$I(k,\omega )[A^{}(k,\omega ,R_k)f(\omega )]R_\omega .$$
(7)
To speed up the calculations, we combine the spectral function with the momentum resolution $`R_k`$ in

$$A^{}(k,\omega ,R_k)\frac{\sqrt{\mathrm{\Sigma }^{\prime \prime 2}+R_k^2}}{(\omega ϵ_k)^2+\mathrm{\Sigma }^{\prime \prime 2}+R_k^2}$$
(8)
The absence of a strong asymmetry in the $`\mathrm{\Gamma }`$X MDC’s results in a straightforward influence of the momentum resolution, which is in contrast to the influence of the energy resolution,  which is taken into account via convolution with the resolution function $`R_\omega (\omega )=(R_\omega \sqrt{\pi })^1exp(\omega ^2/R_\omega ^2)`$. For the imaginary part of the self-energy we use the following approximation: $`\mathrm{\Sigma }^{\prime \prime }(\omega ,T)=\sqrt{(\alpha \omega )^2+(\beta T)^2}`$ with $`\alpha =1`$ and $`\beta =2`$ ($`\omega `$ and $`T`$ in energy units), which, as can be seen from Fig. 7, gives a resonable fit to the experimental data.

Fig. 7 shows typical experimental Pb-BSCCO $`\mathrm{\Gamma }`$X EDM’s (for 30 K and 300 K) as contour plots (left panels), together with the results of the simulation (right panels). The quasiparticle dispersion $`ϵ_k`$ includes the effect of the real part of self-energy $`\mathrm{\Sigma }^{}(\omega ,T)`$, but in the region of interest near to the Fermi level we consider $`ϵ_k=v_Fk`$, where $`v_F`$ is simply the renormalized Fermi velocity at $`\omega =0`$. We took $`v_F`$= 2 eVÅ from the experimental data. Having described the basis of our simulations, in the following we analyse how accurate the kF determined by the different methods is.

### A Maximum MDC method

To evaluate the precision of this method quantitatively we simulate the $`\mathrm{\Gamma }`$X $`E_F`$-MDC according to (7) and (8). The results of this simulation are presented in Fig. 8.

The first observation is that the error in determining kF is very small using this method. For example, for room temperature the error is less than 0.001Å-1 for an energy resolution of 19 meV as was used in the experiment. Even for a resolution of 50 meV, the error is maximally 0.007Å-1. Secondly, the shift of the ’observed’ kF from the true value is only weakly temperature dependent, which therefore cannot be considered as an obstacle to the use of the ’maximum MDC’ method.

### B $`\mathrm{\Delta }T`$ method

The original $`\mathrm{\Delta }T`$ proposal  is based upon there being a temperature dependence of the position of the MDC maxima. Thus, the applicability or otherwise of the $`\mathrm{\Delta }T`$ method to the HTSC can also be judged from Fig. 8. The first point is that, as discussed above, the T-dependent shifts of the MDC maxima are very small in BSCCO, in contrast to the case in TiTe2. Furthermore, Fig. 8(a) also shows clearly that there is no common crossing point on the right flank of the $`E_F`$-MDC’s, which is a result of the temperature dependence of the width of the MDC’s, thus making the $`\mathrm{\Delta }`$T method inapplicable. Finally, even if one assumes a temperature-independent width of the $`E_F`$-MDC, the accuracy of the $`\mathrm{\Delta }T`$ method, $`\delta k`$, is related to the uncertainty in the determination of the relative intensities of each of the MDC pairs, $`\delta I`$: $`\delta k\delta I/(dI/dk)_{k=k_F}`$. In our case, to reach an accuracy of $`3\times 10^3`$ Å-1, the $`\delta I`$ would have to be less than 3% of $`I`$, which is beyond most present experimental capabilities.

### C Gradient $`I_{int}(k)`$

Fig. 9(a) shows the results of the simulation as regards gradient $`n(k)`$, in which

$$\frac{dI_{int}(k)}{dk}_{\omega _{min}}^{\omega _{max}}\left[\frac{dA^{}}{dk}f(\omega )\right]R_\omega (\omega )𝑑\omega $$
(9)
for $`\omega _{min}=0.1`$ eV, $`\omega _{max}`$ = 0.6 eV is plotted for four different temperatures. None of the maxima $`dI_{int}(k)/dk`$ lie on the $`\mathrm{\Delta }`$k=0 line, indicating a systematic error in the determination of kF. Fig. 9(b) shows the temperature dependence of this shift away from the true kF, $`\mathrm{\Delta }k`$, plotted for different $`R_k`$ = (10, 30, 60, 100) meV which are equivalent to (3, 10, 20, 33)$`\times 10^3`$ Å-1 or (0.09, 0.27, 0.54, 0.90) of angular resolution. Since our currently best intrumental angular resolution is 0.2, we discuss the curve for an angular contribution of 30 meV. Here the error at low temperatures is between 0.002 and 0.003 Å-1, which is as good as the maximum MDC method at these temperatures. For higher temperatures (e.g. for T$`>`$T\* in the HTSC for which the Fermi surface is not gapped), the error from gradient $`I_{int}(k)`$ has risen to 0.008 Å-1, some 8 times higher than the corresponding value for the maximum-MDC method.