text:
# Confinement Phenomenology in the Bethe-Salpeter Equation

## I Introduction

QCD models based on the Dyson–Schwinger equations \[DSEs\] provide an excellent tool in the study of nonperturbative aspects of hadron properties in QCD Roberts and Schmidt (2000). Such models can implement quark and gluon confinement Roberts and Schmidt (2000); Burden et al. (1992); Krein et al. (1992); Maris (1995), generate dynamical chiral symmetry breaking Atkinson and Johnson (1988); Roberts and McKellar (1990), and maintain Poincaré covariance. It is straightforward to implement the correct one-loop renormalization group behavior of QCD Maris and Roberts (1997), and obtain agreement with perturbation theory in the ultraviolet region. Provided that the relevant Ward identities are preserved in the truncation of the DSEs, the corresponding currents are conserved. Axial current conservation induces the Goldstone nature of the pions and kaons Maris et al. (1998); electromagnetic current conservation produces the correct hadronic charge without fine-tuning Oettel et al. (2000).

Previous work Maris and Roberts (1997); Maris and Tandy (1999) has shown this model to provide an efficient description of the light-quark pseudoscalar and vector mesons via two infrared parameters within the rainbow truncation of the DSE for solution of the dressed quark propagators coupled to the ladder approximation for the Bethe-Salpeter equation \[BSE\]. Furthermore, in impulse approximation, the elastic charge form factors of the pseudoscalars Maris and Tandy (2000) and the electroweak transition form factors of the pseudoscalars and vectors Maris and Tandy (2002); Ji and Maris (2001) are in excellent agreement with data. As a collorary, the strong decays of the vector mesons are well-described in impulse approximation without parameter adjustment Jarecke et al. (2002).

Nonperturbative solutions of the DSEs are almost always implemented in the Euclidean metric for technical reasons and this work is no exception. A resulting complication is that the solution of the BSE for meson bound states requires an analytic continuation in the meson momentum to reach the on-mass-shell point $`P^2=M^2`$. This causes the quark $`p^2`$ in the BSE to vary throughout a domain of the complex plane bounded by a parabola which widens in proportion to the mass $`M`$ of the meson. For low mass mesons, such as the $`\pi `$ and $`K`$ mesons, this excursion into the complex plane is simple to handle. Although the difficulties encountered in studies of ground state vector mesons are greater, they too may be overcome directly Maris and Tandy (1999). For mesons with masses greater than about $`1\mathrm{GeV}`$, these excursions into the complex plane are deep enough to encounter singularities in the rainbow DSE solutions for the quark propagators. As a result, little is known about the implementation of the BSE in this situation or even whether a solution is well-defined. For this reason, DSE-BSE model studies of mesons have for the most part been restricted to light mesons. An exception is the study of Ref. Jain and Munczek (1993) in which the Ansatz adopted for the quark propagators in the non-Euclidean domain precluded singular behavior.

It has long been known that the rainbow approximation to the fermion DSE produces complex conjugate singularities on the timelike half of the complex-momentum plane ($`\mathrm{}(p^2)<0`$) in QED Fukuda and Kugo (1976); Atkinson and Blatt (1979). The phenomenon has been studied in detail within QED Maris and Holties (1992); Maris (1994) and in an Abelian-like model of QCD Stainsby and Cahill (1992). In these early studies the singularities were attributed to artifacts of the rainbow approximation. A connection between the generation of complex conjugate singularities and the possible absence of real mass poles was noticed in a QED study Fukuda and Kugo (1976). Later studies within QCD Burden et al. (1992); Krein et al. (1992); Gribov (1999) raised suggestions that complex singularities may not be artifacts of the approximation at all, but rather may be a property of the full theory, possibly related to confinement and the absence of real quark mass poles. Beyond the rainbow approximation, very little is known about DSE solutions and their singularity structure. These questions and the influence of approximations have yet to be completely resolved. However they are being addressed; for example, in one simple model calculation, the inclusion of the full dressing of the quark-gluon vertex in the DSE leads to the disappearence of all quark singularities in the complex plane save essential singularities at infinity Burden et al. (1992).

It has been well-established by Dyson-Schwinger studies that in order to obtain an acceptable amount of dynamical chiral symmetry breaking in agreement with empirical observations, the effective rainbow DSE kernel must have a significant infrared strength. This in turn tends to move mass poles away from the real timelike $`p^2`$ axis, at least for a domain characterized by several GeV. This suggests that quark propagators may be successfully modeled in terms of parametrizations based on entire functions Roberts (1996) or at least based on functions with singularities away from the real axis such as complex-conjugate poles.

If such an approach to modeling QCD is to be extended to include $`\overline{q}q`$ and $`\overline{q}Q`$ meson spectra beyond 1 or 2 GeV, then it is necessary to determine how solutions of the BSE are defined when complex singularities in the quark propagators may be encountered. This issue is addressed in the following using a simple parmetrization of the quark propagator. Although the method developed herein is first applied to the BSE, it has applications beyond meson spectroscopy. For example, in production processes at moderate momentum transfers, the use of complex mass pole representations for confined quarks has been suggested as a means to handle the sizable timelike momentum carried by the quarks in these processes Ahlig et al. (2001). One deals with these issues also in applications of the DSE-BSE approach to electromagnetic elastic Maris and Tandy (2000) and transition form factors Maris and Tandy (2002); Ji and Maris (2001) of light mesons at high-momentum transfer.

To simplify the analysis, we adopt a parameterization of the dressed quark propagator in terms of complex conjugate mass poles. In the standard approach to the meson BSE, two of the four integration variables are integrated numerically and it is during the integration of these variables that poles may be encountered as a result of analytic continuation. In this article, a clear implementation, based on Feynman integral techniques, is introduced to map the four-dimensional integration in the BSE onto a one-dimensional domain where numerical methods of evaluating Cauchy integrals are well-established. We find that when a complex conjugate mass pole parametrization of the quark propagators is employed, the would-be $`\overline{q}q`$ decay width of the meson generated from a particular mass pole is exactly cancelled by the contribution of the conjugate partner pole. As a result, the meson mass $`M`$ remains real valued, even above the pseudo-threshold $`M^2>\mathrm{}(4m_1^2)`$ where $`m_1`$ is the lowest complex mass in the propagator representation.

The primary limitation of the method developed here is that it relies on an explicit parameterization of dressed quark propagators in terms of pairs of complex conjugate mass poles. However the limitation is offset by the finding that such a representation can provide excellent fits to both DSE solutions (as is shown in Sec. V) and also to recent lattice-QCD data Detmold (2002). Of course, use of such a phenomenological representation will not exactly preserve the dynamical relation between vertices and propagators necessary to satisfy the Ward Takahashi identities of QCD. In particular, without the axial vector Ward Takahashi identity, the chiral-limit pseudoscalar meson solutions will not obey Goldstone’s theorem Maris et al. (1998). As a result, the toy model introduced in Sec. V would not be recommended for studies of light pseudoscalar mesons; chiral-symmetry-preserving DSE approaches have already been well-established Maris and Roberts (1997); Maris and Tandy (1999) for that sector. Rather, the present work explores techniques for extending the present DSE-based approach to mesons with masses above 1 GeV where one expects to encounter complex conjugate singularities.

The article is organized as follows. In Sec. II, our notation for the ladder BSE for meson states is introduced, as is the complex mass pole representation of the quark propagators which will be employed. In Sec. III, a simplified BSE integral equation is used to illustrate the technique for analytic continuation of the BSE. The fact that the resulting bound state mass will be real is demonstrated explicitly. In Sec. IV, the method is extended to the general BSE for a vector meson. A simple model for the ladder BSE kernel is introduced in Sec. V and numerical results are presented. A discussion of this approach and outlook is presented in Sec. VI.

## II The Bethe-Salpeter Equation

In the following the Euclidean metric is employed. The scalar product of two four vectors is then $`ab=_{\mu =1}^4a_\mu b_\mu `$, and $`a^2>0`$ for a spacelike 4-vector $`a_\mu `$. The Dirac $`\gamma `$-matrices are Hermitian $`\gamma _\mu ^{}=\gamma _\mu `$ and obey the anti-commutation relation $`\{\gamma _\mu ,\gamma _\nu \}=2\delta _{\mu \nu }`$. The dressed quark propagator $`S(p)`$ and meson Bethe-Salpeter (BS) amplitude $`\mathrm{\Gamma }(p;P)`$ are solutions of the renormalized DSE,

$`S(p)^1`$ $`=`$ $`Z_2i/p+Z_4m_q(\mu )`$ (1)
$`+Z_1{\displaystyle _q^\mathrm{\Lambda }}g^2D_{\mu \nu }(pq){\displaystyle \frac{\lambda ^i}{2}}\gamma _\mu S(q)\mathrm{\Gamma }_\nu ^i(q,p),`$
and the BSE,

$`\mathrm{\Gamma }(p;P)`$ $`=`$ $`{\displaystyle _q^\mathrm{\Lambda }}K(p,q;P)\left(S(q_+)\mathrm{\Gamma }(q;P)S(q_{})\right),`$ (2)
Here $`D_{\mu \nu }(k)`$ is the renormalized dressed-gluon propagator, $`\mathrm{\Gamma }_\nu ^i(q,p)`$ is the renormalized dressed quark-gluon vertex, $`\mathrm{\Gamma }(p;P)`$ is the BS amplitude for a quark-antiquark bound-state meson, and $`K(q,p;P)`$ is the renormalized two-particle irreducible $`q\overline{q}`$ scattering kernel. We consider equal mass quarks with momenta $`q_\pm =q\pm P/2`$ where $`q`$ is the quark-antiquark relative momentum, and $`P=q_+q_{}`$ is the meson momentum which satisfies $`P^2=M^2`$ where $`M`$ is the meson mass.

In Eq. (2) the double contraction of Dirac indices has been denoted by $``$. If the Dirac indices of the elements were reinstated in the BSE and written as parenthetic Roman letters $`(a,b,\mathrm{})`$, then Eq. (2) would appear as

$`\mathrm{\Gamma }_{(ab)}(p;P)`$ $`=`$ $`{\displaystyle _q^\mathrm{\Lambda }}K_{(ab;cd)}(p,q;P)S_{(de)}(q_+)`$ (3)
$`\times \mathrm{\Gamma }_{(ef)}(q;P)S_{(fc)}(q_{}).`$
Use of the double contraction operator $``$ allows us to suppress Dirac indices, and simplify the forms of the equations appearing herein.

