text:
# Structure of Bi nanolines: using tight-binding to search parameter space

## I Introduction

When confronted with a new structure observed in experiment, or when trying to find possible pathways for reactions or diffusion events, there is often a need for performing electronic structure calculations on many different structures rather quickly, in order to search through parameter space efficiently. However, the need is not just for speed, but also for reasonable accuracy, as possible candidates must not be discarded out of hand. In this paper, we explore the use of tight binding as a tool to perform a search for the structure of Bi nanolines on Si(001), and consider future directions for these types of search using techniques developed over the last few years.

If a Bi-covered Si(001) surface is annealed at 570-600C, most of the Bi desorbs, but a small amount self-assembles to form nanolines, 1.5 nm wide and often over 400 nm longK.Miki et al. (1999a, b). These nanolines have great promise as templates for nanowiresJ.H.G.Owen et al. (2002a), when metals such as Al or Ga are deposited on them. However, to understand and control these processes and the growth of the nanolines, it is important to know their atomic structure; earlier proposals for the structureK.Miki et al. (1999a); D.R.Bowler (2000); M.Naitoh et al. (2000) have been shown to be wrong by recent experimental results, and a valid structure has only recently been found using a combination of STM and electronic structure modellingJ.H.G.Owen et al. (2002b). We describe the search for the structure, and the tools used for that search, in this paper.

Tight bindingC.M.Goringe et al. (1997) is a compromise between accuracy (it makes approximations, but preserves quantum mechanics) and computational efficiency (simulations of several hundred atoms can be run on a typical desktop workstation in a few hours). On semiconductor surfaces, where a rich variety of reconstructions exist, depending on a subtle interplay of electronic and geometric effects, the use of a quantum mechanical technique is vital if correct answers are to be found. We have developed parameterisations for Si(001), H on Si(001)D.R.Bowler et al. (1998a) and Bi on Si(001)D.R.Bowler (2000), and performed simulations on hydrogen diffusionJ.H.G.Owen et al. (1996); D.R.Bowler et al. (1998b, 2000a), gas-source growth of Si(001)J.H.G.Owen et al. (1997a, b, c); D.R.Bowler and M.G.Bowler (1998); D.R.Bowler and C.M.Goringe (1998) and Bi on Si(001)K.Miki et al. (1999a, b); D.R.Bowler (2000); J.H.G.Owen et al. (2002b) and have found that tight binding produces results which are in good agreement with both experiment and first principles techniques. In this paper, we further explore the use of tight binding as a tool to search through parameter space.

The problem we were faced with was that of finding the structure of the Bi nanolines. Various criteria had been found from experiment which possible candidates needed to fulfil: the lines must have a high kinking energy; the lines must occupy four Si dimers but lie between the substrate dimers; they must be more stable than the Bi-(2$`\times `$n) high coverage phase; and they must repel missing dimer defects and steps down. We needed an accurate, quantum mechanical technique that would allow us to test many possible structures quickly, allow us to have confidence in new reconstructions and distorted bonds, and that would also allow us to model large unit cells (to test kinking and defect repulsion) – a linear scaling or $`𝒪(N)`$ tight binding technique was found to be perfect. The rest of the paper describes the search: we start by describing the technical details of the modelling; we then consider the problem of modelling Bi on Si(001) and comparing structures with different amounts of Bi, and look at the stability of the Bi-(2$`\times `$n) high coverage phase; we present the details of various more or less successful models, and compare their energies to first principles calculations, and then we conclude.

## II Technical Details

We performed tight binding simulations using a linear scaling technique (where the computational effort scales linearly with the number of atoms, not with its cube as is true for traditional techniquesS.Goedecker (1999)) as large cells, of up to several thousand atoms, were necessary to test the kinking and dimer repulsion properties of the nanoline structures. Specifically, we used an implementation of the Density Matrix Method (DMM)X.-P.Li et al. (1993) by GoringeC.M.Goringe (1995). We used tight binding parameterisations for Si-Si bonds developed specifically for Si(001)D.R.Bowler et al. (1998a) and for Bi on Si(001)D.R.Bowler (2000). The unit cells were always ten layers deep, with the bottom two layers constrained to lie in bulk-like positions, and terminated in hydrogen. For simple structure calculations, each cell was one dimer row wide (using the p($`2\times 2`$) reconstruction) and sixteen dimers long (i.e. 384 atoms). When calculating kinking energies we used the same size cell, but four dimer rows wide (i.e. 1536 atoms). These methods have been used before with great success for calculations of many thousands of atoms (e.g. for calculation of step kinking energies on Si(001)D.R.Bowler and M.G.Bowler (1998)).

For an accurate assessment of promising candidates, we performed first principles density functional (DFT) simulations within the local density approximation (LDA) using the VASP codeG.Kresse and J.Furthmüller (1996). We used ultrasoft pseudopotentials, a plane wave cutoff of 150 eV (which gave energy difference and force convergence) and a Monkhorst-PackH.J.Monkhorst and J.D.Pack (1976) grid of $`4\times 4\times 1`$ k-points. Each unit cell was the same depth and width as the tight binding cells, but only eight dimers long, again with the bottom two layers held fixed in bulk-like positions and terminated in hydrogen. The change in energy caused by the shorter cell was found to be less than 0.1 eV; the longer tight binding cells were used to investigate nanoline-dimer repulsion, which could not be done using LDA.

A key question to consider when modelling the stability of structures with different chemical elements is how to compare structures with different stoichiometries. As the experiments we are considering were performed at high temperature on a Si(001) substrate, we can use the energy of bulk silicon to correct for different numbers of silicon atoms (though this raises questions about k-point sampling and error cancellation for the LDA calculations that will be considered below). For bismuth, however, we are faced with a harder problem. In the next section, we detail the approach used for calculating the periodicity of the high coverage Bi-($`2\times `$n) phase, and extract a sensible reference energy for Bi dimers on Si(001).

## III Energetics of Bi on Si(001)

At high coverages, Bi forms rows of dimers on Si(001), relieving the strain due to the large difference in lattice constants by having missing dimer defects at regular intervals, forming a high coverage, Bi-(2$`\times `$n) reconstruction. As we see Bi nanolines forming within this (2$`\times `$n) surface (and we also see that the nanolines are more stable than it – they form while it is desorbing) we start by describing it. We approach this problem, which will involve different amounts of Bi, and hence different stoichiometries, in the same way that we have approached modelling of the Ge/Si(001) surfaceJ.Oviedo et al. (2002). It is perfectly easy to calculate the energy of formation of trenches of missing dimers simply by removing Bi dimers, but to make contact with experimental results we need to find an energy for these Bi dimers that we have removed. We shall assume that the Bi-covered surface is in thermal equilibrium (which is easy to arrange experimentally), so that for a given number of Bi atoms we can write the probability of any given arrangement of these atoms as proportional to $`\text{exp}\left(E_a/kT\right)`$, where $`E_a`$ is the energy of the arrangement. The assumption of thermal equilibrium greatly simplifies the system that we are considering, while being perfectly legitimate at a number of different temperatures.

Now, as we are in thermal equilibrium, we only need know the relative energetics of different arrangements of the given number of Bi atoms. It will be useful to define a reference arrangement which we take to be a non-defective, lattice-matched layer of Bi covering a certain area of the Si(001) surface (with a boundary which we assume to be rectangular, but whose shape is arbitrary). Now we consider arrangements where we form regularly spaced trenches of missing dimers; conceptually we start by adding a number of dimers to the reference system equal to the number that we will remove and then forming the missing dimer trenches and relaxing. This conceptual approach allows us to compare different spacings by using the energy of Bi dimers in the perfect monolayer to correct for different amounts of Bi, which we find by calculating the difference in energy between a clean Si(001) surface and the same surface covered by a monolayer of Bi. Using our tight binding parameterisation we find this number to be -11.41 eV per Bi dimer.

We show the difference in energy per Bi dimer between a full, strained monolayer and the reconstructed surface for various values of period, n, in Fig. 1. There is a well defined minimum for periods of four and five (i.e. three and four Bi dimers) which is in excellent agreement with experimental observationsK.Miki et al. (1999a); Park et al. (1993). The modelling of this reconstruction gives us further confidence in our parameterisation, and presents a limit for the energetics of the Bi nanoline (which must be more stable than the Bi-($`2\times `$n) surface). Our tight binding simulations find that the energy per Bi dimer (adsorption energy and surface energy) in this reconstruction is -11.6 eV at n$`=4`$ and n$`=5`$.

When considering structures with different amounts of Bi, we can use the fact that the nanolines are observed to form in the high coverage Bi-(2$`\times `$n) phase to extract a sensible reference energy for Bi dimers. We have chosen the energy to be -11.4 eV per dimer, the energy of the reference phase, and used this where necessary. As mentioned above, we use the bulk Si energy to correct for differing numbers of Si atoms.

## IV Details of Successful Models

In this section, we present details of various of the models considered, and discuss how well tight binding has performed as a tool to discriminate between these models in the light of subsequent DFT modelling. First, we start by summarising the requirements placed on the models by experimental observations.

### IV.1 Experimental Constraints

We have been able to distil a number of requirements for the nanoline structure from our experimentsK.Miki et al. (1999a, b); J.H.G.Owen et al. (2002a, b):

* Stability: The lines form at a temperature where the high coverage (2$`\times `$n) Bi islands are evaporating, which implies that the energy per Bi dimer must be larger than these islands (in tight binding, more than -11.6 eV/dimer as given above). In RHEED experiments, the difference in desorption barrier for Bi from the islands and from the nanolines was determined to be 0.25 eVK.Miki et al. (1999a).
* Width, Registration and Appearance: The nanolines run perpendicular to the Si(001) surface dimer rows, occupying the space of four Si dimers (15.36 Å). They consist of two dimer-like features separated by $``$6.3 Åwhich lie between the Si(001) dimers, but are in line with the dimer rows.
* Straightness: Over many hundred observations of many lines (for instance in previous papersK.Miki et al. (1999a, b)), a kink has never been seen. Given the high temperature at which these nanolines form, it seems likely that a large thermodynamic kink energy is required, though a kinetic barrier cannot be ruled out.
* Defect Exclusion Zone: The nanolines are observedJ.H.G.Owen et al. (2002b) to repel defects and down step edges to a distance of 3-4 nm. As the defects and steps have a tensile stress field, we expect that the nanoline model should also have a tensile stress field.

### IV.2 Energetics from tight binding and DFT

Over the course of approximately a year, we have considered many different models. In general, we have first calculated the energy and kinking energy of possible models using tight binding. Only if these looked promising were they modelled with the more expensive DFT.

Within this section we present details of five unsuccessful models and what we believe to be the correct modelJ.H.G.Owen et al. (2002b), giving energies found with both tight binding and DFT. These results show that the tight binding calculations certainly ordered the structures correctly in terms of energy compared to DFT calculations (to within 0.1 eV) even though the absolute energies were sometimes rather too large, and thus confirm the use of this method as extremely valuable when searching a large configuration space. It is also worth noting that the tight binding simulations were computationally cheap (running in a matter of hours on a desktop workstation) and also that they deepened the collaboration between experiment and theory (many tight binding simulations were run by an experimentalist (JHGO), enabling detailed discussions of the remaining promising candidates). As an illustration, we found that tight binding relaxation calculations on a large unit cell used for kinking calculation (1536 atoms) took 2 hours and 4 minutes on a PowerMac G4/867MHz (using open-source g77 and ATLAS BLAS routines optimised for AltiVec) and 1 hour and 45 minutes on a Pentium 3/800MHz (using the Intel Fortran compiler for Linux and either ATLAS BLAS routines or the Intel Math Kernel library).

In Fig. 2 we show the ball-and-stick models of five unsuccessful candidates, and in Fig. 3 we show the final candidateJ.H.G.Owen et al. (2002b) (shown in off-axis views as the structures are rather complex). The five unsuccessful models shown are: the original structure proposedK.Miki et al. (1999a, b); D.R.Bowler (2000), which fitted original STM data but is only 3 Si dimers wide and has a very low kinking energy (labelled “Original”); the model proposed by Naitoh et al.M.Naitoh et al. (2000), which has the wrong spacing between features, is energetically poor and has a very low kinking energy (labelled “Two1DVs”); a structure made from back-to-back rebonded B-type steps, with Bi substituted in the second layer, which has the wrong appearance and is not energetically good enough (labelled “Btype4DV”); two Bi ad-dimers, which have no kinking energy and are energetically poor, and have no incentive to remain in this configuration (labelled “TwoAdDimers”); and two squares of Bi on top of a reconstruction similar to the final structureJ.H.G.Owen et al. (2002b) and to the structure of steps on As-covered Ge(001)S.B.Zhang et al. (2001).

We give the energies for these models, along with a variation on the Btype4DV structure (the same structure but with only a 3 dimer width, which is energetically rather better, but has the wrong spacing and a low kinking energy, labelled “Btype3DV” and not illustrated) and the final candidate (labelled “DoubleCore”) in Table 1. Where the number of silicon atoms has differed from the clean surface, we have corrected using an appropriate value for the bulk energy for silicon. In the DFT calculations, it is questionable to use simply an energy calculated from a bulk Si cell, as the k-point sampling will differ unless fully converged in both cells; we also expect error cancellation to be better in the same sized unit cells. We approached the problem from two directions: first, we calculated the total energy of a clean Si(001) surface using both an eight layer and a ten layer slab, and used the energy difference to give us a value for the energy of bulk Si. We compared this to calculations of a bulk Si cell which were fully converged with respect to k-point sampling (we used a $`4\times 4\times 4`$ mesh, and checked with an $`8\times 8\times 8`$ mesh), and found that the energy difference less than 0.05 eV per atom. The energy that we used was 5.96 eV/atom.

From these energies, we can see that the tight binding calculations have correctly ordered almost all of the structures, though the energies are sometimes too large. It is extremely pleasing that it has performed so well for such a large variety of structures, especially the final candidate which contains a large stress and significant reconstruction. In particular, it can be seen that there are a wide range of bond lengths (Bi-Bi are generally about 2.9 Å, but Bi-Si vary from 2.60 to 2.82 Å) and angles (from $`84^{}`$ to $`106^{}`$) which the tight binding calculations have modelled correctly. Bi seems to favour $`p^3`$ bonding, with 90 bond angles and a pair of electrons in an s orbital – and the tight binding calculations have modelled this well. The only point where the ordering is wrong is for the TwoSquares structure, which tight binding found to be more stable than it should have been. The Bi-Bi bonds are generally shorter than found in DFT calculations ($``$3.0-3.1 Å), suggesting that the Bi-Bi bond in TB might be too strong, accounting for the increased stability of this structure. The magnitudes of the TB energies are too large, probably because the model was fitted to bulk Bi-Bi and Bi-Si bonds, and this environment is at the surface, allowing Bi to form bonds with angles close to $`90^{}`$ (its preferred bonding angle).

The unsuccessful models were rejected for a number of reasons: they were energetically poor (Two1DVs, TwoAdDimers, TwoSquares); they had the wrong appearance (Btype4DV and Btype3DV); they had a low kinking energy (Original, Two1DVs, TwoAdDimers, Btype3DV); they had the wrong spacing of Bi dimers (Two1DVs, Btype3DV, Btype4DV); they had the wrong registry of Bi dimers relative to the substrate (Original, Two1DVs); and they were less stable than the high coverage Bi-(2$`\times `$n) surface (Two1DVs, TwoAdDimers, Btype3DV, Btype4DV, TwoSquares).

The final candidate (DoubleCore) has the correct appearance, a high kinking energy (3.75 eV from tight binding calculations), the correct spacing of features (6.3 Å), a high repulsion for defects and down step edges (out to 6 or 8 dimers or $``$3 nm), and is more stable than the high coverage Bi-($`2\times `$n) surface. We have discussed this structure in more detail elsewhereJ.H.G.Owen et al. (2002b).

## V Discussion and Conclusions

The key drawback with tight binding is that a parameterisation must be created for every interaction, and this is frequently a time-consuming procedure. Another extremely important problem for tight binding is that of transferrability: how well a parameterisation describes interactions that are dissimilar to the environment used for fitting. The Bi-Bi and Bi-Si parameterisation used was created simply and quicklyD.R.Bowler (2000), using the universal tight binding parameters of Harrison for Bi-Bi hoppingW.A.Harrison (1980) and fitting the scaling to the bulk modulus only; the parameterisations have only been used in environments close to those used for fitting. This simple procedure (which has also been applied to Ge-Ge and Ge-SiD.R.Bowler (2002)) has been shown here to be remarkably effective.

However, there is clearly a need for more reliable ways to explore parameter space in situations such as that encountered in this paper, and developments in the field known as ab initio tight binding are providing the answers. Starting with the Harris-Foulkes theoryJ.Harris (1985); W.Foulkes and R.Haydock (1989), the link between DFT methods using local, atomic-like orbitals and tight binding was shown. This was quantified in an elegant demonstration that tight binding could be derived from DFT via a well-defined series of approximationsA.P.Sutton et al. (1988). Ideas such as these have led to different approaches to improving and generalising tight binding.

Typically, confined atomic orbitals (found by solving the Schrödinger equation within a confining potential)O.F.Sankey and D.J.Niklewski (1989) are used as a basis set, adding orbitals beyond the minimal set as necessary for completeness. The problem with any approach like this is to decide where to make approximations. There are two key strands to these approximations: to neglect self consistency (but to retain three-centre integrals)O.F.Sankey and D.J.Niklewski (1989); and to add self-consistency but to neglect three-centre termsD.Porezag et al. (1995). Of course, there have been developments such as adding self-consistency to the three-centre approachA.A.Demkov et al. (1995).

The ultimate generalisation of these ideas lies in local orbital-based DFT methods such as SiestaD.Sánchez-Portal et al. (1997), PlatoA.P.Horsfield (1997); S.D.Kenny et al. (2000) and ConquestD.R.Bowler et al. (2000b), which point the way forward, as they enable the creation of the matrix elements (equivalent to a parameterisation) accurately, both in advance and in real time: the link between tight binding methods and DFT-based methods is becoming ever stronger. Clearly the level of approximation chosen needs to trade off against the size and complexity of the system being modelled.

In conclusion, we have presented tight binding and DFT calculations for the high coverage Bi-($`2\times `$n) surface which are in good agreement with experiment, and we have explored a variety of possible structures for Bi nanolines. Tight binding calculations have been shown to agree well with DFT calculations in terms of ordering structures, and enabled us to perform quick, efficient searches of parameter space. In particular, it has enabled us to find an extraordinary structure for the Bi nanoline, which in turn sheds light on the role of surface stress in self-assembly of nanostructures on semiconductor surfaces.

###### Acknowledgements.