In Eqs. (1) and (2) $`_k^\mathrm{\Lambda }^\mathrm{\Lambda }d^4k/(2\pi )^4`$ denotes a translationally invariant ultraviolet regularization of the momentum space integral with mass-scale $`\mathrm{\Lambda }`$. Lorentz covariance entails that the solution of Eq. (1) is of the form $`S(p)^1=i/pA(p^2)+`$ $`B(p^2)`$. The Lorenz scalars $`A(p^2)`$ and $`B(p^2)`$ are renormalized at spacelike $`p^2=\mu ^2`$ such that $`A(\mu ^2)=1`$ and $`B(\mu ^2)=m_q(\mu )`$, where $`m_q(\mu )`$ is the renormalized current quark mass. After renormalization, one is free to remove the regularization by taking the limit $`\mathrm{\Lambda }\mathrm{}`$.

Recently, a successful model has been developed that provides an excellent description of the masses, decays and other properties of the light pseudoscalar and vector mesons Maris and Roberts (1997); Maris and Tandy (1999). The model consists of the rainbow-ladder truncations for the quark DSE (1) and meson BSE (2) and the use of an effective $`\overline{q}q`$ interaction constrained to coincide with perturbative QCD in the ultraviolet domain and containing a phenomenological infrared behavior. In this particular model, the DSE kernel in rainbow truncation is

$$Z_1g^2D_{\mu \nu }(k)\mathrm{\Gamma }_\nu ^i(q,p)4\pi \alpha _{\mathrm{eff}}(k^2)D_{\mu \nu }^{\mathrm{free}}(k)\gamma _\nu \frac{\lambda ^i}{2},$$
(4)
where $`D_{\mu \nu }^{\mathrm{free}}(k=pq)`$ is the free gluon propagator in Landau gauge, and $`\alpha _{\mathrm{eff}}(k^2)`$ is the effective running coupling. The ladder truncation of the BSE kernel is

$$K(p,q;P)4\pi \alpha _{\mathrm{eff}}(k^2)D_{\mu \nu }^{\mathrm{free}}(k)\frac{\lambda ^i}{2}\gamma _\mu \frac{\lambda ^i}{2}\gamma _\nu ,$$
(5)
where $`k=pq`$ and $``$ denotes the direct (or Cartesian) product of matrices. This truncation scheme is self-consistent in that it ensures that the dressing of the quark-antiquark vector and axial-vector vertices generated by Eqs. (4) and (5) satisfy their respective Ward-Takahashi identities. This feature is important since the axial-vector Ward-Takahashi identity guarantees that, in the chiral limit of zero current quark mass, the ground state pseudoscalar meson bound states are the massless Goldstone bosons arising from dynamical chiral symmetry breaking in accordance with the Goldstone’s theorem Maris and Roberts (1997); Maris et al. (1998). The vector Ward-Takahashi identity guarantees electromagnetic current conservation for mesons and nucleons if the impulse approximation is used to describe the current in terms of quarks Maris and Tandy (2000); Oettel et al. (2000). The rainbow-ladder truncation is particularly suitable for the flavor octets of light pseudoscalar and vector mesons where higher-order contributions to the quark-gluon skeleton-graph expansion have significant cancellations Bender et al. (1996, 2002) and so may be safely neglected.

With the ladder truncation from Eq. (5), the BSE for a light-quark vector meson is

$`\lambda (P^2)\mathrm{\Gamma }_\rho (p;P)=`$
$`{\displaystyle _q^\mathrm{\Lambda }}K_{\mu \nu }^L(k)\gamma _\mu S(q_+)\mathrm{\Gamma }_\rho (q;P)S(q_{})\gamma _\nu ,`$
with $`k=pq`$ and

$$K_{\mu \nu }^L(k)=\frac{16\pi }{3}\alpha _{\mathrm{eff}}(k^2)D_{\mu \nu }^{\mathrm{free}}(k).$$
(7)
Here $`\mathrm{\Gamma }_\rho `$ is a transverse four vector. We will refer to the vector meson under study here as the $`\rho `$ meson although, due to a number of simplifications made for illustration purposes, the particular model solutions considered here should not be taken as physical representations of the $`\rho `$ meson.

A linear eigenvalue $`\lambda (P^2)`$ has been introduced into Eq. (LABEL:BSE\_ladder) so that it will yield solutions over a continuous range of $`P^2`$. A physical bound state of mass $`M_\rho `$ corresponds to an eigenvalue of $`\lambda (P^2=M_\rho ^2)=1`$. To find such a solution, it is clear that Eq. (LABEL:BSE\_ladder) must be analytically continued from Euclidean space where $`P^2>0`$ to Minkowski space where $`P^2<0`$ and the on-mass-shell condition $`P^2=M_\rho ^2`$ may be realized. In terms of the Euclidean four momentum $`P_\mu `$, such an analytic continuation would correspond to $`P_4=`$ $`i\sqrt{\stackrel{}{P}^2+M_\rho ^2}`$, from which it follows that quark momenta $`q_\pm ^2=q^2`$ $`M_\rho ^2/4\pm i\sqrt{q^2}M_\rho z`$ that enter the BSE are complex functions of the integration variables $`q^2>0`$ and the direction cosine $`1z+1`$. Therefore, in order to carryout such an analytic continuation of the BSE (LABEL:BSE\_ladder), the dressed quark propagators $`S(q_\pm ^2)`$ must be known in the complex-momentum plane within a parabolic region having the negative real (timelike) point $`M_\rho ^2/4`$ as the apex and extending symmetrically about the real axis.

Studies within the rainbow-ladder truncation Maris and Tandy (1999) find that the non-analytic points of $`u`$\- and $`d`$-quark propagators that occur nearest to the origin $`q^2=0`$ are located at $`q^2=0.207\pm i0.331\mathrm{GeV}^2`$ while those of the $`s`$-quark propagator are found at $`q^2=0.376\pm i0.602\mathrm{GeV}^2`$ Jarecke et al. (2002). This finding also agrees with qualitative observations from an earlier study Frank and Roberts (1996). From the positions of these non-analytic points one can define a critical mass $`M_c=\sqrt{P^2}>0`$ above which these singular points enter the domain of integration for the BSE. For the model of LABEL:Maris:1999nt, these critical masses are $`M_c=1.09`$ GeV for $`\overline{u}u`$ mesons and $`M_c=1.47`$ GeV for $`\overline{s}s`$ mesons. Clearly, the non-analytic points lie outside the domains of integration for ladder BSE calculations of ground state pseudoscalar and vector mesons, but may lie within the integration domain for heavier mesons. This observation provides one of the main reasons that many BSE studies of mesons focus primarily on light meson states.

The appearence of poles in the domain of integration of the BSE raises two questions: 1) does the analytic continuation of the BSE to timelike momenta $`P^2<0`$ yield physically meaningful quantities, and if so, 2) by what methods does one obtain the solutions of this singular integral equation? The goal of this article is to explore these questions. To facilitate an initial investigation, it is convenient to adopt a parametrization of the quark propagator $`S(p)`$ as a sum of $`N`$ pairs of complex conjugate mass poles,

$$S(p)=\underset{n=1}{\overset{N}{}}\left\{\frac{z_n}{i/p+m_n}+\frac{z_n^{}}{i/p+m_n^{}}\right\},$$
(8)
where $`m_n`$ are complex-valued mass scales and $`z_n`$ are complex coefficients. Comparison with the general form $`S(p)^1=i/pA(p^2)+`$ $`B(p^2)`$, allows identification of the amplitudes $`A(p^2)`$, $`B(p^2)`$, and the quark running-mass $`M(p^2)=B(p^2)/A(p^2)`$. Equivalently, comparison with the general form $`S(p)=i/p\sigma _V(p^2)+`$ $`\sigma _S(p^2)`$ allows identification of the scalar and vector amplitudes of the propagator, $`\sigma _S(p^2)`$ and $`\sigma _V(p^2)`$ respectively. These are needed for the BSE. For example, the quark propagator of Eq. (8) has a scalar amplitude of the form

$$\sigma _S(p^2)=\underset{n=1}{\overset{N}{}}\left\{\frac{z_nm_n}{p^2+m_n^2}+\frac{z_n^{}m_n^{}}{p^2+m_{n}^{}{}_{}{}^{2}}\right\},$$
(9)
and a vector amplitude $`\sigma _V(p^2)`$ which is obtained by removing each of the $`m_n`$ and $`m_n^{}`$ from the numerators.

It is shown in Sec. V that this form of quark propagator with $`N=3`$ provides an excellent fit to DSE solutions that employ realistic effective $`q\overline{q}`$ interactions. Furthermore, this form has been used recently to provide a parametric fit of the chiral-quark propagator Detmold (2002) obtained from lattice-QCD simulationsBowman et al. (2002); Bowman (2002).

## III Analytic Continuation: A Simple Example

A general study of the analytic continuation of the BSE is quite involved and complicated by the Dirac algebra necessary to describe quark-antiquark states. It is prudent to first proceed by simplifying several of these complications and focus on the issue of analytic continuation of the eigenvalue $`\lambda (P^2)`$ and the associated integrations over the complex-conjugate poles in the quark propagator. These restrictions are lifted in Sec. IV where the approach introduced here is extended to the general BSE.

In the following, assume that the quark propagator is well-represented by a single pair of complex conjugate poles (that is, $`N=1`$ in Eq. (8)). Furthermore, assume that the vector-meson BS amplitude may be represented by the single Dirac amplitude $`\mathrm{\Gamma }_\rho (p;P)=(\delta _{\rho \alpha }P_\rho P_\alpha /P^2)\gamma _\alpha V(q;P)`$, where $`P_\rho `$ is the four momentum of the bound-state vector meson. After these simplifications, the BSE (LABEL:BSE\_ladder) may be converted to a scalar eigenvalue equation for the amplitude $`V(q;P)`$ by multiplication from the left by $`\gamma _\rho `$ followed by a trace over the Dirac indices to obtain