DRB thanks the UK Engineering and Physical Sciences Research Council and the Royal Society for support through fellowships. JHGO is supported by the Japanese Science and Technology Agency (JST) as an STA Fellow. This study was performed through Special Coordination Funds of the Ministry of Education, Culture, Sports, Science and Technology of the Japanese Government (Research Project on active atom-wire interconnects). We are happy to acknowledge useful discussions with Dr. K.Miki, Prof. G.A.D.Briggs, Dr. W.McMahon and Prof. M.J.Gillan. DFT calculations were carried out in the HiPerSPACE Centre at UCL.


format: text

text:
# EVIDENCE FOR AN OUTER DISK IN THE PROTOTYPE ‘COMPACT ELLIPTICAL’ GALAXY M32

## 1 Introduction

The 32nd object in the catalog of Messier (1850) has come to be known as the archetype of high surface brightness, low-luminosity, compact elliptical galaxies (de Vaucouleurs 1961). It has been proposed that they may be the dense cores of tidally truncated, or at least modified, ordinary elliptical galaxies (King 1962; Faber 1973; Nieto & Prugniel 1987; Choi, Guhathakurta, & Johnston 2002, and references therein). It has also been suggested that M32 (NGC 221) may in fact never have been an elliptical galaxy, but is instead the bulge of a (partially) stripped disk galaxy (Bekki et al. 2001; see also Nieto 1990).

Using the tight ($`r_s=0.91`$) correlation between the central concentration index $`C_{r_\mathrm{e}}(1/3)`$ of a bulge111By the term bulge it is meant both an elliptical galaxy and the bulge of a disk galaxy. and the mass of its central supermassive black hole (SMBH; Graham et al. 2001a), provides a new constraint capable of determining which proposition, if either, is correct. Simply by modelling M32’s surface brightness profile as either a one-component ($`r^{1/n}`$) elliptical, or as a 2-component ($`r^{1/n}`$ bulge + exponential) lenticular or spiral galaxy, and knowing its SMBH mass, allows one to decipher which scenario is more probable.

The paper is laid out as follows. Section 2 introduces M32’s surface brightness profile which is modelled in Section 2.1. Section 3 introduces M32’s velocity dispersion measurements and SMBH mass estimates into the discussion, supporting the notion that the bulge and disk decomposition is the correct one.

## 2 The surface brightness profile of M32

M32’s major-axis, $`R`$-band surface brightness profile presented in Kent (1987) is reanalyzed here. The data was obtained using an RCA CCD attached, at different times, to three telescopes with different fields of view at the Whipple Observatory on Mount Hopkins. The observing procedures are described in Kent (1983) and the (remarkably standard) reduction procedure given in Kent (1987).

The main complication with this galaxy is its proximity to M31 (Andromeda). It resides (in projection) within the outer disk of M31, and so the disk of M31 had to first be modelled and then subtracted. Kent wrote, “The light from M31 was removed approximately by fitting second order polynomials to the background light in the frame of M32 obtained with the Bausch & Lomb 8000 \[a 20 cm aperture telescope with a $`15\mathrm{}\times 25\mathrm{}`$ field of view\] excluding an area about M32 itself”. Ultimately, uncertainties in the sky-background level resulted in the termination of the profile at $`\mu (r=280`$$`)`$24 $`R`$-mag arcsec-2.

Further complications and solutions are described in Kent (1987), and won’t be repeated here for the following reason. After commencing this work, Choi et al. (2002) presented new $`B`$\- and $`I`$-band surface brightness profiles for M32. These show very good agreement over the radial range in common with Kent (1987), and extend to 420$`\mathrm{}`$. This lends confidence that the publicly available data from Kent (1987) is reliable. The profiles of Choi et al. (2002) will however be referred to in a qualitative manner.

### 2.1 Modelling the surface brightness profile

Inner components such as nuclear disks, star clusters, flattened cores, etc., are known to reside within the central $`1\mathrm{}`$ of many galaxy’s (e.g., Rest et al. 2001; Ravindranath et al. 2001). We wish to avoid such features here as we are presently only concerned with the bulge (and outer disk, if one exists) of M32. However, typically observed galaxies are considerably more distant than M32. At a distance of 0.8 Mpc, with the exception of the Milky Way, M32 is some 10 times closer than any other galaxy with a positive SMBH detection (Kormendy & Gebhardt 2001; Merritt & Ferrarese 2001b). Its nuclear region is therefore very well resolved: 1$`\mathrm{}`$ is equivalent to 3.87 pc. If M32 was located at the distance of the Virgo cluster, exclusion of its central arcsecond would translate to the removal of the inner $``$ 60-80 pc ($``$ 15-20$`\mathrm{}`$ at its present location).

Böker et al. (2001) chose to model the central excesses (above that of the bulge) found in HST images of spiral galaxies with a power-law. Tonry (1984) noted that the inner 10$`\mathrm{}`$ of M32 displays a power-law shape. This may therefore be a somewhat common feature of bulges222In passing, it is noted that this may be a feature that results in the over-estimation of the Sérsic index $`n`$ when modelling ground-based images which have not resolved, or avoided, such central excesses.. Michard & Nieto (1991) presented evidence that the inner $`5\mathrm{}`$ of M32 contains a nuclear disk, however this has subsequently been refuted by Lauer et al. (1998) who found no photometric evidence for separate nuclear components in HST images. However, the presence of an apparent excess central flux did result in Kent (1987) excluding the inner 15$`\mathrm{}`$ from his model fitting, and Choi et al (2002) excluding the inner 10$`\mathrm{}`$, which is also done here. Apart from this, the entire surface brightness profile beyond 10$`\mathrm{}`$ is modelled here. The data from Table 3 of Kent (1987), and also a resampling of Kent’s Figure 1 (to give an equal spacing in radius, effectively providing a more even weighting \[radially\] to the data), are modelled.

Both a seeing-convolved $`r^{1/n}`$ model and a seeing-convolved ($`r^{1/n}`$ \+ exponential) model were fitted. In both cases, all model parameters were simultaneously fitted using the quasi-Newton, non-linear least-squares algorithm UNCMND (Kahaner, Moler, & Nash 1988) which was iterated until convergence on the optimal solution giving the smallest $`\chi ^2`$ value. All parameters were allowed to range freely in the fitting process, the only constraint was that they must be positive real numbers. The seeing was reported by Kent (1987) to have a FWHM of 1.3$`\mathrm{}`$, and therefore, due to the exclusion of the inner 10$`\mathrm{}`$, has little effect on the results which are shown in Figure 1. Several truncations of the outer profile were explored, but the overall conclusion was always the same: M32 cannot be modelled as, that is to say it is not (structurally), a single component system. The curvature in the residual profile of Figure 1a is classic evidence of this. Importantly, M32 can be described exceptionally well with the standard two-component model used these days to model disk galaxies with bulges333Fitting an ($`r^{1/4}`$ \+ exponential) model resulted in 86% more scatter than using an ($`r^{1/n}`$ \+ exponential) model, and a rather unimpressive fit.. The notably small residuals are compelling evidence that this model is likely to be correct. This is, however, not to say that M32 has a rotationally supported disk of stars, only that is has an outer exponential distribution of stars. To help evaluate if the two-component model is indeed correct, let us look at the resulting structural parameters to see if they are consistent with those of known bulge/disk systems.

The bulge-to-disk size ratio ($`r_\mathrm{e}/h`$) is 0.20, in good agreement with that of normal disk galaxies (Graham 2001, and references therein) and suggests nothing unusual. The logarithm of the bulge-to-disk luminosity ratio ($`\mathrm{log}B/D`$) is 0.22 and typical of an S0 galaxy. The effective half-light radius of the bulge is 27$`\mathrm{}`$, in reasonable agreement with the value of 32$`\mathrm{}`$ derived by Kent when modelling the radial interval $`15\mathrm{}<r<100\mathrm{}`$, and even agrees well with the value of 30$`\mathrm{}`$ found by de Vaucouleurs (1953). The effective bulge surface brightness is 18.23 $`R`$-mag arcsec-2, the central bulge surface brightness is 15.31 $`R`$-mag arcsec-2, and the absolute magnitude of the bulge is 16.34 $`R`$-mag.

Further support for the above bulge/disk decomposition comes from the new location of M32 in several structural parameter diagrams for bulges. M32 is a distant outlier in the insightful $`\mathrm{log}n`$-$`B_T`$ diagram presented in Jerjen & Binggeli (1997; their figure 2). Noting that their value of $`n`$ corresponds to our value of $`1/n`$, when $`n=1.5`$ (in our notation) M32 moves up into the very center of points, and indeed the very center of the relation defined by the dwarf galaxies (see also Graham 2001, his figure 14). The reason for this is that the disk had biased Jerjen & Binggeli’s one-component $`r^{1/n}`$ fit to M32’s light profile, resulting in a value of $`n`$ which is $``$5 in our notation. This additionally explains the deviant nature of M32 in the $`\mathrm{log}r_0`$-$`B_T`$ diagram, and most of the discrepancy in the $`\mu _{0,\mathrm{bulge}}`$-$`B_T`$ diagram. The suggestion by Wirth & Gallagher (1984) that compact elliptical galaxies may be the extension of brighter elliptical galaxies, based on color, luminosity, and velocity dispersion (and the existence of isolated compact elliptical galaxy candidates) would appear to be correct in the sense that the bulge of M32 appears to be the faint extension of bulges and elliptical galaxies in general.