$`\lambda (P^2)V(p;P)=`$
$`{\displaystyle \underset{a,b=\pm 1}{}}{\displaystyle _q^\mathrm{\Lambda }}{\displaystyle \frac{f_{ab}(p,q,P)V(q;P)}{(q_+^2+\mu ^2+ia\mathrm{\Delta }^2)(q_{}^2+\mu ^2+ib\mathrm{\Delta }^2)}}.`$ (10)
Observe that the kernel of this equation is bilinear in the propagator amplitudes $`\sigma _S(q_\pm ^2)`$ and $`\sigma _V(q_\pm ^2)`$ but only the denominators of these functions are displayed explicity in Eq. (10). The remaining factors from the quark propagators and effective $`q\overline{q}`$ interaction have been collected into the scalar function $`f_{ab}`$ whose detailed form is unimportant to the present discussion but may be derived following the discussion of Sec. IV. In addition, the denominators have been written in terms of $`\mu ^2=m_R^2m_I^2`$ and $`\mathrm{\Delta }^2=2m_Rm_I`$, where $`m=m_R+im_I`$ is the complex mass scale parameter in the quark propagator in Eq. (9). This form allows a parametrization of the square of the quark mass $`m^2`$ and its complex conjugate $`m^2`$ in terms of the values $`a=\pm 1`$.

In principle, the $`q`$-dependence of the integrand in Eq. (10) is only known once the solution $`V(q;P)`$ has been obtained. For the present, consider the artificial situation in which the solution of the BSE is known to be of the form $`V=g\varphi (q^2)U(z_q)`$ where $`z_q=qP/\sqrt{q^2P^2}`$ and $`\varphi (q^2)`$ and $`U(z)`$ are known functions. In this case, the constant $`g`$ may be eliminated from the BSE and following an integration over $`d^4p`$ one obtains,

$$\lambda (P^2)=\underset{a,b=\pm 1}{}\lambda _{ab}(P^2),$$
(11)
and

$$\lambda _{ab}(P^2)=\frac{d^4q}{(2\pi )^4}\frac{F_{ab}(q,P)}{(q_+^2+\mu ^2+ia\mathrm{\Delta }^2)(q_{}^2+\mu ^2+ib\mathrm{\Delta }^2)},$$
(12)
where $`F_{ab}(q,P)`$ is a non-singular, Lorentz invariant function of $`q`$ and $`P`$. The objective is then to evaluate the four integrals $`\lambda _{ab}(P^2)`$ and show how they may be analytical continued to negative values of $`P^2`$. This is carried out explicitly in the remainder of this section.

When $`P^2`$ is continued to negative values, the quark momenta in Eq. (12) satisfy $`q_\pm ^2=(q_{}^2)^{}`$. The two factors in the denominator of $`\lambda _{a,a}(P^2)`$ (that is, $`a=b`$ in Eq. (12)) will thus be complex conjugates. As a result, for $`P^2`$ sufficiently timelike, there will always be some value of the integration variable $`q_\mu `$ for which both factors in the denominator of $`\lambda _{a,a}(P^2)`$ vanish simultaneously. The result is a non-analytic point in $`\lambda _+(P^2)`$ and $`\lambda _+(P^2)`$. Such behavior would usually be indentified with the decay threshold for the vector meson into an asymptotically-free quark and antiquark pair. However, in the following, it is shown that when the integrals in Eq. (12) are added together the imaginary parts cancel exactly thereby eliminating the possibility of quark-antiquark production thresholds. In contrast to this, the two factors in the denominator of $`\lambda _{aa}(P^2)`$ never vanish simultaneously and consequently there are no corresponding non-analytic points in either $`\lambda _{++}(P^2)`$ or $`\lambda _{}(P^2)`$.

To analyze the behavior in detail, it is advantageous to restructure the integral in Eq. (12) so that the analytic continuation from the Euclidean domain $`P^2>0`$ to the timelike domain $`P^2<0`$ is straightforward. For values of $`P^2>0`$, one may combine the two factors in the denominator using the Feynman parametrization,

$$\frac{1}{AB}=_0^1𝑑\alpha 𝑑\beta \frac{\delta (\alpha +\beta 1)}{(\alpha A+\beta B)^2}$$
(13)
and obtain

$`\lambda _{ab}(P^2)`$ $`=`$ $`{\displaystyle _1^{+1}}{\displaystyle \frac{dy}{2}}{\displaystyle \frac{d^4q}{(2\pi )^4}\frac{F_{ab}(q,P)}{D(q,P,y)^2}},`$ (14)
where $`y=\alpha \beta `$, and

$$D(q,P,y)=y(qP+i\frac{ab}{2}\mathrm{\Delta }^2)+(q^2+\frac{P^2}{4}+\mu ^2+i\frac{a+b}{2}\mathrm{\Delta }^2).$$
(15)
One then eliminates the term linear in $`q`$ in the denominator by a shift of the origin $`q_\nu q_\nu yP_\nu /2`$. The resulting denominator depends on $`P^2`$, $`q^2`$, and $`y`$ but not on $`qP`$. After this shift, the function $`F_{ab}(q,P)`$ in the numerator still depends $`qP`$. All other variables save $`q^2`$ and $`y`$ can be integrated out yielding,

$`\lambda _{ab}(P^2)`$ $`=`$ $`{\displaystyle _1^{+1}}𝑑y{\displaystyle _0^{\mathrm{}}}𝑑q^2{\displaystyle \frac{q^2\stackrel{~}{F}_{ab}(q^2,y;P^2)}{[q^2q_{ab}^2(P^2,y)]^2}},`$ (16)
where

$$q_{ab}^2(P^2,y)=\frac{P^2}{4}(1y^2)\mu ^2\frac{i}{2}\left(y(ab)+a+b\right)\mathrm{\Delta }^2,$$
(17)
is independent of $`q^2`$, and

$$\stackrel{~}{F}_{ab}(q^2,y;P^2)=\frac{1}{16\pi ^3}_1^{+1}𝑑z_q\sqrt{1z_q^2}F_{ab}(q\frac{y}{2}P,P).$$
(18)
Here $`z_q`$ is the direction cosine defined by $`qP=\sqrt{q^2P^2}z_q`$. The singular behavior due to the propagators is now a single second-order pole in the variable $`q^2`$. Its location $`q_{ab}^2(P^2,y)`$ is given in Eq. (17) and is not in the path of integration for spacelike $`P^2>0`$; the integral is regular in Euclidean space.

The $`y`$ integration is performed after the integration over $`q^2`$. Hence, for a given value of $`y`$, the $`q^2`$-integration domain is the positive real axis and the integrand has a pole in the complex $`q^2`$ plane whose location $`q_{ab}^2(P^2,y)`$ moves with $`P^2`$. The location has a negative real part for $`P^2>0`$ and remains negative as $`P^2`$ passes through zero and becomes moderately timelike. Only when $`P^2<4\mu ^2`$ will the real part of the pole location become positive and have a chance of impinging upon the domain of the $`q^2`$-integration. For that to happen, the imaginary part of $`q_{ab}^2(y)`$ would have to vanish. As anticipated above, Eq. (17) shows that this never happens when $`a=b`$; for $`\lambda _{++}(P^2)`$ and $`\lambda _{}(P^2)`$ the two factors in the denominator of Eq. (12) are never zero simultaneously. It is therefore trivial to carryout these integrations. However, for $`a=b`$, the imaginary part of $`q_{ab}^2(y)`$ can become zero when the Feynman variable $`y=0`$ and care must be taken to properly handle such singular integrals.

As a specific example, the trajectory of the pole location $`q_+^2(y)`$ is shown in Fig. 1. It begins at $`(\mu ^2,\mathrm{\Delta }^2)`$ when $`y=1`$, moves along the dashed curve as $`y`$ increases, crossing the real axis when $`y=0`$, and terminates at $`(\mu ^2,\mathrm{\Delta }^2)`$ when $`y=1`$. For $`y<0`$, the calculation of the integral is straightforward. However for $`y>0`$, the pole in $`q^2`$ has crossed the real axis and the analytic continuation of the integral $`\lambda _+(P^2)`$ is defined by continuous deformation of the contour in a manner that ensures the pole does not cross the contour of integration. This is illustrated in Fig. 1, for the deformed contour corresponding to $`y+1/2`$. Deformation of a integration contour in this manner is equivalent to adding the pole residue to the result of integration along the real $`q^2`$ axis. If the residue term were neglected, one would obtain a discontinuous result due to the crossing of a branch cut onto the wrong sheet.

For the integral $`\lambda _+(P^2)`$, the pole trajectory traces out the same path as in Fig. 1 but in the opposite direction; that is, $`y>0`$ corresponds to the upper half $`q^2`$-plane and the pole residue will be of opposite sign from $`\lambda _+(P^2)`$. After analytic continuation to $`P^2<0`$, the results for the four integrals can be expressed as

$`\lambda _{ab}(P^2)=`$ $`{\displaystyle _1^{+1}}𝑑y{\displaystyle _0^{\mathrm{}}}𝑑s{\displaystyle \frac{s\stackrel{~}{F}_{ab}(s,y;P^2)}{[sq_{ab}^2(y,P^2)]^2}}`$ (19)
$`+{\displaystyle \frac{ab}{2}}\mathrm{\hspace{0.33em}2}\pi i\mathrm{\Theta }({\displaystyle \frac{P^2}{4}}\mu ^2)`$
$`\times {\displaystyle _0^1}dy{\displaystyle \frac{d}{ds}}\left(s\stackrel{~}{F}_{ab}(s,y;P^2)\right)|_{s=q_{ab}^2},`$
where $`\mathrm{\Theta }(x)`$ is the Heaviside step function.

In the timelike region, inspection of Eq. (19) shows that individual integrals $`\lambda _{ab}(P^2)`$ are complex valued and satisfy

$`[\lambda _{a,b}(P^2)]^{}`$ $`=`$ $`\lambda _{a,b}(P^2).`$ (20)
The imaginary parts cancel in pairs thereby producing a real result for the eigenvalue $`\lambda (P^2)=_{a,b=\pm }\lambda _{ab}(P^2)`$. Notice that this does not mean that one may neglect the addition of residue terms required by the analytic continuation of the integrals. Rather, when $`\mathrm{\Delta }^2>0`$ the residue (second) term of Eq. (19) is not purely imaginary and the sum of the real parts of the residues are required for the final result for $`\lambda (P^2)`$. Of course, these residue terms appear only for $`P^2`$ more timelike than the pseudo-threshold point $`P^2=4\mu ^2`$, that is, for bound state masses $`M_\rho >2\mu `$.