What of the disk? Kent (1987) remarked that there “seems to be an excess of light at large radii, with the excess having an exponential profile”. The ellipticity profile of Kent (1987) also suggests the presence of a distinct outer component in M32, rising from $`ϵ`$=0.11 at 150$`\mathrm{}`$ to $`ϵ`$=0.19 by $`200\mathrm{}`$ (where the ellipticity was then held constant). This feature is even more clearly evident in the ellipticity profile of Choi et al. (2002), rising steadily from $`ϵ0.14`$ at 100$`\mathrm{}`$ to $`0.35`$ at 250$`\mathrm{}`$. The position angle also changes notably at 100$`\mathrm{}`$-150$`\mathrm{}`$.

The central disk surface brightness is 21.28 $`R`$-mag arcsec-2 (Figure 1b). Applying the standard disk inclination correction $`2.5C\mathrm{log}(1ϵ)`$, where $`ϵ`$ is the ellipticity of the disk ($`ϵ`$0.3 from Choi et al. 2002) and $`C`$=0.5 in the $`R`$-band (Tully & Verheijen 1997), would give a face-on central disk surface brightness of 21.48 $`R`$-mag arcsec-2. The definition of a low surface brightness (LSB) galaxy, is one in which the central disk surface brightness is more than one magnitude fainter than the canonical Freeman (1970) value of 21.65 $`B`$-mag arcsec-2 (which is about 20.5 $`R`$-mag arcsec-2). Thus, M32 would just about qualify as a bulge-dominated LSB disk galaxy (Beijersbergen, de Blok, & van der Hulst 1999; Impey & Bothun 1997). It is, however, because of it’s size and magnitude, a dwarf galaxy.

The scale-length of the disk is 130$`\mathrm{}`$ or 0.5 kpc. There is, however, possible evidence of a disturbance in the very outer profile; it turns downward from an exponential profile at around 250$`\mathrm{}`$. That is, there appears to be a lack of light, relative to the exponential part of the profile, beyond $`250\mathrm{}`$. This behavior is visible in both the data of Kent (1987) and also that of Choi et al. (2002). It is accompanied by a marked change in the behavior of the ellipticity profile of Choi et al. (2002), flattening (or even decreasing slightly) beyond 250$`\mathrm{}`$. This may be a sign that material has been stripped away from the outer disk, although it is stressed that this conclusion is largely speculative. However, if true, this turnover would be biasing the surface brightness profile fit, making the disk scale-length appear shorter than it was before tidal-stripping commenced. Truncating the (equally spaced) profile at 220$`\mathrm{}`$, to avoid the potentially stripped outer disk, and thereby (possibly) sampling only the original undisturbed profile gives $`h=336\mathrm{}`$ (1.27 kpc) and $`\mu _0=22.42`$ for the disk. For the bulge, $`r_\mathrm{e}=29\mathrm{}`$ and $`n=1.99`$ (see Figure 2; using the profile which is logarithmically spaced in radius gives $`n=2.08`$). Here the value of $`\mathrm{log}B/D`$=-0.09, and $`r_\mathrm{e}/h`$=0.09, which are again not unreasonable. An $`r^{1/n}`$-only model fails to provide a convincing fit to this truncated radial range.

## 3 Discussion

The surface brightness profile of M32 can be modelled remarkably well as a combination of an $`r^{1/n}`$ bulge and an exponential disk of stars. The ellipticity profile additionally supports this interpretation of the data. However, could, instead, M32’s surface brightness profile be so disturbed that the outer exponential envelope is actually due to material pulled off from what was once a one-component $`r^{1/4}`$ (i.e. $`n=4`$) elliptical galaxy? This scenario appears unlikely for the following reasons.

Firstly, somewhat persuasive evidence comes from the the central velocity dispersion of M32 which has been measured to be 76$`\pm `$10 km s-1 (van der Marel et al. 1998). This figure is in good agreement with the average value of 74 km s-1 obtained from numerous estimates listed in Hypercat.444Hypercat can be reached at http://www-obs.univ-lyon1.fr/hypercat/ It also agrees with Gebhardt et al.’s (2000) estimate of 75 km s-1 for the luminosity weighted velocity dispersion within one effective radius. Recently, it has been discovered that the central velocity dispersion of a bulge correlates strongly ($`r=0.8`$) with the shape of the bulge light profile (as measured with the Sérsic index: Graham et al. 2001b; Graham 2002). Galaxies with measured velocity dispersions less than 100 km s-1 are observed to have values of $`1<n<2`$. The measured value of $`n=1.5`$ for the bulge of M32 is thus exactly what one would expect from its dynamics, not $`n=4`$. Bulges (which includes elliptical galaxies) with values of $`n4`$ have velocity dispersions typically greater than 100 km s-1.

The second, related, line of reasoning comes from M32’s SMBH mass. From ground based kinematical data, Tonry (1984,7) predicted a SMBH mass of (3-10)$`\times 10^6M_{\mathrm{}}`$ at the center of M32. This pioneering work has been confirmed with HST data, from which van der Marel et al. (1998) derived a SMBH mass of (3.4$`\pm `$0.7)$`\times 10^6M_{\mathrm{}}`$, and more recently Joseph et al. (2001) find a mass of (2-4)$`\times 10^6M_{\mathrm{}}`$. The velocity dispersion ($`\sigma `$) and SMBH mass ($`M_{\mathrm{bh}}`$) combination of M32 is known to fall on the $`\mathrm{log}M_{\mathrm{bh}}`$-$`\mathrm{log}\sigma `$ relation for (non-disturbed) ellipticals and bulges (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Merritt & Ferrarese 2001a). Thus, unless the process of tidal stripping modifies both the central velocity dispersion and the SMBH mass in such a way that it preserves the $`\mathrm{log}M_{\mathrm{bh}}`$-$`\mathrm{log}\sigma `$ relation, the central structure and dynamics of M32 are likely to be that of the original (i.e. undisturbed) galaxy. Indeed, due to the higher densities at smaller radii, the central velocity dispersion and hence inner mass distribution and therefore inner light profile are expected to remain largely unaffected by the outer stripping process. Taken with the result in the previous paragraph, this strongly suggests that the $`n=1.52.0`$ bulge profile is the original shape of the bulge.

Figure 3 shows the $`\mathrm{log}M_{\mathrm{bh}}`$-$`\mathrm{log}n`$ relation for bulges. It is a variant of the relation between SMBH mass and central bulge concentration ($`C_{r_\mathrm{e}}[1/3]`$) shown in Graham et al. (2001a). (The parameter $`C_{r_\mathrm{e}}(1/3)`$ is a monotonicly increasing function of $`n`$, Trujillo et al. 2001.) If M32 did ever have an $`r^{1/4}`$ profile, then (from Figure 3) its SMBH mass should have once been and should still be $`10^8M_{\mathrm{}}`$, and certainly greater than $`10^7M_{\mathrm{}}`$. Given the bulge value of $`n=1.5`$ agrees with the actual SMBH mass estimate (and with M32’s velocity dispersion), the stars composing the outer exponential envelope are almost certainly an excess (relative to the bulge) which have not come from a reshaped bulge555The reason Choi et al. (2002) claimed to be able to fit an $`r^{1/4}`$ bulge to the inner profile of M32 is likely because of the limited radial range $`10\mathrm{}<r<30\mathrm{}`$ they used to do this..

All of this is not to say that the outer profile of M32 has not been disrupted, Indeed, the deficit of stars in the outer profile, causing the downward kink (or break) in the disk, may have been due to gravitational stripping by M31. This deficit is also visible in the $`B`$\- and $`I`$-band profiles of Choi et al. (2002) and hence less likely to be a systematic error in the profile extraction technique of Kent (1987). Indeed, the recently reported tidal stream of metal-rich stars around M31 is thought to have likely come from M32 and/or NGC 205 (Ibata et al. 2001). The presence of a disk may also explain the intermediate-age (5-15 Gyr) stellar population in M32 (e.g. Grillmair et al. 1996; Davidge 2000).

Perhaps a comment on nomenclature would be beneficial here. Kormendy & Gebhardt (2001) used the words ‘compact’ and ‘fluffy’ when referring to bulges. It should be noted that this has no reference at all to varying profile ‘shape’ (or equivalently ‘concentration’, according to the mathematical definition given in Trujillo et al. 2001, their equation 5 and 6). That is, no reference to departures from the $`r^{1/4}`$ law, or structural homology, are implied by their terms ‘compact’ and ‘fluffy’. The $`r^{1/4}`$ law only has a horizontal scale term ($`r_\mathrm{e}`$) and a vertical scale term ($`\mu _\mathrm{e}`$ or $`\mu _0`$), the shape, or concentration, is exactly the same for all $`r^{1/4}`$ profiles. Kormendy & Gebhardt (2001) discussed variations in $`\mu _e`$ and $`r_e`$ when they referred to ‘compact’ and ‘fluffy’. Hence, although M32 is regarded as compact, it’s central concentration $`C_{r_\mathrm{e}}(1/3)`$ is actually rather low.

Bekki et al.’s (2001) N-body/SPH simulations of tidal interactions between M31 and an orbiting early-type spiral galaxy predict either a complete stripping of the disk, or at least a vertical heating of the satellite’s disk to create a thick disk. Clearly a faint disk, with very little gas (Welch & Sage 2001), still surrounds M32. Much of the gas and stars may indeed have been stripped away, resulting in the low surface brightness disk. Also possible, is the suggestion by Bekki et al. (2001) that tidal interactions with M31 funnelled some of M32’s gas to its center, forming a massive starburst (see also Noguchi & Ishibashi 1986). This could account for the excess central flux within the inner $`10\mathrm{}`$ of M32 having an age of $``$4 Gyrs (Vazdekis & Arimoto 1999; del Burgo et al. 2001).