Before extension of this treatment to the general problem, consider the connection between the above discussion of analytic continuation and the usual case encountered above physical two-body thresholds. In particular, consider the limiting case of infinitesimal $`\mathrm{\Delta }`$ for the single integral $`\lambda _{}(P^2)`$. This corresponds to the standard case of employing elementary propagators with the conventional boundary conditions of a negative infinitesimal $`iϵ`$ added to the real mass. Above the two-body threshold $`P^2<4\mu ^2`$, there would be $`\overline{q}q`$ production described by the residue term in $`\lambda _{}(P^2)`$. It follows that a bound state calculation, based on such an integral alone, would result in a complex mass describing an unstable meson with a finite $`\overline{q}q`$ decay width. However, in the present case the presence of conjugate quark poles ensures that the analytic continuation introduces both positive and negative would-be widths that cancel exactly. The resulting BSE eigenvalue $`\lambda (P^2)`$ in Eq. (10) will be purely real and therefore bound state masses defined by $`\lambda (P^2=M_\rho ^2)=1`$ will be real and will describe a meson with no $`\overline{q}q`$ decay width.

## IV Generalization to the Bethe-Salpeter Equation

The technique described in the previous section is generalized to allow solutions to the homogeneous BSE for an arbitrary quark-antiquark bound state. To be specific, we consider an equal quark-antiquark vector meson bound state. One can write the most-general BS amplitude for a massive vector meson as a sum of eight Dirac covariants and Lorentz-scalar amplitudes $`F_\alpha (q;P)`$,

$`\mathrm{\Gamma }_\rho (q;P)`$ $`=`$ $`(\delta _{\rho \mu }{\displaystyle \frac{P_\rho P_\mu }{P^2}})(\gamma _\mu F_1+q_\mu \overline{)}qF_2`$ (21)
$`+q_\mu \overline{)}PF_3+\gamma _5ϵ_{\mu \alpha \nu \beta }\gamma _\alpha q_\nu P_\beta F_4`$
$`+q_\mu F_5+\sigma _{\mu \nu }q_\nu F_6+\sigma _{\mu \nu }P_\nu F_7`$
$`+q_\mu \sigma _{\alpha \beta }q_\alpha P_\beta F_8).`$
In a more concise form this is

$$\mathrm{\Gamma }_\rho (q;P)=\underset{\alpha =1}{\overset{8}{}}\mathrm{\Lambda }_\rho ^\alpha (q,P)F_\alpha (q;P),$$
(22)
where $`\mathrm{\Lambda }^\alpha `$ are the Dirac covariants in Eq. (21); for example, $`\mathrm{\Lambda }_\rho ^{\alpha =1}=`$ $`(\delta _{\rho \mu }P_\rho P_\mu /P^2)\gamma _\mu `$. When all indices of a product of two Dirac covariants are contracted in the manner corresponding to the usual trace operation, one obtains

$`\mathrm{tr}_\mathrm{D}\left[\mathrm{\Lambda }_\mu ^\alpha (q,P)\mathrm{\Lambda }_\mu ^\beta (q,P)\right]`$ $`=`$ $`\mathrm{\Lambda }_\mu ^\alpha (q,P)\mathrm{\Lambda }_\mu ^\beta (q,P),`$ (23)
$`=`$ $`t_{\alpha \beta }(q,P),`$
where $`t_{\beta \alpha }(q,P)`$ is a Lorentz scalar function of $`q^2`$, $`P^2`$ and $`qP`$, and we have used the double-contraction product $``$ notation introduced in Eq. (2).

Of course, one is free to choose any other set of Dirac covariants $`\mathrm{\Lambda }_\rho ^\alpha (q,P)`$ as long as they are complete. In particular, one is free to choose a set of covariants for which $`t_{\alpha \beta }\delta _{\alpha \beta }`$ in which case the resulting set of coupled integral equations may be of a simpler form than those obtained herein. However, in practice the behavior of such covariants can be poorly behaved near $`q=0`$, and can lead to numerical difficulties Maris and Tandy (1999). Therefore, the set of covariants given by Eq. (21) is preferable and the lack of orthogonality is easily accommodated by $`t_{\alpha \beta }(q,P)`$.

In the previous section, it was shown that when the $`q^2`$ and $`qP`$ dependence of the BS amplitude $`\mathrm{\Gamma }_\mu (q;P)`$ is known, the analytic continuation of the BSE may be carried out in a straight-forward manner. In practice, the functional form of the BS amplitude is known only after the solution to the BSE has been obtained. The approach outlined in the previous section may be generalized by expansion of the functional dependence of the BS amplitudes $`F_\alpha (q;P)`$ in terms of complete sets of known functions of $`q^2`$ and $`z_q=qP/\sqrt{q^2P^2}`$. Consider the expansion

$$F_\alpha (q;P)=\underset{s,i=0}{\overset{\mathrm{}}{}}U_i(z_q)\varphi _s(q^2)g_{si}^\alpha (P^2),$$
(24)
where $`U_i(z_q)`$ are Chebyshev polynomials of the second kind satisfying the ortho-normality condition

$$\frac{2}{\pi }_1^1𝑑z\sqrt{1z^2}U_i(z)U_j(z)=\delta _{ij},$$
(25)
and $`\varphi _s(q^2)`$ are a set of ortho-normal functions, that are specified later, and that satisfy

$$_0^{\mathrm{}}𝑑q^2\varphi _s(q^2)\varphi _t(q^2)=\delta _{st}.$$
(26)
The required number of basis terms will be determined by convergence of the BSE solution with respect to the meson observables of interest. For example, one particular choice of basis $`\{\varphi _s(q^2)\}`$ which proves to be efficient for observables dominated by infrared physics may be inefficient for observables sensitive to the ultraviolet. In the present work, we take convergence of the meson mass as the criterion.

For brevity, we define the operator $`𝒪_{si}(p)`$ which projects out the expansion coefficients $`g_{si}^\alpha (P^2)`$ from the BS amplitude in Eq. (24). Using Eqs. (22) and (23), one finds

$`g_{si}^\alpha (P^2)`$ $`=`$ $`𝒪_{si}(p)t_{\alpha \alpha ^{}}^1(p,P)\mathrm{\Lambda }_\rho ^\alpha ^{}\mathrm{\Gamma }_\rho (p;P)`$ (27)
$`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑p^2{\displaystyle _1^{+1}}𝑑z_p\sqrt{1z_p^2}U_i(z_p)\varphi _s(p^2)`$
$`\times t_{\alpha \alpha ^{}}^1(p,P)\mathrm{tr}_\mathrm{D}\left[\mathrm{\Lambda }_\rho ^\alpha ^{}(p,P)\mathrm{\Gamma }_\rho (p;P)\right].`$
The set of coefficients $`g_{si}^\alpha (P^2)`$ will be refered to as BS amplitudes since they are equivalent to the BS amplitude $`\mathrm{\Gamma }_\rho (p;P)`$ in that they contain all the dynamical information regarding the solution of the BSE. The expansion coefficients $`g_{si}^\alpha (P^2)`$ are functions of $`P^2`$ only. The operator $`𝒪_{si}(p)`$ is an integral operator that is introduced for brevity and notational convenience.

When both sides of the BSE (LABEL:BSE\_ladder) are projected according to Eq. (27) one obtains

$`\lambda (P^2)g_{si}^\alpha `$ $`=`$ $`𝒪_{si}(p){\displaystyle }{\displaystyle \frac{d^4q}{(2\pi )^4}}t_{\alpha \alpha ^{}}^1(p,P)[\mathrm{\Lambda }_\rho ^\alpha ^{}(p)`$ (28)
$`K(p,q;P)\left(S(q_+)\mathrm{\Gamma }_\rho (q;P)S(q_{})\right)],`$
$`=`$ $`{\displaystyle \underset{AB\beta tj}{}}𝒪_{si}(p){\displaystyle }{\displaystyle \frac{d^4q}{(2\pi )^4}}t_{\alpha \alpha ^{}}^1(p,P)[\mathrm{\Lambda }_\rho ^\alpha ^{}(p)`$
$`K(p,q;P)\left(\stackrel{~}{\mathrm{\Lambda }}^A(q_+)\mathrm{\Lambda }_\rho ^\beta (q;P)\stackrel{~}{\mathrm{\Lambda }}^B(q_{})\right)]`$
$`\times \sigma _A(q_+^2)\sigma _B(q_{}^2)\varphi _t(q^2)U_j(z_q)g_{tj}^\beta ,`$
where the sums are over $`\beta \{1,\mathrm{},8\}`$, and $`t,j\{0,1,\mathrm{}\}`$, and $`A,B\{V,S\}`$. The second equality follows by expanding the BS amplitude $`\mathrm{\Gamma }_\rho (q;P)`$ according to Eqs. (22) and (24). In Eq. (28) another set of Dirac covariants, $`\stackrel{~}{\mathrm{\Lambda }}^A(q)=\{i\gamma q,\mathrm{\hspace{0.33em}1}\}`$ with $`A=V,S`$, have been introduced to allow the propagator representation $`S(q)=_{A=V,S}\stackrel{~}{\mathrm{\Lambda }}^A(q)\sigma _A(q^2)`$.

At this point, it is helpful to specify a coordinate system. We take $`P_\mu =`$ $`\sqrt{P^2}(0,0,0,1)`$, $`p_\mu =`$ $`\sqrt{p^2}(0,0,\sqrt{1z_p^2},z_p)`$, and the integration momentum $`q_\mu `$ is represented as

$$q_\mu =\sqrt{q^2}\left(\begin{array}{c}\mathrm{cos}\varphi \sqrt{1x_q^2}\sqrt{1z_q^2}\\ \mathrm{sin}\varphi \sqrt{1x_q^2}\sqrt{1z_q^2}\\ x_q\sqrt{1z_q^2}\\ z_q\end{array}\right).$$
(29)
The integrand in Eq. (28) has no dependence upon angle $`\varphi `$; the only dependence on direction cosine $`x_q`$ is through the $`pq`$ dependence of the BS kernel $`K(p,q;P)`$. In Eq. (28) we may therefore substitute

$`{\displaystyle d^4q}=\pi {\displaystyle _0^{\mathrm{}}}𝑑q^2q^2{\displaystyle _1^1}𝑑z_q\sqrt{1z_q^2}{\displaystyle _1^1}𝑑x_q.`$ (30)
The integration over $`x_q`$ in Eq. (28) produces a function of $`p^2,q^2,z_p,z_q`$ and $`P^2`$ which can be expanded in the basis functions. That is, we define a quantity $`𝒱`$ by

$`{\displaystyle _1^{+1}}dx_qt_{\alpha \alpha ^{}}^1(p,P)[\mathrm{\Lambda }_\rho ^\alpha ^{}(p)K(p,q;P)`$ (31)
$`\left(\stackrel{~}{\mathrm{\Lambda }}^A(q_+)\mathrm{\Lambda }_\rho ^\beta (q;P)\stackrel{~}{\mathrm{\Lambda }}^B(q_{})\right)]`$
$`=`$ $`{\displaystyle \underset{stij}{}}\varphi _s(p^2)U_i(z_p)𝒱_{si;tj}^{\alpha A;\beta B}(P^2)\varphi _t(q^2)U_j(z_q).`$
Upon substituting this back into Eq. (28) above, one obtains a discretized form of the BSE,

$`\lambda (P^2)g_{si}^\alpha (P^2)=𝒱_{si;w\mathrm{}}^{\alpha A;\beta B}(P^2)𝒢_{w\mathrm{},tj}^{AB}(P^2)g_{tj}^\beta (P^2),`$ (32)
where repeated indices are understood as being summed over. Here we have introduced the projection of the product of propagator amplitudes in the form

$`𝒢_{si,tj}^{AB}(P^2)`$ $`=`$ $`{\displaystyle \frac{\pi }{(2\pi )^4}}{\displaystyle _0^{\mathrm{}}}𝑑q^2q^2{\displaystyle _1^1}𝑑z_q\sqrt{1z_q^2}\varphi _s(q^2)U_i(z_q)`$ (33)
$`\times \sigma _A(q_+^2)\sigma _B(q_{}^2)\varphi _t(q^2)U_j(z_q).`$
Thus far, the BSE has been reduced to the discrete eigenvalue problem given in Eq. (32). Physical solutions are at $`\lambda (P^2=M_\rho ^2)=1`$ from which one identifies the existence of a bound state vector meson of mass $`M_\rho `$. The only task remaining is the determination of the elements in Eq. (32) that make up the kernel of the eigenvalue problem. These are $`𝒱`$ and $`𝒢`$. The generalized interaction $`𝒱`$ is determined in terms of the kernel by the projection in Eq. (31).

Consider the explicit form of $`𝒢`$. Clearly there is a strong similarity between the form of Eq. (33) and the simple integral whose analytic continuation was considered in Sec. III. This is clarified by the expansion of the quark propagator amplitudes in terms of the mass pole terms from Eq. (9) in the form

$$\sigma _A(q^2)=\underset{n=1}{\overset{N}{}}\underset{a=\pm 1}{}\frac{Z_n^{Aa}}{q^2+\mu _n^2+ia\mathrm{\Delta }_n^2},$$
(34)
where $`Z_n^{Aa=+1}=z_n`$ or $`z_nm_n`$ for $`A=V`$ or $`S`$, respectively, and $`Z_n^{Aa=1}=z_n^{}`$ or $`z_n^{}m_n^{}`$, respectively. Then, Eq. (33) can be written as

$`𝒢_{st;ij}^{AB}(P^2)={\displaystyle \frac{1}{2}}{\displaystyle \underset{n,m=1}{\overset{N}{}}}{\displaystyle \underset{a,b=\pm 1}{}}Z_n^{Aa}Z_m^{Bb}`$
$`{\displaystyle \frac{d^4q}{(2\pi )^4}\frac{\varphi _s(q^2)U_i(z_q)\varphi _t(q^2)U_j(z_q)}{(q_{}^2+\mu _n^2+ia\mathrm{\Delta }_n^2)(q_+^2+\mu _m^2+ib\mathrm{\Delta }_m^2)}}.`$ (35)
These integrals in Eq. (35) are of the type $`\lambda _{ab}(P^2)`$ already explored in our discussion of analytic continuation in Sec. III. The only difference is that the number of integrals that have to be performed can be large since it depends on the number of basis terms for representation of the BS amplitudes via Eq. (24), and the number of conjugate pairs of mass poles for representation of the propagators.

The same methods employed in Sec. III may be used here to calculate the integrals $`𝒢_{st;ij}^{AB}(P^2)`$ for all values of $`P^2`$ while accounting properly for the analytic continuation. In particular, the Feynman integral method used here will lead to a denominator $`[q^2q_{nmab}^2(P^2,y)]^2`$, where the pole location is

$`q_{nmab}^2(y)`$ $`=`$ $`{\displaystyle \frac{P^2}{4}}(1y^2)+y{\displaystyle \frac{\mu _n^2\mu _m^2}{2}}+{\displaystyle \frac{\mu _n^2+\mu _m^2}{2}}`$ (36)
$`+iy{\displaystyle \frac{a\mathrm{\Delta }_n^2b\mathrm{\Delta }_m^2}{2}}+i{\displaystyle \frac{a\mathrm{\Delta }_n^2+b\mathrm{\Delta }_m^2}{2}}.`$
The integration over $`q^2`$ is carried out first, followed by the integration over the Feynman parameter $`y`$, and finally the integration over $`z_q`$. Integrations over the two remaining angles in Eq. (35) are trivial. In general, for a quark propagator parametrized in terms of $`N`$ pairs of complex-conjugate poles, there are $`(2N)^2`$ pole trajectories, half of which impinge on the $`q^2`$ integration domain and lead to residue additions. The imaginary parts of the integrals in Eq. (35) cancel in pairs as described in the previous section.

In ladder approximation, the determination of $`𝒱`$ is straight-forward. There are never any poles encountered when $`P^2`$ is analytically continued to negative values. Of course, singularities may be encountered when contributions beyond ladder are maintained in the BS kernel, and these must be treated in a way analogous to $`𝒢`$. This possibility is not addressed in the present work.

## V Model Calculations

### V.1 Ladder-Rainbow Kernel

Here we use a particular model for the DSE-BSE ladde-rainbow kernel to numerically implement the preceeding developments. The model is specified by the “effective coupling” $`\alpha _{\mathrm{eff}}(k^2)`$ and we employ the Ansatz Maris and Tandy (1999)

$`{\displaystyle \frac{4\pi \alpha _{\mathrm{eff}}(k^2)}{k^2}}`$ $`=`$ $`{\displaystyle \frac{4\pi ^2Dk^2}{\omega ^6}}\mathrm{e}^{k^2/\omega ^2}`$ (37)
$`+{\displaystyle \frac{4\pi ^2\gamma _m(k^2)}{\frac{1}{2}\mathrm{ln}\left[\tau +\left(1+k^2/\mathrm{\Lambda }_{\mathrm{QCD}}^2\right)^2\right]}},`$
with $`\gamma _m=12/(332N_f)`$ and $`(s)=(1\mathrm{exp}\frac{s}{4m_t^2})/s`$. The ultraviolet behavior matches that of the 1-loop QCD running coupling $`\alpha _s(k^2)`$. The resulting solution of the ladder-rainbow DSE-BSE system of equations in the UV region generates the correct 1-loop perturbative QCD structure. The first term implements the strong infrared enhancement in the region $`0<k^2<1\mathrm{GeV}^2`$ phenomenologically required Hawes et al. (1998) to produce a realistic value for the chiral condensate. With $`m_t=0.5\mathrm{GeV}`$, $`\tau =\mathrm{e}^21`$, $`N_f=4`$, $`\mathrm{\Lambda }_{\mathrm{QCD}}=0.234\mathrm{GeV}`$, and a renormalization scale $`\mu =19\mathrm{GeV}`$, it has been found that $`\omega =0.4\mathrm{GeV}`$ and $`D=0.93\mathrm{GeV}^2`$ give a good description of $`\overline{q}q`$, $`m_{\pi /K}`$ and $`f_\pi `$ with physically acceptable current quark masses, Maris and Roberts (1997); Maris and Tandy (1999). The propagator amplitudes from this rainbow DSE model have recently been shown Jarecke et al. (2002) to have the same qualitative features as lattice QCD simulations Bowman et al. (2002); Bowman (2002).

The vector meson masses and electroweak decay constants produced by this model are in good agreement with experiments Maris and Tandy (1999). Without any readjustment of the parameters, this model agrees remarkably well with the most recent Jlab data Volmer et al. (2001) for the pion charge form factor $`F_\pi (Q^2)`$. Also the kaon charge radii and electromagnetic form factors are well described Maris and Tandy (2000); Maris (2000). The strong decays of the vector mesons into a pair of pseudoscalar mesons are also well-described within this model Maris and Tandy (2001); Jarecke et al. (2002).

### V.2 Results

The DSE solution for the $`u/d`$ quark propagator from the model interaction given in Eq. (37) can be well fit with $`N=3`$ pairs of complex conjugate poles in the representation given in Eq. (8). The result displayed in Fig. 2 corresponds to the parameter set

$`m_1`$ $`=`$ $`0.547+i0.303\mathrm{GeV},`$
$`z_1`$ $`=`$ $`0.200+i0.475,`$
$`m_2`$ $`=`$ $`1.262+i0.570\mathrm{GeV},`$
$`z_2`$ $`=`$ $`0.142+i0.045,`$
$`m_3`$ $`=`$ $`1.560+i0.564\mathrm{GeV},`$
$`z_3`$ $`=`$ $`0.160+i0.015.`$ (38)
This parameterization was chosen so that $`Z(p^2)`$ and $`M(p^2)`$ were reproduced well at both $`p^2=0`$ and $`p^2=100\mathrm{GeV}^2`$ while the main features of the momentum dependence are preserved. We imposed the constraint that $`M(p^2)`$ should approach its UV limit from above. The pole locations are $`p^2=m_i^2`$ where $`m_i^2=`$ $`\mu _i^2+i\mathrm{\Delta }_i^2`$. Thus we have $`\mu _1=`$ $`0.455\mathrm{GeV}`$, $`\mathrm{\Delta }_1=`$ $`0.576\mathrm{GeV}`$, $`\mu _2=`$ $`1.13\mathrm{GeV}`$, $`\mathrm{\Delta }_2=`$ $`1.12\mathrm{GeV}`$, and $`\mu _3=`$ $`1.46\mathrm{GeV}`$, $`\mathrm{\Delta }_3=`$ $`1.33\mathrm{GeV}`$.

With use of this propagator representation in the BSE, a $`u/d`$ quark meson bound state with mass greater than the lowest pseudo-threshold $`2\mu _1=`$0.910 GeV would be needed to test the method we have described. The model-exact $`\rho `$ meson mass from the kernel Eq. (37) is 0.742 GeV Maris and Tandy (1999) and is therefore not suitable. Heavier $`u/d`$ mesons have not been carefully studied within this model interaction due to the very singularity issue we are concerned with. In fact the evidence from exploratory studies is that, for example, $`m_{a_1}`$ in this ladder-rainbow model is only about 200 MeV above $`m_\rho `$ and also does not provide a clear test.