Compact elliptical galaxies are a rare class of objects. A closer inspection of such objects seems warranted in order to inspect whether the species is indeed real, or simply a case of misclassification.

I wish to thank Peter Erwin for providing me with Kent’s (1987) surface brightness profile of M32, resampled with equal spacing in radius, and for useful discussions which helped to shape this paper. I am also grateful to Carme Gallart for her comments on this work.


text:
# Mesons in a Poincare Covariant Bethe–Salpeter Approach

## I Introduction

Recently, the scalar mesons have attracted a lot of interest as the reanalysis of the pseudoscalar meson scattering data indicated the existence of a flavor SU(3) nonet in this channel . It is therefore desirable to gain deeper understanding of the constituent structure of the scalar mesons together with a comprehensive description of the meson states in the other spin–parity channels. The ultimate goal would be to understand all low–lying meson states and resonances as non–perturbative bound states in Quantum Chromo Dynamics (QCD).

A relativistic framework for analyzing mesons as composite objects is provided by the Bethe–Salpeter equations that extract poles in the quark–antiquark scattering kernel. The attraction needed to bind quarks and antiquarks emerges from dressed multiple gluon exchange. Thus the essential ingredients to these equations are the quark and gluon propagators as well as the quark–gluon vertex. In addition, these $`n`$-point Green’s functions are related by their Schwinger–Dyson equations which are part of an infinite tower of non–linear integral equations. There has been some progress in the understanding of the infrared behavior of the gluon propagator from recent Yang–Mills lattice calculations  as well as from studies of the coupled system of gluon and ghost Schwinger–Dyson equations . Nevertheless, for phenomenological applications the frequently adopted strategy is to model the gluon propagator as well as the quark–gluon vertex and consistently derive the quark propagator from its Schwinger–Dyson equation because it facilitates the continuation to complex Euclidean momenta.

These types of calculations have a long history, for reviews see Refs. . Early versions adopted pointlike gluon propagators in coordinate space that eventually lead to Nambu–Jona–Lasinio (NJL) type models , pointlike propagators in momentum space were also considered . These models are particularly simple because either solving the Schwinger–Dyson equation yields a free quark propagator or the Bethe–Salpeter integral equations reduce to algebraic equations. The main target particularly of the NJL–model studies have been the pseudoscalar mesons. It turned out that they can be adequately described once the important feature of dynamical chiral symmetry breaking is incorporated, i.e. the interaction is strong enough so that the resulting quark propagator signals a non–vanishing quark condensate reflected by a non–zero constituent quark mass. Then the pseudoscalar mesons can be understood as the would–be Goldstone bosons of chiral symmetry breaking. However, as the Bethe–Salpeter equations involve the dressed (constituent) quark propagators, binding can only be achieved kinematically and meson states with masses larger than twice the constituent quark mass cannot be described consistently. For that reason, model gluon propagators have been developed that yield quark propagators without poles for real momenta as an attempt to include the confinement phenomena . Again these studies focused on pseudoscalar mesons  while a comprehensive investigation for the scalar, pseudoscalar, vector and axialvector mesons has not been carried out so far. Other studies  made contact with perturbative QCD by considering a model gluon propagator that matches the pertinent anomalous dimension. This contribution has negligible effect on the meson properties, but its inclusion makes cumbersome the extraction of the solutions to the Schwinger–Dyson equations for the large time–like momenta that enter the Bethe–Salpeter equations with meson states other than the pseudoscalars. We consider this attempt interesting, but as we do not want to focus on the ultraviolet properties of mesons here, we regard it an unnecessary technical complication. Therefore we focus rather on establishing a model as simple as possible that we consider a pertinent starting point to study the structure and properties of low–lying mesons in a relativistic framework. Our model interaction is parameterized in form of a non–trivial gluon propagator that contains sufficient strength to cause dynamical chiral symmetry breaking. For technical reasons it turns out that a Gaussian shape function for the propagator in momentum space is most suitable. Essentially we consider this model propagator as an effective interaction that relativistically describes the binding of quarks and antiquarks to mesons. Furthermore, we take the quark–gluon vertex function to be the tree level one since this procedure provides a framework that is consistent with chiral symmetry when the ladder approximation for the Bethe–Salpeter equation is employed . For approaches going beyond ladder approximation see e.g. Ref. .

This paper is organized as follows. In Section II we will introduce the effective interaction and solve the Schwinger–Dyson equation for the quarks. We will put particular emphasis on the analytic continuation of the resulting quark propagator to time–like momenta that enter the Bethe–Salpeter equations. We will discuss the structure of the Bethe–Salpeter equations and present the solutions in Section III. In Section IV we will present the numerical results in the sector of the three light quarks (up, down and strange) and compare them to the empirical data. In Section V we will conclude and suggest a possible extension of the current approach in particular with regard to the possibility that the scalar meson might have to be considered as two–quark – two–antiquark bound states . We devote an appendix to discuss the numerical stability of our results.

## II The Quark Schwinger-Dyson Equation

The purpose of this section is twofold. First we will explain our model gluon propagator whose dressing parameterizes the effective interaction and then solve the corresponding Schwinger–Dyson equation for the quark propagator.

As already discussed in the introduction we take a Gaussian form for dressing the model gluon propagator. This follows the work of Ref.  with the exception that we omit the logarithmic tail that matches perturbative QCD because we are not interested in the ultraviolet properties of mesons. We therefore write,

$$g^2G_{\mu \nu }^{ab}(q)=4\pi ^2D\delta ^{ab}t_{\mu \nu }(q)\frac{q^2}{\omega ^2}\mathrm{exp}\left(\frac{q^2}{\omega ^2}\right)$$
(1)
where $`\mu ,\nu `$ are Lorentz indices, $`t_{\mu \nu }(q)`$ is the transverse momentum projector and $`a,b`$ label color. While the prefactors in eq. (1) are chosen to make subsequent equations more concise, $`D`$ and $`\omega `$ are dimensionful parameters that we will determine from fitting empirical data. The coefficient $`D`$ sets the strength of the interaction and $`\omega `$ is the value at which the scalar function in the parameterization is maximal. Hence $`\omega `$ sets the interaction scale. The dressed gluon propagator (1) is supposed to represent a sensible hadron model and hence one can envisage that $`\omega `$ will have a value of several hundred $`\mathrm{MeV}`$.

We interpret the effective interaction (1) as the propagator (in Landau gauge) of a gluon that gets absorbed and emitted by the quarks that eventually get bound to form mesons. To completely define the interaction, we need to parameterize the quark–gluon coupling. To establish chiral symmetry we apply the rainbow–ladder approximation to the system of Schwinger–Dyson and Bethe–Salpeter equations. This implies that the quark–gluon coupling is given by the tree level interaction vertex, $`ig\gamma _\mu \frac{\lambda ^a}{2}`$, where $`\lambda ^a`$ is a Gell–Mann matrix acting in color space. Note, that we have already included the coupling constant $`g`$ in the definition of the effective interaction (1).

Then the Schwinger–Dyson equation for the (inverse) quark propagator becomes ($`\overline{}d^4k=d^4k/(2\pi )^4`$)

$$S^1(p)=i\text{ / }p+m_0+\overline{}d^4k\gamma _\mu S(k)\gamma _\nu g^2\frac{\lambda ^a}{2}\frac{\lambda ^b}{2}G_{\mu \nu }^{ab}(kp)$$
(2)
where $`m_0`$ is the current mass of the considered quark. This contribution represents the only explicit distinction between quarks of different flavors. Of course, its effects will implicitly propagate through the whole calculation. However, for notational simplicity we will continue to suppress flavor labels. A suitable parameterization of the quark propagator is inspired by the form of a free fermion propagator

$$S(p)=\left[\frac{1}{i\text{ / }pA(p^2)+B(p^2)}\right].$$
(3)
In solving the Schwinger–Dyson equation (2) we have to find the scalar functions $`A(p^2)`$ and $`B(p^2)`$. It is also very instructive to define a mass function via $`M(p^2)=B(p^2)/A(p^2)`$. In particular $`M(p^2=0)`$ plays the role of a constituent quark mass and a large value thereof signals dynamical chiral symmetry breaking.

As usual we work in Euclidean space with Hermitian Dirac matrices***This can be related to the standard Minkowski space Dirac matrices via $`\gamma _4=\gamma _0^M,\gamma _j=i\gamma _j^M`$. that obey $`\{\gamma _\mu ,\gamma _\nu \}=2\delta _{\mu \nu }`$ and $`\gamma _5=\gamma _1\gamma _2\gamma _3\gamma _4`$. Inserting the effective interaction (1) and performing the standard trace algebra, we then deduce the following coupled equations for the propagator functions

$`A(x)`$ $`=`$ $`1+{\displaystyle \frac{D}{\omega ^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dyyA(y)}{(yA^2(y)+B^2(y))}}`$ (5)
$`\times {\displaystyle \frac{2}{\pi }}{\displaystyle _1^1}𝑑z\sqrt{1z^2}\left[{\displaystyle \frac{2}{3}}y+\left(1+{\displaystyle \frac{y}{x}}\right)\sqrt{xy}z{\displaystyle \frac{4}{3}}yz^2\right]\mathrm{exp}\left\{{\displaystyle \frac{x+y2\sqrt{xy}z}{\omega ^2}}\right\},`$
$`B(x)`$ $`=`$ $`m_0+{\displaystyle \frac{D}{\omega ^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dyyB(y)}{(yA^2(y)+B^2(y))}}`$ (7)
$`\times {\displaystyle \frac{2}{\pi }}{\displaystyle _1^1}𝑑z\sqrt{1z^2}\left[x+y2\sqrt{xy}z\right]\mathrm{exp}\left\{{\displaystyle \frac{x+y2\sqrt{xy}z}{\omega ^2}}\right\},`$
where the four dimensional integral measure has been expanded such that $`x=p^2,y=k^2`$ and $`z=pk/\sqrt{p^2k^2}`$.

In a first step we solve eqs. (5) and (7) for spacelike momenta, i.e. for real positive $`x`$. Then we observe that the integrals on the $`RHS`$ of these equations only involve the propagator functions at real arguments $`y`$ and we can use them to numerically compute the propagator functions for arbitrary complex $`x`$. At first sight, it appears that $`A(x)`$ and $`B(x)`$ could not be consistently continued because the cut along the negative $`x`$–axis (associated with $`\sqrt{x}`$) would yield different results when continuing in the upper or the lower half–plane and it would be impossible to resolve the ambiguity in $`\sqrt{x}\pm i\sqrt{\xi }`$ when continuing $`x\xi `$. Fortunately, this is not an obstacle. By expanding the exponential functions in eqs. (5) and (7) it becomes obvious that all terms that contain $`\sqrt{x}`$ are odd in $`z`$ and thus vanish when integrating over this angular variable. Thus we are free to choose either of the two signs above. For definiteness we work in the upper half–plane with $`\sqrt{x}i\sqrt{\xi }`$ along the negative half–line.

The observation that it is sufficient to know the propagator functions along the real positive axis to compute them for all complex arguments of interest is a major point of this present study. It is an especially important issue to the study of the Bethe–Salpeter equation, since the quark propagator must be sampled in the complex plane of Euclidean momenta. In other words, the quark functions are calculated at precisely the momenta for which they are used, with no fitting functions, interpolation or extrapolation. However, in order to do so numerically, the quark propagator functions must be known very accurately along the real positive axis. This is achieved by noting that in eqs. (5) and (7) the angular ($`z`$) integrals can be done analytically for real and positive $`x`$. This property of the angular integrals is a particular feature of the Gaussian dressing function (1). The resulting equations read

$`A(x)`$ $`=`$ $`1+D{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dyyA(y)}{(yA^2(y)+B^2(y))}}\mathrm{exp}\left\{{\displaystyle \frac{x+y}{\omega ^2}}\right\}`$ (9)
$`\times \left\{\left(1+{\displaystyle \frac{y}{x}}+2{\displaystyle \frac{\omega ^2}{x}}\right)I_2\left({\displaystyle \frac{2\sqrt{xy}}{\omega ^2}}\right)2{\displaystyle \frac{\sqrt{y}}{\sqrt{x}}}I_1\left({\displaystyle \frac{2\sqrt{xy}}{\omega ^2}}\right)\right\},`$
$`B(x)`$ $`=`$ $`m_0+D{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dyyB(y)}{(yA^2(y)+B^2(y))}}\mathrm{exp}\left\{{\displaystyle \frac{x+y}{\omega ^2}}\right\}`$ (11)
$`\times \left\{\left({\displaystyle \frac{\sqrt{y}}{\sqrt{x}}}+{\displaystyle \frac{\sqrt{x}}{\sqrt{y}}}\right)I_1\left({\displaystyle \frac{2\sqrt{xy}}{\omega ^2}}\right)2I_2\left({\displaystyle \frac{2\sqrt{xy}}{\omega ^2}}\right)\right\},`$
where $`I_n`$ are modified Bessel functions. These equations are one–dimensional coupled non–linear integral equations which can be straightforwardly solved numerically with a high degree of precision. The solutions for real (Euclidean) momenta are subsequently substituted into eqs. (5) and (7) to yield the propagator functions for complex momenta.

## III The Bethe-Salpeter Equation

Having obtained the quark propagators in the complex plane from the Schwinger–Dyson equations we have collected all ingredients for the Bethe–Salpeter integral equations. They will ultimately yield the quark meson vertex functions, $`\mathrm{\Gamma }`$, that describe mesons as bound quark–antiquark pairs. Strictly speaking, the Bethe–Salpeter equation is an eigenvalue problem, valid only at the resonance pole $`P^2=M^2`$, where $`M`$ is the mass of the resonance. It is derived from considering the four–point quark Green’s function that involves exchanges of resonance mesons. Such an exchange is characterized by a pole in that four–point function, and the homogeneous Bethe–Salpeter equation below (12) determines the position of that pole. All other regular terms in the vicinity of this pole are neglected. Demanding furthermore that the residue of this pole is unity yields the normalization condition (23) for the vertex functions.

### A Bethe–Salpeter vertex functions

The vertex functions resulting from the Bethe–Salpeter equation are characterized by three momenta out of which only two are linearly independent due to momentum conservation at the vertex. If we denote the meson momentum $`P`$ and the momentum of the incoming quark $`p+\xi P`$ then the momentum of the outgoing quark (= incoming antiquark) is $`p+(\xi 1)P`$. This suggests to label the vertex functions by $`p`$ and $`P`$: $`\mathrm{\Gamma }(p,P)`$. We have also introduced the arbitrary momentum partition parameter $`\xi [0,1]`$. Due to strict relativistic covariance the results for physical observables do not depend on $`\xi `$. Unfortunately we will have to assume approximations within the full numerical computation (see Section 4) that violate covariance to some extent. We will study the $`\xi `$–dependence of our results and verify that within a wide range $`\xi `$–independence is maintained, see Appendix A. This will represent an a posteriori validation for the relativistic covariance of our computations.

We now turn to the main target of our studies, the Bethe–Salpeter integral equations for the vertex function $`\mathrm{\Gamma }(p,P)`$ in ladder approximation :

$$\mathrm{\Gamma }(p;P)=\frac{4}{3}\overline{}d^4k\left[\gamma _\nu S(k+\xi P)\mathrm{\Gamma }(k;P)S(k+(\xi 1)P)\gamma _\mu \right]g^2G_{\mu \nu }(kp).$$
(12)
Here we have factorized the color factors in the effective interaction, $`G_{\mu \nu }^{ab}(q)=\delta ^{ab}G_{\mu \nu }(q)`$ and performed the corresponding trace. The flavor content of the meson is not made explicit in eq. (12) as we have suppressed the flavor labels in the quark propagators. It is understood that the two propagators in eq. (12) are taken such as to account for the flavor quantum numbers of the considered meson. In the model that we will consider, the up and down quarks will be assumed to have equal current masses ($`m_0`$ in eq. (7)) and thus also identical propagator functions $`A(x)`$ and $`B(x)`$. For the light quarks, which should give rise to the familiar $`SU(3)_f`$ nonet, we are thus left with three representatives of each of the multiplets that are distinguished by their isospin number, $`I=0,\frac{1}{2},1`$ . We must also specify the meson angular momentum and parity. This is reflected by the Dirac and Lorentz decomposition of the meson vertex functions. This decomposition is known in the literature and here we follow Ref. . For the pseudoscalar channel ($`J^P=0^{}`$) we take

$$\mathrm{\Gamma }^{(P)}(p;P)=\gamma _5\left[\mathrm{\Gamma }_0^{(P)}(p;P)i\text{ / }P\mathrm{\Gamma }_1^{(P)}(p;P)i\text{ / }p\mathrm{\Gamma }_2^{(P)}(p;P)[\text{ / }P,\text{ / }p]\mathrm{\Gamma }_3^{(P)}(p;P)\right].$$
(13)
The decomposition for a scalar ($`J^P=0^+`$) meson reads

$$\mathrm{\Gamma }^{(S)}(p;P)=\mathrm{\Gamma }_0^{(S)}(p;P)i\text{ / }P\mathrm{\Gamma }_1^{(S)}(p;P)i\text{ / }p\mathrm{\Gamma }_2^{(S)}(p;P)[\text{ / }P,\text{ / }p]\mathrm{\Gamma }_3^{(S)}(p;P).$$
(14)
The vector ($`J^P=1^{}`$) channel involves eight scalar functions

$`\mathrm{\Gamma }_\mu ^{(V)}(p;P)`$ $`=`$ $`\left[\gamma _\mu {\displaystyle \frac{P_\mu \text{ / }P}{P^2}}\right]\left[i\mathrm{\Gamma }_0^{(V)}(p;P)+\text{ / }P\mathrm{\Gamma }_1^{(V)}(p;P)\text{ / }p\mathrm{\Gamma }_2^{(V)}(p;P)+i[\text{ / }P,\text{ / }p]\mathrm{\Gamma }_3^{(V)}(p;P)\right]`$ (17)
$`+\left[p_\mu {\displaystyle \frac{P_\mu pP}{P^2}}\right]\left[\mathrm{\Gamma }_2^{(V)}(p;P)+2i\text{ / }P\mathrm{\Gamma }_3^{(V)}(p;P)\right]`$
$`+\left[p_\mu {\displaystyle \frac{P_\mu pP}{P^2}}\right]\left[\mathrm{\Gamma }_4^{(V)}(p;P)+i\text{ / }P\mathrm{\Gamma }_5^{(V)}(p;P)i\text{ / }p\mathrm{\Gamma }_6^{(V)}(p;P)+[\text{ / }P,\text{ / }p]\mathrm{\Gamma }_7^{(V)}(p;P)\right].`$
In the axialvector channel we have two modes that are distinguished by their charge conjugation properties . For $`J^{PC}=1^{++}`$ the decomposition is