To test the present method, we artificially modify the pole parameterization of the quark propagator so that the lowest pseudo-threshold moves below the vector meson mass, and the latter moves to 0.770 GeV. We achieve this with the following parameterization in terms of $`N=2`$ pairs of conjugate mass poles

$`m_1`$ $`=`$ $`0.400+i0.200\mathrm{GeV},`$
$`z_1`$ $`=`$ $`0.15,`$
$`m_2`$ $`=`$ $`0.550+i0.350\mathrm{GeV},`$
$`z_2`$ $`=`$ $`0.35+i0.37.`$ (39)
The resulting propagator amplitudes are also shown in Fig. 2. This parameterization corresponds to $`\mu _1=`$ $`0.346\mathrm{GeV}`$, $`\mathrm{\Delta }_1=`$ $`0.4\mathrm{GeV}`$, and $`\mu _2=`$ $`0.424\mathrm{GeV}`$, $`\mathrm{\Delta }_2=`$ $`0.621\mathrm{GeV}`$. With this propagator, the lowest pseudo-threshold in the BSE is at $`2\mu _1=`$0.692 GeV. The subsequent BSE solution is not to be taken as a physical represenation of the $`\rho `$ meson; its purpose is to provide a test calculation where the propagator singularities are within the domain of integation.

The basis functions used to expand the $`q^2`$ dependence of the BS amplitudes are

$$\varphi _s(q^2)=\frac{N_s}{\mathrm{\Lambda }}L_s(w)e^{\frac{w^2}{2}}$$
(40)
where $`w=q^2/\mathrm{\Lambda }^2`$, the typical hadron scale is $`\mathrm{\Lambda }=`$ 0.800 GeV, $`L_s(w)`$ are Laguere polynomials, and $`N_s`$ is fixed by the normalization condition given in Eq. (26). To obtain better than 1% accuracy for the eigenvalue $`\lambda (P^2)`$, we found that 7 terms in this basis were sufficient. The number of Chebyshev terms used for the angle basis was one. The first covariant from the set given in Eq. (21) is known to be dominant for the vector meson solution Maris and Tandy (1999) and this was the only one used here.

The resulting BSE eigenvalue $`\lambda (P^2)`$ is given in Fig. 3 as the solid curve. In the spacelike region ($`P^2>0`$) the eigenvalue is a monotonically, decreasing function of $`P^2`$. When $`P^2`$ decreases below zero and becomes timelike, the eigenvalue $`\lambda (P^2)`$ increases for a while then crosses unity, at $`P^2=M_\rho ^2`$, signifiying a bound state at $`M_\rho =0.770`$ GeV. In Fig. 3 this value of $`P^2`$ is shown as a thin, vertical line. On the timelike side of this line, the eigenvalue continues to increase, then begins to decrease and ultimately goes to zero for ($`P^20`$). This is because the eigenvalue is a measure of the magnitude of the BSE kernel and the propagators therein falloff as the difference between the momentum arguments and the pole positions becomes larger than any mass scale in the interaction.

The solid curve in Fig. 3 represents the proper analytic continuation of the BSE from $`P^2>0`$ to $`P^2<0`$. The dashed line represents the result from a direct evaluation of the BSE integral without regard to the possibility of singularities and their movement as a function of $`P^2`$ in relation to the integration domain. That is, only a transcription from $`P^2>0`$ to $`P^2<0`$ has been made in the integrand of Eq. (35); the Feynman technique for combining the denominators and shifting the 4-momentum variable to complete the square has not been made. The integration variables are $`q^2`$ and $`z_q=\widehat{q}\widehat{P}`$. It is clear that for the spacelike domain ($`P^2>0`$), and for the limited timelike domain for which no poles move into the domain of integration ($`4\mu _1^2<P^2<0`$, i.e., below the lowest pseudo-threshold), the two approaches give identical results. The dashed line represents the approach to the BSE that has usually been implemented; it has been limited to states light enough to avoid the propagator singularities.

With the parameters given above, the quark poles are first encountered within the BSE integration when the analytic continuation reaches $`P^2=4\mu _1^2`$, that is, at $`P^2=0.48\mathrm{GeV}^2`$. Although the contibutions of the conjugate pair of poles to the imaginary part of the BS eigenvalue $`\lambda (P^2)`$ cancel, their real contribution produces a discontinuity in the derivative $`\lambda ^{}(P^2)`$ at that point. Careful inspection of the solid curve in Fig. 3 reveals that this cusp occurs at the precise value where the full calculation and the naïve calculation diverge.

## VI Discussion

We have provided a method to obtain meson bound state solutions from the Bethe-Salpeter equation in Euclidean metric when the dressed quark propagators have time-like complex conjugate mass poles within the integration domain. This approximates features encountered in recent QCD modeling via the Dyson-Schwinger equations; the absence of real mass poles simulates quark confinement. For this exploratory study we represent the quark propagators as a sum of complex-conjugate mass poles, and project the BSE on to complete basis sets to represent the momentum and angle dependence. We use Feynman integral techniques to combine the propagator denominators and map the integration to a one-dimensional domain. This allows a clear analysis of the analytic continuation in total meson momentum needed to reach the mass shell and which causes the singularities to impinge upon the integration domain. The BSE linear eigenvalue remains real; in other words the eigenmass remains real. The would-be decay width from one pole is exactly cancelled by the effect of the partner pole; the meson is stable against decay into a quark-antiquark pair. This describes the confinement of quark and antiquark within the meson bound state, even though the meson mass is “above threshold”.

One of the limitations of the method presented here is that it relies upon the projection of the BSE on to a complete basis set of functions to represent the momentum dependence of the BS amplitudes. In general the behavior of the BS amplitudes is not known until after the BSE is solved. If the chosen basis set is an inefficient representation, one would expect this to result in a lack of convergence. The (vector) meson mass that we have studied here is an integrated quantity dominated by infrared physics; it is not surprizing that convergence is easily achieved with the basis states of Eq. (40) which fall off exponentially in the ultraviolet region. However, we know from studies such as Ref. Maris and Tandy (1999) that the leading ultraviolet power law behavior of the vector meson BS amplitudes is $`1/q^2`$. Clearly the present basis would be inefficient for applications that depend strongly on the BS amplitude in this domain, such as the asymptotic behavior of the pion charge form factor Maris and Roberts (1998). For such studies, appropriate basis sets would have to be utilized.

The principal reason that the BS amplitudes are expanded in terms of known basis functions is that the Feynman integral technique for handling the two propagator denominators requires a complex variable shift for the rest of the integrand. In principle, the analytic continuation of the BS amplitude into the complex plane is determined by the dynamics of the BSE and can be determined only through the solution. The implementation of the present method of solution requires a known analytic behavior for the BS amplitudes as expressed through the basis functions. One would expect an inadequacy of the basis in this respect to show up as poor convergence.

Future work using the BSE solution method presented here will include extension of the DSE modeling approach to SU(3) flavor meson states above 1 GeV and meson form factors for $`Q^2>3\mathrm{GeV}^2`$.

###### Acknowledgements.

This work is supported by the National Science Foundation under grant Nos. PHY-0071361 and INT-0129236.


format: text

text:
# Temperature Dependent Magnetic Anisotropy in (Ga,Ma)As Layers

## Abstract

It is demonstrated by SQUID magnetization measurements that (Ga,Mn)As films can exhibit rich characteristics of magnetic anisotropy depending not only to the epitaxial strain but being strongly influenced by the hole and Mn concentration, and temperature. This behavior reflects the spin anisotropy of the valence subbands and corroborates predictions of the mean field Zener model of the carrier mediated ferromagnetism in III-V diluted magnetic semiconductors with Mn. At the same time the existence of in-plane uniaxial anisotropy with  the easy axis is evidenced. This is related to the top/bottom symmetry breaking, resulting in the lowering of point symmetry of (Ga,Mn)As to the C2v symmetry group. The latter mechanism coexists with the hole-induced cubic anisotropy, but takes over close to $`T_\text{C}`$.

The discovery of carrier-mediated ferromagnetism at temperatures in excess of 100 K in (III,Mn)V dilute magnetic semiconductors (DMS) grown by molecular beam epitaxy (MBE) has made it possible to combine complementary properties of semiconductor quantum structures and ferromagnetic systems in single devices, paving the way for the development of functional semiconductor spintronics.Spintronics  Therefore, the understanding of magnetic anisotropy in these systems and the demonstration of methods for its control is timely and important. It has been known, since the pioneering works of Munekata et al.Mune93  and Ohno et al.,Ohno96Shen97  that ferromagnetic (In,Mn)As and (Ga,Mn)As films are characterized by a substantial magnetic anisotropy. Remarkably, due to magnetic dilution, the ordinary shape anisotropy plays here only a marginal role and, accordingly, explains neither direction nor large magnitude of the observed anisotropy field $`H_{un}`$.

It has been found by anomalous Hall effect studies,Ohno96Shen97  that the direction of the easy axis is mainly controlled by epitaxial strain in these systems. Generally, for layers under tensile biaxial strain \[like (Ga,Mn)As on a (In,Ga)As buffer\] perpendicular-to-plane magnetic easy axis has been observed (perpendicular magnetic anisotropy, PMA). In contrast, layers under compressive biaxial strain \[as canonical (Ga,Mn)As on a GaAs substrate\] have been found to develop in-plane magnetic easy axis (in-plane magnetic anisotropy, IMA). At first glance, this sensitivity to strain appears surprising, as the Mn ions are in the orbital singlet state $`{}_{}{}^{6}A_{1}^{}`$. For such a state, the strain-induced single ion anisotropy is expected to be rather small and, indeed, electron paramagnetic resonance (EPR) studies of Mn in GaAs yielded relevant spin Hamiltonian parameters by two orders of magnitude too small to explain the observed values of $`H_{un}`$.Fedo02