$`\mathrm{\Gamma }_\mu ^{(A)}(p;P)`$ $`=`$ $`\gamma _5\left[\gamma _\mu {\displaystyle \frac{P_\mu \text{ / }P}{P^2}}\right]\left[i\mathrm{\Gamma }_0^{(A)}(p;P)+\text{ / }P\mathrm{\Gamma }_1^{(A)}(p;P)\text{ / }p\mathrm{\Gamma }_2^{(A)}(p;P)+i[\text{ / }P,\text{ / }p]\mathrm{\Gamma }_3^{(A)}(p;P)\right]`$ (19)
$`+\gamma _5\left[p_\mu {\displaystyle \frac{P_\mu pP}{P^2}}\right]\left[\mathrm{\Gamma }_2^{(A)}(p;P)+2i\text{ / }P\mathrm{\Gamma }_3^{(A)}(p;P)\right],`$
while the $`J^{PC}=1^+`$ mode is decomposed as

$$\mathrm{\Gamma }_\mu ^{(\stackrel{~}{A})}(p;P)=\gamma _5\left[p_\mu \frac{P_\mu pP}{P^2}\right]\left[\mathrm{\Gamma }_1^{(\stackrel{~}{A})}(p;P)+i\text{ / }P\mathrm{\Gamma }_2^{(\stackrel{~}{A})}(p;P)i\text{ / }p\mathrm{\Gamma }_3^{(\stackrel{~}{A})}(p;P)+[\text{ / }P,\text{ / }p]\mathrm{\Gamma }_4^{(\stackrel{~}{A})}(p;P)\right].$$
(20)
In what follows we will omit the superscripts that label the spin and parity channels because these channels do not mix and there should hence be no confusion.

The primary object of this paper is to extract the bound state masses for the various flavor combinations and angular momentum channels. The traceology that is involved to project the Bethe–Salpeter equation (12) onto the Dirac components $`\mathrm{\Gamma }_i`$ is lengthy, but straightforwardly performed using the algebraic manipulation package FORM . The corresponding projection results in sets of coupled equations for the $`\mathrm{\Gamma }_i`$. After carrying out two of the three angular integrals analytically we are left with functions of the squared momenta $`p^2`$ and $`P^2`$ as well as the angle between $`p`$ and $`P`$: $`z=pP/\sqrt{p^2P^2}`$. The $`z`$–dependence is analyzed by an expansion in Chebyshev polynomials $`T_k`$

$$\mathrm{\Gamma }_i(p;P)=\underset{k=0}{\overset{N_{ch}1}{}}(i)^k\mathrm{\Gamma }_i^k(p^2;P^2)T_k(z).$$
(21)
Since the $`T_k`$ form an orthonormal set, we can project the equations for the Dirac components onto $`\mathrm{\Gamma }_i^k(p^2;P^2)`$. Finally, the $`k^2`$–integral in the Bethe–Salpeter equation (12) is implemented numerically as a matrix equation for the unknown $`\mathrm{\Gamma }_i^k(p_j^2;P^2)`$, $`p_j^2`$ being the discrete values of the momentum squared. The kernel, $`K`$ of that matrix parametrically depends on the meson momentum $`P^2`$. We solve that matrix equation as an eigenvalue problem by tuning the meson momentum to $`P^2=M^2`$, such that $`\mathrm{Det}(1K)=0`$. This yields the desired meson mass $`M`$.

For practical reasons we need to truncate the expansion (21) at a certain order, $`N_{ch}`$. As noted before this violates relativistic covariance. Fortunately the momentum partition dependence provides a measure for the degree at which covariance is violated. We can mitigate this degree by increasing $`N_{ch}`$. Numerically we furthermore scan the momentum partition dependence of the computed meson masses to a posteriori verify covariance for our solutions to the Bethe–Salpeter equations, see Appendix A. Unless otherwise stated we use $`N_{ch}=4`$.

### B Pseudoscalar leptonic decay constant

The solution to the Bethe–Salpeter equation not only yields the meson masses but also the meson quark vertex functions that can be used to compute meson properties. Here we will focus on the pseudoscalar decay constants $`f_\pi `$ and $`f_K`$. In order to calculate these, we first have to normalize the vertex functions $`\mathrm{\Gamma }(p;P)`$. The Bethe–Salpeter equation is a homogeneous equation, and thus needs an additional normalization condition. As mentioned previously, that condition is obtained from demanding the pole in the four–quark Green’s function to be unity. For equal momentum partitioning, (i.e. for $`\xi =1/2`$ only) it reads

$`2P_\mu =3{\displaystyle \overline{}d^4k}`$ $`\mathrm{Tr}\{\overline{\mathrm{\Gamma }}(k,P){\displaystyle \frac{S(k+P/2)}{P_\mu }}\mathrm{\Gamma }(k,P)S(kP/2)`$ (23)
$`+\overline{\mathrm{\Gamma }}(k,P)S(k+P/2)\mathrm{\Gamma }(k,P){\displaystyle \frac{S(kP/2)}{P_\mu }}\}`$
where the trace is over Dirac matrices. The conjugate vertex function $`\overline{\mathrm{\Gamma }}`$ is defined as

$$\overline{\mathrm{\Gamma }}(p,P)=C\mathrm{\Gamma }^T(p,P)C^1,$$
(24)
where use is made of the charge conjugation matrix, $`C=\gamma _2\gamma _4`$. The quark propagator derivatives are calculated by differentiating the quark Schwinger–Dyson equations (5) and (7) analytically and then numerically integrating the corresponding expressions.

The decay constants are finally obtained from the coupling of the axial current to the quark loop

$$f=\frac{3}{M^2}\overline{}d^4k\mathrm{Tr}\left\{\mathrm{\Gamma }(k,P)S(k+P/2)\gamma _5\text{ / }PS(kP/2)\right\},$$
(25)
where again the trace is over the Dirac matrices and the vertex functions are normalized according to eq. (23). Again we have not made explicit the flavor quantum numbers.

## IV Numerical Results and Discussion

We are now in the position to present our numerical results to the system of Schwinger–Dyson and Bethe–Salpeter equations. Even though the meson results are the primary target of our investigation it is fruitful to first consider the quark propagator functions obtained from the Schwinger–Dyson equations, as they enter the kernel of the Bethe–Salpeter equations.

### A Quark propagator functions

In Fig. 1 we show the quark propagator functions $`A(p^2)`$, $`B(p^2)`$ and $`M(p^2)=B(p^2)/A(p^2)`$ along the positive, real spacelike axis, $`p^2>0`$.

This is the numerical solution to the coupled equations (11), which we emphasize is the basis for the quark solutions for complex momenta. The solution clearly shows that dynamical chiral symmetry breaking is occurring: the mass function $`M(p^2)`$ attains a sizable non–zero value, even in the case that the bare quark mass $`m_0`$ is zero. This phenomenon has been extensively studied (see for example Refs. ), and is an example of genuinely non–perturbative behavior as dynamical mass generation cannot occur at any order in perturbation theory. Recall that the effective interaction (1) that enters the Schwinger–Dyson equations does not contain the perturbative UV behavior, rather it has an exponential damping at high momenta. This is manifested in the quark propagator functions as a sharp transition from the low momentum behavior to the bare values in the high momentum region. This transition occurs at about $`1\mathrm{G}\mathrm{e}\mathrm{V}`$.

The quark mass function $`M(p^2)`$ is also a useful quantity to discuss the changes of the quark propagator with the input parameters. Table I shows $`M(p^2=0)`$ for those parameter sets of $`\omega `$, $`D`$ and $`m_0`$ that we will later use to compute the meson masses.

As already mentioned the quark mass function indicates the extent to which dynamical chiral symmetry breaking occurs and plays the role of the mass of the constituent quarks within hadrons. We see that the quark mass function decreases as the effective scale $`\omega `$ increases and the strength $`D`$ decreases.

Since the main thrust of the paper is the use of the quark propagator functions evaluated at precisely the complex momenta squared for which they enter the Bethe–Salpeter equation, it is important to study the behavior of the mass function at least along the timelike axis, i.e. for negative real momenta squared. In Fig. 2 we present typical results for the mass function squared $`M^2(p^2)`$ in that regime.

Although the effective interaction (1) is simple, the timelike functions are distinctly non–trivial. Already from Fig. 2 we observe a significant variation with the current quark mass. This indicates that the specific forms of timelike quark propagator functions possess pronounced model and parameter dependences. Investigating the model dependence requires changing the effective interaction (1) and is not subject of the present study. For the parameter dependence, we fortunately find (see below) that it does not significantly effect the model predictions for the meson masses. There are divergences whenever $`A`$ and $`B`$ change sign at different $`p^2`$. Note that there are several particle–like poles ($`p^2=M^2(p^2)`$) where the quark goes on–shell. This is an important observation, since it is not at all clear whether a confined quark can be on–shell. Indeed, one might be tempted to conclude that with on–shell singularities and a subsequent lack of confinement, the model will not be able to give meson masses at all. As we will see, this is fortunately not the case.

### B Mesons