In the system in question, however, the ferromagnetic spin-spin exchange interaction is mediated by the band holes, whose Kohn-Luttinger amplitudes are primarily built up of anion p orbitals in tetrahedrally coordinated semiconductors. Furthermore, in semiconductors, in contrast to metals, the Fermi energy is usually smaller than the atomic spin-orbit energy. Hence, as noted by some of the present authors and co-workers,Diet97 ; Diet00  the strain or confinement-induced anisotropy of the valence band can result in a sizable anisotropy of spin properties. Indeed, the quantitative calculation within the mean-field Zener model,Diet00 ; Diet01a  in which the valence band is represented by the 6$`\times `$6 Luttinger hamiltonian, explains the experimental values of $`H_{un}`$ in (Ga,Mn)As with the accuracy better than a factor of two. Moreover, by combining theories of magnetic anisotropyDiet00 ; Diet01a ; Abol01  and of magnetic stiffness,Koni01  it has been possible to describe the width of stripe domains in (Ga,Mn)As PMA films.Diet01b  However, the theories in questionDiet01a ; Abol01  contain a number of predictions that call for a detail experimental verification.

In this communication we present magnetic anisotropy studies carried out by direct magnetization measurements in a dedicated SQUID magnetometer. Our results demonstrate that the temperature-induced reorientation of the easy axis from  to , and then to  or equivalent directions occurs in films of (001) (Ga,Mn)As on GaAs with appropriately low values of $`p`$. We argue that these findings corroborate the above mentioned theoretical expectations. However, the data collected near $`T_\text{C}`$ reveal a non-equivalence of the $`[110]`$ and $`[\overline{1}10]`$ crystal directions, the latter corresponding to the hard axis. This points to the existence of a mechanism lowering the point symmetry from D2d to C2v. We find that this particular anisotropy develops independently of $`p`$. Such a symmetry breaking has already been inferred from transport data Kats98 ; Gall02  and Kerr rotation.Hrab02  Finally, we find that in another class of samples the hard axis remains oriented along the $`[\overline{1}10]`$ direction but the three other main directions (,  and ) become equally easy. We discuss reasons and consequences of these unusual characteristics of magnetic anisotropy in (Ga,Mn)As.

Two series of (Ga,Mn)As films with the thickness of 350 and 500 nm were coherently deposited by MBE onto (100) GaAs substrates under As4/Ga beam equivalent pressure (BEP) ratio of 5 and 25, and growth temperature $`T_{sub}`$ = 220 or 270$`^\text{o}`$C. The Mn composition $`x`$ varies from 0.007 to 0.067, as established by x-ray diffraction measurements.Scho01  High resolution x-ray diffraction shows high crystal quality of (Ga,Mn)As with the rocking curve widths comparable to those of the GaAs substrate and pronounced finite-thickness fringes, indicating flat interfaces and surfaces. The layers are pseudomorphic with respect to the GaAs substrate, as verified by reciprocal space maps around the asymmetric  reflex.

Magnetization measurements are carried out in a home made superconducting quantum interference device (SQUID) magnetometer down to 5 K and in the magnetic fields up to 3 kOe. By design, the signal detection axis of the magnetometer (denoted in Fig. 1 as am), is aligned with the direction of the magnetic field *H* (vertically). As a result the system is sensitive only to the vertical component of the magnetization vector. A special care is put forward to screen the sample from external magnetic fields and to keep the parasite residual fields generated by the magnet at the lowest possible level. The Meissner effect of a pure lead sample proves that even after excursion to the field of 1 kOe, the residual field remains below 60 mOe. Such a low value is essential for successful magnetic remanence studies or for measurements at temperatures close to $`T_\text{C}`$, where the coercivity of (Ga,Mn)As drops into a sub-Oersted range. For studies of in-plane anisotropy, (Ga,Mn)As layers are shaped by chemical etching into circles of about 5 mm in diameter.

Figures 1(a) and 1(b) present typical hysteresis loops for a Ga0.977Mn0.023As sample, which exhibits the temperature driven change from PMA to IMA. At low temperatures, Fig. 1(a), a perfect square hysteresis is obtained when the magnetic field is perpendicular to the film surface, while an elongated loop is seen when the in-plane orientation is probed. This clearly evidences that the easy axis is perpendicular. Remarkably, a reverse behavior is observed at higher temperatures, Fig. 1(b). These findings demonstrate that the easy axis flips from the perpendicular to in-plane orientation when the temperature increased.

In order to trace directly the reorientation of spontaneous magnetization, we have examined the temperature dependence of remanent magnetization $`M_{\text{REM}}`$ along selected crystallographic directions according to the following procedure. After being oriented on the holder, the sample is cooled down through $`T_\text{C}`$ in the field *H*$`{}_{\text{FC}}{}^{}=1`$ kOe, which is at least a few times larger than the coercive field $`H_c`$. Then, the field is removed at 5 K, and the measurement of the $`M_{\text{REM}}`$ component along the direction of *H*$`_{\text{FC}}`$ commences on increasing temperature in the residual field *H*$`{}_{r}{}^{}60`$ mOe.

As shown in Fig. 1(c), the $`M_{\text{REM}}`$ signal at low temperatures is tenfold larger for the perpendicular orientation, *H*$`{}_{\text{FC}}{}^{}[001]`$, comparing to the case of parallel arrangement, *H*$`{}_{\text{FC}}{}^{}[100]`$. This reconfirms the appearance of PMA for large values of spontaneous magnetization $`M_s`$ in this film despite the presence of a sizable compressive strain. The PMA holds until a certain temperature $`T_\text{L}`$ ($`T_\text{L}16`$ K in this case), at which a sudden drop of the signal is detected in the perpendicular orientation. Nevertheless, the $`M_{\text{REM}}`$ signal for the parallel case remains small. This seems to contradict our scenario, according to which the easy axis is in-plane for $`T>T_\text{L}`$. We note, however, that the spontaneous flip of magnetization to the in-plane configuration results in a demagnetized state. Actually, a visible maximum of the parallel $`M_{\text{REM}}`$ between $`T_\text{L}`$ and $`T_\text{C}`$ reflects the presence of the field *H*r, which magnetizes the film once the in-plane anisotropy and coersive fields become sufficiently small. Indeed, a temperature cycling of about $`\mathrm{\Delta }T=2`$ K just below this maximum increases the magnitude of the in-plane magnetization by a factor of three, as shown in Fig. 2.

We have found the temperature-induced cross over from PMA to IMA described above in five other samples of (Ga,Mn)As.Sawi02  Recently, this effect has been detected in (Al,Ga,Mn)As under similar conditions.Taka02  According to low temperature resistance measurements, such samples exhibit an insulating character, which points to rather small hole concentrations $`p`$, a conclusion consistent with relatively low magnitudes of the Curie temperature in these films, $`T_\text{C}<30`$ K. Remarkably, the opposite behavior occurs in (In,Mn)As under a tensile strain, where the easy axis flips from the in-plane to the out of plane on warming.Endo01  We assign the temperature dependent anisotropy to the repopulation of the valence band subbands, which occurs when the Mn magnetization changes. For a qualitative analysis we refer to numerical predictions depicted in Fig. 10 of Ref. Diet01a, . It is expected that for compressive (tensile) strain the easy axis changes direction from the in-plane to perpendicular orientation on lowering (raising) $`p`$, the critical value being of the order of $`p_c=5\times 10^{19}`$ cm-3 for $`M_s`$ corresponding to the saturation value of magnetization for $`x=0.02`$. We conclude that the theoretically expected value of $`p_c`$ is consistent with our results for (Ga,Mn)As/GaAs, which demonstrate the existence of PMA only in samples with sufficiently small values of $`p`$. Furthermore, at given $`p`$, a cross-over from the PMA to IMA is expected on decreasing $`M_s`$, that is on increasing $`T`$, as observed. This temperature-driven reorientation of the easy axis reflects the equipartition of the valence subband population at $`T`$ approaching $`T_\text{C}`$ and, therefore, reconfirms the crucial role of the valence band holes in the ferromagnetism of (III,Mn)V systems. Unfortunately, however, a more quantitative comparison is hampered by the well-known difficulties in the accurate determination of $`p`$ in films of (Ga,Mn)As, particularly near the metal-insulator transition.

We now turn to (Ga,Mn)As films which, owing to sufficiently large $`p`$, exhibit IMA in the whole temperature range above 5 K. In this case, according to the theoretical predictions presented in Fig. 9 of Ref. Diet01a,  and in Fig. 6 of Ref. Abol01, , the easy axis is expected to switch between $`[100]`$ and $`[110]`$ or equivalent in-plane directions as a function of $`p`$. Alternatively, at given $`p`$, the easy direction may change when $`M_s`$ is varied, say, by $`x`$ or $`T`$. In fact, there exists a class of samples, in which we observe the temperature-induced change of the IMA. Data collected in Fig. 3, by presenting $`M_{\text{REM}}`$ measured along the four principal in-plane crystallographic directions $`[100]`$, $`[110]`$, $`[010]`$, and $`[\overline{1}10]`$, exemplify this case. We note that at low temperatures the highest values of $`M_{\text{REM}}`$ are collected along $`100`$ directions. However, this alters at about 10 K, above which the  becomes the most preferred direction in the system. It is not only that $`M_{\text{REM}}`$ is the stronger for this direction but also this remanence is by a factor of $`\sqrt{2}`$ larger than that for the $`100`$ directions, as is expected from simple geometric considerations.

Surprisingly, however, the data reveals the existence of an uniaxial anisotropy along $`[110]`$. As shown in Fig. 3, $`M_{\text{REM}}`$ measured along $`[\overline{1}10]`$ direction vanishes completely above 15 K indicating that this is the hard direction of the system. At the same time, the magnitude of $`M_{\text{REM}}`$ is the same for both $`100`$ directions. We argue that this is the result of the lack of the inversion symmetry in (Ga,Mn)As (001) films. Due to the biaxial strain, the initial Td point symmetry of zinc-blende structure is lowered to D2d. But the bottom (Ga,Mn)As/GaAs interface is different than the top (Ga,Mn)As/vacuum one. This further lowers the symmetry to C2v, where the three principal directions are: ,  and $`[\overline{1}10]`$. They are not equivalent, in C2v the $`[110]`$ $`[\overline{1}10]`$ symmetry gets broken, while the $`[100][010]`$ one is maintained through reflectivity operations by the $`(110)`$ and $`(\overline{1}10)`$ planes, which conforms with the presented results. Since the cubic-like anisotropy energy is proportional to $`M_s^4`$ whereas the uniaxial one to $`M_s^2`$, the latter is dominating at high temperatures, where $`M_s`$ is small. We also note that when $`M_{\text{REM}}^{[\overline{1}10]}`$ vanishes the $`M_{\text{REM}}^{100}/M_{\text{REM}}^{[110]}`$ ratio drops to $`\sqrt{2}/2`$, the value expected for the easy axis along $`[110]`$.

Another class of samples are those, in which relatively high values of $`p`$ and/or $`x`$ make the cubic easy axis to be oriented along $`100`$ in the whole relevant temperature range. In such samples, the coexistence of the cubic and uniaxial anisotropy leads to the case ”not $`[\overline{1}10]`$”, as presented in Fig. 4(a). This is generally the $`[100]`$ easy axis samples in which the  direction gets equally easy as the $`[100]`$ directions at high temperatures. So, close to $`T_\text{C}`$ these three major directions are virtually equivalent, while the $`[\overline{1}10]`$ axis stays hard. As shown in Fig. 4(b), such an anisotropy changes considerably the shape of the hysteresis loop, yielding $`M^{[\overline{1}10]}(H=0)=\sqrt{2}M_r/4`$, instead of zero, as observed in the previous case, for which both cubic and uniaxial easy axes were along the $`[110]`$ direction.

In summary, our studies have demonstrated the rich characteristics of magnetic anisotropies in (Ga,Mn)As/GaAs, which – in addition to epitaxial strain – vary with the hole and Mn concentrations as well as with the temperature. According to theory,Diet01a ; Abol01  these reflect spin anisotropy of the valence band subbands whose shape depends on strain, while their splitting and population on magnetization and hole concentration. At the same time, our findings have provided the magnetic corroboration for the existence of uniaxial in-plane anisotropy. By group theoretical considerations we have linked this unexpected anisotropy to the top/bottom symmetry breaking, a suggestion calling for a microscopic modelling. Furthermore, in view of our results, a complex temperature dependence of spontaneous magnetization and of related properties, suggesting — for instance — the presence of phases with different Curie temperatures, can actually stem from the temperature-induced switching of the easy axis direction. Since this orientation depends on the hole concentration which, in turn, can be varied by the electric field,Ohno00  it appears possible to switch magnetization direction in a field-effect transistor with the channel of a ferromagnetic semiconductor.

###### Acknowledgements.

The authors thank J. Ferré, H. Ohno and W. Van Roy for valuable discussions. Support of the FENIKS project (EC:G5RD-CT-2001-0535) is gratefully acknowledged.


text:
# 1 Introduction

## 1 Introduction

Because a background of relic particles gives a snapshot of the state of the universe at the time when these particles decoupled from the primordial plasma, relic gravitational waves (gravitons) are a potential source of informations on the state of the very early universe and the physics at correspondingly high energies, which can not be accessed experimentally in any other way?. The photons of the cosmic microwave background (CMB) give us a snapshot of the state of the universe at $`t3\times 10^5`$ years after the big-bang, while the gravitons of the stochastic gravitational wave background (SGWB) encode in its spectrum the information of the universe at $`t10^{44}`$ seconds. These properties make the detection of SGWB very interesting, which is the main purpose of the proposal discussed here.

A SGWB can be characterized by its energy distribution in frequency, especially the dimensionless quantity

$$\mathrm{\Omega }_{\mathrm{gw}}(f):=\frac{1}{\rho _{\mathrm{critical}}}\frac{d\rho _{\mathrm{gw}}}{d\mathrm{ln}f},$$
(1)
where $`\rho _{\mathrm{critical}}=3c^2H_0^2/8\pi G`$ is the critical energy density required (today) to close the universe (the Hubble constant $`H_0=100h_{100}`$ km sec-1 Mpc-1 = $`3.2\times 10^{18}h_{100}`$ Hz ). Up to now, there have been three observational constraints on SGWB spectrum?,?: (i) The strongest observational constraint on $`\mathrm{\Omega }_{\mathrm{gw}}(f)`$ comes from the high degree of isotropy observed in the CMB, $`\mathrm{\Omega }_{\mathrm{gw}}(f)h_{100}^2<7\times 10^{11}\left(H_0/f\right)^2\mathrm{for}H_0<f<30H_0.`$ (ii) The second comes from the standard model of big-bang nucleosynthesis, $`_{f>10^8\mathrm{Hz}}d\mathrm{ln}f\mathrm{\Omega }_{\mathrm{gw}}(f)h_{100}^2<10^5.`$ (iii) The last observational constraint comes from monitoring the radio pulses emitted by a number of stable millisecond pulsars for a decade, $`\mathrm{\Omega }_{\mathrm{gw}}(f=10^8\mathrm{Hz})h_{100}^2<10^8.`$

In order to get a lower constraint on SGWB spectrum or even directly detect SGWB within the very low frequency band (VLF), $`10^7`$ to $`10^9`$ Hz, we propose to construct a special 50 m radio telescope at National Astronomical Observatories, Chinese Academy of Sciences (NAOC) to monitor about 10 strong millisecond pulsars. After reviewing the basics of pulsars as probes for gravitational wave (section 2), we present our research plan in details (section 3).

## 2 Pulsar timing array as efficient gravitational wave detectors

Pulsars have been believed to be extremely stable natural clocks. Detweiler? had shown that measurements of pulse arrival time (TOA) may be used to search for gravitational waves with a period on the order of years. For SGWB, his found

$$\rho _{\mathrm{gw}}=\frac{243\pi ^3f^4R^2(t)}{208G}.$$
(2)
where the residual $`R(t)`$ is the difference between observed and expected TOA. The span of measurements T should be longer than $`1/f`$ (with $`1/f=`$ a few years, the measurements may be taken at intervals of say, 10-20 days).

The observed $`R(t)_{\mathrm{obs}}`$ consists of deviations caused by various effects. Some of them are systematic, such as the uncertainties in the position and proper motion of pulsars, errors in the positions of barycenter. While others are random, for example, propagation delays caused by interstellar turbulence, relativistic effects by intervening massive objects, irregularities occurred in pulsar parameters, or in reference clock. Some of these effects could be modeled or be estimated, and hence be corrected. The resulted $`R^2(t)`$ is therefore the mixture of three components: (1) the remaining effects left over after the correction; (2) pulse timing noise which is related with (a) signal-to-noise ratio of pulsar observation, (b) pulse shape and the definition of arrival time, and (c) method in obtaining the arrival time; (3) background gravitational waves. Thus the quantity $`R^2(t)`$ from pulsar timing observations produces an upper limit on the energy density of SGWB.

Following Detweiler’s work, a number of observations had been made over a decade due to the advent of the discovery of millisecond pulsars. Millisecond pulsars have high stability and sharp pulses, therefore the precision of timing measurements is extraordinary high, reaching $`0.1\mu \mathrm{s}`$ sometimes. A millisecond pulsar, besides its extraordinary high inertia which warrants high stability, the shape of its pulses being so sharp has caused a great increase in the precision of timing measurements. As a result, the energy density of SGWB in the frequency range around 0.5 cycle/year was found to be $`\rho <10^6\rho _{\mathrm{critical}}`$ (see eg. Ref 5).

This result is encouraging. In particular, a much better sensitivity will be obtained if one can use two or more pulsars as coincident detectors?,?,?,?. Such a “pulsar timing array” could be a very efficient SGWB detector. To quantify this, consider the Doppler shift from the $`i`$th pulsar in the array?

$$\frac{\mathrm{\Delta }\nu _i(t)}{\nu _i}=\alpha _is(t)+n_i(t)$$
(3)
where $`s(t)`$ is the GW signal, $`\alpha _i`$ a geometrical factor depending on the line-of-sight direction, and $`n_i(t)`$ contains the noise intrinsic to the pulsar timing. The cross-correlation between two pulsars then gives?

$$\alpha _i\alpha _js^2(t)+\alpha _isn_i+\alpha _jsn_j+n_in_j$$
(4)
and, since the timing noises from different pulsars are uncorrelated, increasing the observation time will limit the cross-correlation towards $`\alpha _i\alpha _js^2(t)`$.

Unfortunately, the galactic millisecond pulsars known around 1990, for which timing data have been taken already, were concentrated in almost the same region of the sky, and correlation analysis of these data does not improve much the result. Now the known millisecond pulsars are much more uniformly distributed on the sky. It is the right time to construct such a “pulsar timing array”?. The Berkeley pulsar group led by Don Backer has been working on this direction for a few years.

## 3 A special 50 m radio telescope at NAOC

We propose to construct a special 50m radio telescope at NAOC, devoted to pulsar timing measurement. It will work on low frequencies, with 1400 MHz as the highest (so that the instruments will not be very expensive). At least 10 strong millisecond pulsars can be chosen from the 27 candidates listed in Table 1 for regular timing measurements. Each object should be observed 1 hour every day. With the best receiving system today, accumulation of 15 day’s measurements will have a sensitivity of about 0.05 mJy, so that a signal-to-noise ratio above 10 is expected for these pulsars.

In addition to what mentioned above, multi-object observation strategy is of importance for the detection of SGWB in the following senses: Consider a number of pulsars separated at different parts in the sky, gravitational wave from any direction sweeping across the Earth will be marked by the change of the rates of local clock at a frequency $`f`$. Owing to the quardrupole character of the gravitational radiation, such effect will exactly in two pulsar data if their directions are $`180^{}`$ apart. So that $`R(t)`$, in the form of random fluctuations, contained in the record of timing measurements of a pulsar would have an exactly similar shape with the $`R(t)`$ record of a pulsar. As the period of the gravitational wave $`1/f`$ being very long, observations of two pulsars within a time interval of one day can be thought as nearly simultaneous. $`R(t)_{\mathrm{obs}}`$ of the two records, when cross-correlated, will reveal the existence quantitatively of the background gravitational waves. In practice, we will have a number of pairs of pulsars at different separation angles, by the same principle, cross-correlations performed to the $`R(t)_{\mathrm{obs}}`$ of all the pairs will produce a better result than that from a single pair.

Finally, we would like to remark that, among these sources of perturbations contained in $`R(t)`$, many of them are of interest in their own right. The timing measurement data of these objects may be used to, for examples, (1)improve the ephemeris, (2)monitor the intrinsic perturbations of pulsars, (3)establish a time standard for long time scales.

Acknowledgements

The authors are grateful to Prof. JinLin Han for providing us the data in Table 1 and careful reading and improving the manuscript, and also to Prof. Wei-Tou Ni for helpful comments and suggestions. Z.H.Zhu is supported by the Special Funds for Major State Basic Research Projects, the National Natural Science Foundation of China (19903002) and the JSPS fellowship of Japan.

References