Turning to the mesons, we note that there are four model parameters, $`\omega ,D,m_u`$ and $`m_s`$ that we first have to fit to empirical data. To this end, we initially choose the pseudoscalar meson observables $`M_\pi ,M_K`$ and $`f_\pi `$. Then one parameter remains unconstrained by the pseudoscalar sector alone. However, the condition that the quark propagator function reflects dynamical chiral symmetry breaking, i.e. that $`M(p^2=0)0.5\mathrm{GeV}`$ leaves only a small window for the remaining choice, cf. table I. All other masses and decay constants are subsequently model predictions. The resulting model parameter and the predicted kaon decay constant $`f_K`$, that is unexpectedly well reproduced, are shown in Table II.

The subsequently predicted meson masses are shown in Tables III-VI. In all cases we have the inequalities $`M_{u\overline{u}}<M_{u\overline{s}}<M_{s\overline{s}}`$, where the subscript labels the flavor content. These relations just reflect the quark–antiquark picture that is implicit in the present Bethe–Salpeter approach.

Obviously both the pseudoscalar (table II) and vector mesons (table III) can be very well described within our model with the choice $`\omega 0.5\mathrm{GeV}`$. Our results agree with a previous analysis of the vector mesons based on an effective interaction which included the perturbative type term . This shows that such terms do not have a large effect on the meson masses, at least for the pseudoscalar and vector cases. Indeed, in the context of low–energy meson phenomenology we conclude that the logarithmic tail, and its associated renormalization represent an unnecessary obfuscation.

The situation for the scalar mesons (table IV) is not quite that clear. To begin with the particle data group  does not provide a clear picture in this channel but only quotes a wide range for the mass of the lowest scalar (0.4 – 1.2 GeV). More detailed studies of the pseudoscalar scattering amplitudes revealed that the assignment of the scalar meson nonet is not at all established . In particular, these mesons may not be simple quark–antiquark bound states but e.g. might contain sizable admixture of 2quark–2antiquark pairs . In that respect we might interpret our results as a quark–antiquark model prediction for scalar mesons. Our results suggest that such a picture is too simple for these mesons. One might also speculate that the adopted ladder approximation could be insufficient.

For the axialvector mesons we have two channels that are distinguished by their charge conjugation properties, cf. tables V and VI. The quark–antiquark pairs that are bound to axialvector modes with negative charge conjugation eigenvalue tend to be lighter than those with the positive eigenvalue but otherwise equal quantum numbers. Generally we find that our predictions are lower than the assignments made by the particle data group .

We recognize from our results that the model predictions change only slightly within the large range of considered model parameters. This confirms that meson static properties are not too sensitive to the conjectural parameter dependence of the timelike quark propagator functions. Presumably meson properties whose computation involves larger timelike will exhibit a stronger sensitivity.

Already from table V we observe that by increasing $`\omega `$ the predicted mass of the $`J^{PC}=1^+`$ meson with pion flavor quantum numbers approaches the empirical massIn all other channels the dependence on the parameters is much softer within the allowed window.. We therefore further increased $`\omega `$ according to the rules discussed above. For $`\omega 0.8GeV`$ we reproduced the empirical value for the mass in that channel. However, this happened at the expense of significantly lowering $`f_K`$ and loosing the proper description of the vector mesons. We recall that the parameter $`\omega `$ has a physical interpretation as the location of the maximum of the interaction. Thus $`\omega =0.8GeV`$ seems intuitively too large for low–energy hadron physics and an unsatisfactorily description of the $`0^{}`$ and $`1^{}`$ mesons comes without surprise.

For non–diagonal flavor structures such as $`u\overline{s}`$, charge conjugation actually is not a sensible quantum number and the corresponding axial vector mesons $`1^{++}`$ and $`1^+`$ may mix. In table VII we present the results obtained from the full calculation that combines the Dirac decompositions (19) and (20).

Since our Bethe–Salpeter formalism only yields the lowest mass eigenstate within a given channel, the results presented in table VII should be compared to those in table V. The tiny changes for the flavor diagonal mesons are numerical artifacts. Surprisingly the changes for the non–diagonal flavor structure are also only of the order 1%. This suggests only a small mixing between the $`1^{++}`$ and $`1^+`$ states with the flavor structure $`u\overline{s}`$ and the $`1^{++}`$ and $`1^+`$ channels represent good approximations to the actual eigenstates.

We can extend our model beyond the light flavors up, down and strange. The only modification is the increase of the current quark mass, $`m_0`$. In Fig. 3 we show how the meson masses increase as the (equal) quark masses are increased into the charm sector $`m_c=1.125\mathrm{GeV}`$.

Exactly the same numerical code is used to construct these solutions to the Schwinger–Dyson and Bethe–Salpeter equations as for the light flavors. Clearly seen is the smooth way the masses increase from the chiral limit ($`m_0=0`$) into the heavy quark sector ($`m_0=m_c`$). This represents a convincing indicator for the stability of our technique. The $`c\overline{c}`$-meson masses can be loosely extracted (table VIII) and the data are surprisingly well reproduced.

We did not expect to be able to describe such a heavy system with such a simple model, derived (and fitted) as it is from pion physics. The lack of the correct UV behavior for the gluon is seemingly at odds with the scales present. However, the present results suggest that the Bethe–Salpeter equation is capable of describing all the angular momentum states equally well in the charm quark sector.

## V Summary and Outlook

In this paper we have studied the low–lying mesons as quark–antiquark bound states in a covariant approach using an effective interaction. This interaction is characterized by gluon exchange with the gluon propagator being dressed by a Gaussian shape function. The interaction is completed by the quark–gluon vertex that we take to be the tree–level perturbative one. In this manner the rainbow–ladder approximation to the system of Schwinger–Dyson and Bethe–Salpeter equation accounts for chiral symmetry. With this effective interaction, we have then consistently treated this system of integral equations by precisely implementing the quark propagator functions that solve the Schwinger–Dyson equations into the Bethe–Salpeter equations. Once the effective interaction exceeds a certain strength, the Schwinger–Dyson equations exhibit dynamical chiral symmetry breaking and the pseudoscalar mesons emerge as would–be Goldstone bosons. We have then used observed properties of the pseudoscalar mesons to determine the model parameters. The kaon decay constant represents a model prediction. It turned out to be in good agreement with the empirical data. Furthermore our results for the vector meson masses match the experimental data. The situation in the scalar channel is less satisfying. As we solely consider the mesons as bound states of quark–antiquark pairs, it is not surprising that the mass eigenvalues increase with the strangeness content. On the other hand it is astonishing that for current quark masses, $`m_00.2\mathrm{GeV}`$, the lightest scalar mesons turn out to be heavier than the lightest vector mesons. When discussing these results it must be noted that the role of the scalar mesons is still under intense debate. In particular, the question whether they should indeed be considered as quark–antiquark bound states is not yet completely resolved. There are indications, see e.g. Ref.  and references therein, that the scalar meson masses should actually decrease with the strangeness content of these mesons. This can be understood if these mesons are considered as 2quark–2antiquark bound states in the sense of diquark–antidiquark systems .

As an outlook we mention that there is an elegant way to extend the present model to incorporate such degrees of freedom. The Bethe–Salpeter treatment can be straightforwardly extended to study bound states of diquark–antidiquark pairs, once a binding mechanism is established. This could either be achieved by a gluon exchange similar to eq. (1) or by quark exchange between a quark and a diquark. The latter approach has been intensively studied and the corresponding vertex is known from modeling baryon properties . It will also be interesting to see whether these additional degrees of freedom will also affect the mass predictions for the axialvector mesons that currently tend to be on the low side. Investigations in this direction are in progress.

Finally we would like to repeat that the stability of our treatment – measured by the actual momentum partition invariance – allows us to even describe mesons with charm quantum numbers. Without further modifications of the model parameter, except the corresponding current quark mass, our approach reproduces the mass eigenvalues of the flavor neutral states unexpectedly well.

## Acknowledgments

We thank J.C.R. Bloch, C. S. Fischer, D. J. Länge, P. Maris, M. Oettel, H. Reinhardt, S. M. Schmidt, and P. C. Tandy for valuable discussions. This work has been supported by DFG (Al-297/3–3, Al-297/3–4, We-1254/3–2, We-1254/4–2) and COSY (contract no. 41376610).

## A Relativistic covariance

As mentioned in section III the truncation of taking $`N_{\mathrm{ch}}`$ finite in eq. (21) violates relativistic covariance. Studying the dependence of our results on the momentum partitioning parameter $`\xi `$ we can provide a measure of these violations originating from the numerical method. The results for the pseudoscalar case are shown in Fig. 4.

The pion mass is extremely stable over the whole range $`\xi [0,1]`$, whereas for the kaon, the asymmetry between the quark masses must be compensated for by increasing the number of Chebyshev moments, $`N_{\mathrm{ch}}`$. Clearly, for the light pseudoscalars, the technique of inserting the quark propagators directly results in extremely stable results. In the case of the other mesons, the stability becomes less evident as one is studying heavier mass resonances (but with the same quarks). As an example, we display the axial–vector ($`J^{PC}=1^+`$) $`u\overline{s}`$ meson (Fig. 5). There is still a clear range of stability $`\xi [0.4,0.7]`$. Outside this range however, we notice that the technique does not simply break down (a mass can still be found) – it merely gets numerically more difficult to get trustworthy results. We assert the increasing dependence on $`\xi `$ with larger meson masses to the fact that with larger meson masses the calculation becomes more sensitive to the violent behavior seen in the timelike quark propagator functions. The larger the influence of that behavior to the calculation the higher we have to assume the number of Chebyshev moments for an adequate representation. Away from this timelike regime the quark propagator functions are smooth. This results in a good representation of the Bethe–Salpeter kernel and vertex functions with only a few Chebyshev moments.