text:
# A current circuit model of pulsar radio emission

## 1 Introduction

Pulsars are generally accepted to be rotating neutron stars having $`10^410^8`$T magnetic fields that induce electric potential drops along the magnetic field lines which allow efficient particle acceleration of thermally emitted electrons in the region of open field lines (Goldreich and Julian (1969); Sturrock (1971); Ruderman and Sutherland (1975), Arons (1981)). Although the emitted radio luminosity in comparison to the spin-down power of a rotating neutron star is tiny (only the $`10^5`$\- $`10^6`$ part of the spin-down power is emitted in the radio regime) one of the most striking features about pulsars is their very intense pulsed radio emission, the physical origin of which is still poorly understood. The wealth of observational data (intensity, polarization, shape of the radio spectrum) clearly indicates that the pulsar radio emission cannot be of incoherent origin, but must be due to some coherent radiation mechanism (e.g. Lyne and Graham-Smith 1998; Kramer, Wex and Wielebinski (2000) and references therein). Especially at low frequencies (below 1 GHz) a very high degree of coherence is needed to explain the observed intensities in terms of brightness temperatures. Therefore a pulsar emission model has to explain how a large number of phase coupled particles can be stimulated to radiate in a coherent manner.

Currently the most favored models for radio emission mechanism are curvature radiation (e.g. Rankin (1992); Radhakrishnan and Rankin (1990); Gil (1992)) or maser process due to strong turbulence, like the free electron maser, (hereafter FEM) (e.g. Melrose (1978), Rowe (1992a,b)).

Since there are severe difficulties with coherent curvature radiation (Lesch et al. (1998) pointed out that the energetic requirements can hardly be matched, Melrose (1992) discussed the absence of an appropriate bunching mechanism to build up a coherent curvature radiation) we propose a process directly coupled to the most significant eigenfrequency of the magnetospheric plasma: the electron plasma frequency. It has been shown by Melrose (1978) that such an amplified linear acceleration emission (ALAE) process could produce sufficient growth rates i.e. is fast enough to explain the observations. But it needs low Lorentz factors of the emitting particles to work. This model requires low energy particles with energies of a few MeV.

Recently, support for such a low energy scenario came from various criticisms on the so-called standard model with an inner gap where a highly relativistic pair plasma is created:

It was argued by Kunzl et al. (1998a) and confirmed by Melrose and Gedalin (1999, 2000) that the observed radio emission would not escape from a magnetosphere filled with a dense relativistic pair plasma. Additional problems with the production of high energy particles close to the pulsar were discussed by Arendt & Eilek (2000), Lesch et al. (2000) and Jessner et al. (2001). They commonly concluded that the assumption of strong electric fields at the pulsar surface that ultimately result in ultrarelativistic particle energies is at least questionable.

Thus, we introduce an entirely different approach here, which has first been promoted by Shibata (1991) for an analysis of the global magnetospheric structure of rotation-powered neutron stars. In contrast to the standard model assuming a shielding of the electric fields by massive pair production we propose that in the inner region of the magnetosphere (i.e. the open field lines inward of the neutral surface) the situation is similar to the multiple scatterings in a solid state conductor. In other words we have a current whose time average is the same at all cross sections through the open field line region. There is neither acceleration close to the surface (as indicated by our recent findings (Jessner et al. 2001) nor deceleration in the dissipation region (as indicated by the very small dissipation rate given by small the radio luminosity in comparison to the total spin down power).

In this conducting wire model any dissipation region corresponds to a resistance where electric potentials are present. The ”deceleration” caused by the losses is balanced by the ”acceleration” due to the electric field. The most important consequence of this scenario is that the energy loss per particle is no longer limited by the average kinetic energy of the electron in the same way as the kinetic energy of the electrons in a copper wire current circuit is no indicator for the dissipated energy. The one-dimensional approach with a constant voltage applies, because the resistivity perpendicular to the magnetic field is small inside the pulsar (to allow sufficient supply of charged particles for a Goldreich-Julian DC current) and large enough in the magnetosphere (to avoid ”leakage”, meaning current closure across the field lines in the inner magnetosphere) (see Ewart et al. 1975; Itoh 1975).

Our aim is now to obtain the constraints on the dissipation field in such a model and see whether we can explain the luminosities in the radio frequency range at the emission heights given by the observations from a thin radiating layer.

## 2 The model

In a current-carrying collisionless plasma waves are known to introduce resistivity via fluctuating electric fields. These fields present the deviations from the ideal plasma state with an infinitely high electrical conductivity. For an incoming current $`I`$ the locations of a significant high plasma wave intensity act as a resistance $`R`$. The emitted power of such a resistive region is given by $`L=I^2R`$. In a time averaged picture the current has to be divergence-free. So if a DC current flows through a region in which plasma waves generate a resistance the particles lose their kinetic energy meaning that their velocities decline. Although for relativistic particles this effect is small it cannot be entirely neglected.

In a dense pair plasma this energy loss of the particles presents no difficulties since the velocities of the two components can adjust in an adequate way to fulfill both energy and current conservation. However, in a single charge plasma deceleration has to be accompanied by a density rise. So if the particle beam would gradually lose energy the relative density would have to increase all the time. This is totally incompatible with the outer boundary condition of a Goldreich-Julian magnetosphere.

The latter point is the main difference to the argumentation by Shibata (1997) who also assumes a current outflow determined by the system. In Shibata’s model curvature effects causes space charges to accumulate that cannot be removed by acceleration on some field lines as the particle velocity is limited by $`c`$. We assume that these curvature effects do not play a significant role as they will be balanced by particles of the opposite charge drifting inwards from the neutral line and therefore only deceleration by losses will result in a space charge. This would mean that radiation losses produce a global deviation from the GJ case if there was no acceleration of the particles. So such a description would immediately turn out to be inconsistent (non-idealities induce electric fields which accelerates or decelerates particles).

Nevertheless there is a consistent description of a single charge plasma distribution if we recall that the rotating neutron star induces the vacuum electric field. The electric charges on the neutron star surface are displaced in such a way that the field that has to be shielded by the GJ charge and the field is actually anchored in the neutron star. Thus, the vacuum field can be interpreted as a stationary distribution of infinitely heavy positive particles whose strength is given by the entire pulsar system.

This means that electrons that suffer losses will immediately be dragged on by the electric field that builds up in such a way that just an electrostatic wave is produced. A particle transversing this wave is accelerated by this non-vanishing potential drop and looses the gained energy by non-thermal radiation.

In this case the particles have to pile up which means that their velocity has to decrease. Where the density goes down, the velocity has to increase accordingly to keep the current stable ($`\mathrm{div}\stackrel{}{j}=0`$).

Then the resulting electric field is able to accelerate particles and varies with half the wavelength of the density variation (as follows directly from the Poisson equation).

Such a system is similar to a battery whose poles are the region of the open field lines (acting as the negative pole in the parallel rotator) and the edge of the polar cap being the anode. If a particle starts close to the pole and returns on the last closed field line it encounters a voltage drop of typically $`10^{13}`$ V. But, of course, this does not mean that the particle has to be accelerated by the full potential in one step anywhere along its path.

As will be motivated in sect. 7 it is likely that the outer magnetosphere itself provides a constant minimum resistance which clearly dominates all plasma induced resistances so that the latter are only minor corrections. In other words, if the main voltage drop occurs close to the light cylinder almost all of the energy is dissipated there and transformed into high energetic radiation and ultrarelativistic particles. A schematic view of the model is shown in fig. 1..

Of course, there are two basic assumptions in this model, which have to be checked carefully. On the one hand, the pulsar is described as a voltage source. This description requires that the generator’s internal resistance is negligible compared to the total resistance. It will be shown in section 6 that the maximum current allowed would be a highly relativistic Goldreich-Julian current over the entire two polar caps. The corresponding resistance is about $`5\mathrm{\Omega }`$ which is huge in comparison with the transverse resistance of about $`10^8\mathrm{\Omega }`$ (Itoh, 1975) Thus the first assumption appears to be quite reasonable.

The second assumption, namely the one-dimensional current flow in the inner magnetosphere is also very natural. The conductivity can be expressed by a mean collision time $`\tau `$ using the formula

$$\sigma =\frac{ne^2\tau }{m_\mathrm{e}\epsilon _0}.$$
(1)
The transverse collision time $`\tau `$ can be estimated by $`1/\omega _\mathrm{c}`$ (the inverse gyro frequency) whereas the collision time for the parallel current cannot exceed the inverse plasma frequency $`1/\omega _\mathrm{p}`$ (see below). As in the inner magnetosphere $`\omega _\mathrm{p}\omega _\mathrm{c}`$, the one-dimensional approach is appropriate. Another way of confirming this is to calculate the maximum (transverse) drift velocity in the inner magnetosphere. Taking the maximum plasma wave electric field (see below) and typical pulsar parameters ($`P=0.5`$ s, $`B_0=10^8`$ T) we find a drift velocity of

$$v_{\mathrm{drift}}=0.359\frac{\mathrm{m}}{\mathrm{s}}x_{\mathrm{em}}^{3/2}\gamma ^{1/2}$$
(2)
which is far below $`c`$ for typical emission heights and particle energies used in the model ($`x_{\mathrm{em}}50`$, $`\gamma =10`$). Particles move almost perfectly along the field lines as long as relativistic parallel outflow is possible (on the neutral line the situation is a bit different as a pure electron flow cannot cross this border easily).

An important consequence of our description is that the total energy a particle radiates along its path is in no way limited by its own kinetic energy as the actually dissipated energy comes from the inductor, i. e. the battery. Here, the analogy to a battery with a bulb attached through copper wires can be helpful. In a metal conductor the current carrying electrons have velocities of only some millimeters per second which means a kinetic energy of only $`10^{16}`$eV because of the small path length of electrons in a solid state conductor. Nevertheless each of them converts a total of 1.5 eV into heat on its way from the cathode to the anode. The only restriction which is given by the kinetic energy is the potential difference of the electrostatic field in one elementary radiation cell. As we will point out in sect. 4 this limit directly leads to the maximum possible resistivity in the plasma and therefore is never exceeded anyway.

The accelerating fields occurring in the outer magnetosphere do not invalidate the current circuit model, however as the circuit is open (high resistance) near the light cylinder.

Another advantage of this description is that our model uses particle inertia as the source for an accelerating field near the pulsar surface to get relativistic particles. So a very small acceleration region has to occur as the electrons have non-relativistic velocities inside the neutron star; therefore an overdense thin charge layer develops just above the surface where the electrons reach their terminal energy.

Therefore is no need for a ”gap” (meaning a starved region near the neutron star) but inertia naturally provides the fields necessary to reach the mildly relativistic energies, and the terminal energy is determined by the global current instead of some large-scale space charge field (cf. sect. 6).

## 3 Application to radio emission

In the following we derive a quantitative model for the radio emission of pulsars. The particular aim is to explain the energetics of the radiation, from observed or observationally derived pulsar parameters such as

* the period $`P`$
* the dipolar surface magnetic field $`B`$
* the brightness temperature $`T_\mathrm{B}`$ of the low frequency emission
* the emission height $`x_{\mathrm{em}}`$
* the total luminosity in the radio band $`L`$

With these input parameters we derive the radial width of the emission region, which should be small to be consistent with observations. A narrow radiation zone is suggested by both the highly structured emission in space and time as seen by many microstructure experiments (Boriakoff, 1992).

### 3.1 Plasma beams and coherent emission

Emission processes in relativistic plasmas are always connected with the electron plasma frequency, which is (Kunzl et al. 1998a; Melrose and Gedalin 1999)

$$\omega _{\mathrm{pe}}=\sqrt{\frac{ne^2}{\gamma \epsilon _0m_e}}$$

As a very promising candidate for the origin of coherent radio emission of pulsars we consider the collective inverse Compton mechanism, which was detected first in a number of laboratory experiments in which relativistic electron beams with a number density $`n_\mathrm{b}`$ penetrate into a background plasma with density $`n_\mathrm{p}`$ and in which extremely intense coherent microwave emission has been observed (Kato et al. 1983; Levron et al. 1987; Yoshikawa et al. 1993; 1994). The detected radiation has been explained in terms of collective inverse Compton scattering of strong Langmuir turbulence (Benford and Weatherhall 1992; Benford 1992). Relativistic electron beams interact with self-excited strong density fluctuations in the background plasma if the beam-plasma-density ratio $`n_\mathrm{b}/n_\mathrm{p}`$ exceeds 0.1 for a beam Lorentz factor of about 3. The beam electrons are bunched by the electric field of the density fluctuations and radiate coherently in their reference system at $`\gamma \omega _{\mathrm{pe}}`$, transformed into the observer’s frame this gives an emitted frequency of about $`\gamma ^2\omega _{\mathrm{pe}}`$. The efficiency of this energy transfer of beam energy into photon energy is up to 30%.

The general scaling law from laboratory experiments for the relation between beam Lorentz factor and density ratio is (Benford 1992)

$$\gamma _{\mathrm{crit}}=\left(1+\frac{9}{20}\mathrm{lg}\left(\frac{n_\mathrm{b}}{n_\mathrm{p}}\right)\right)^1$$

As shown in Fig. 2 for large beam energies the required density ratios is quite small, whereas for small beam energies the beam density must be comparable to that of the background plasma to drive collective inverse Compton scattering. For a given beam energy one can calculate the required density ratio.

The underlying physics of this mechanism has been investigated in all details and the general consensus is, that a relativistic electron beam with sufficient energy and particle density penetrates into a background plasma the collective inverse Compton scattering is the unavoidable coherent radiative power output (Benford 1992 and references therein). Based on these findings we transfer this process into a neutron star magnetosphere and investigate, whether this process matches the required properties.

Collective inverse Compton scattering is a nonlinear scattering processes whose radiative output in first order does not depend on the origin of the density fluctuations. In contrast, direct excitation mechanisms like maser or beam radiation depend sensitively on the form of the beam energy distribution and the time evolution of the beam (Ursov & Usov 1988; Melrose & Gedalin 1999; Melrose et al. 2001).

For comparison with alternative processes, we give a general scaling of the emitted frequency with the nonrelativistic plasma frequency $`\omega _{\mathrm{pe}}^{(0)}`$ for the different radiation mechanisms

$$\nu =\frac{\omega _{\mathrm{pe}}^{(0)}}{2\pi }\gamma ^{\alpha /2}.$$
(3)
In this notation the value of $`\alpha `$ for the collective inverse Compton scattering is 3, whereas for direct emission it is 1. In the next chapter we will debate an argument why the description of the pulsar system as a resistive current circuit makes $`\alpha =1`$ less likely than a higher value. Inserting the expressions for the Goldreich-Julian density and the plasma frequency, we can solve this equation for $`\gamma `$:

$$\gamma =\left[\frac{\omega _\mathrm{p}(x_{\mathrm{em}})}{2\pi \nu }\right]^{2/\alpha }=\left(\frac{B_0\mathrm{\Omega }e}{2\pi ^2\nu ^2m_\mathrm{e}x_{\mathrm{em}}^3}\right)^{1/\alpha }.$$
(4)
According to the usual estimates, an incoherent process from an electron with a Lorentz factor of $`\gamma `$ cannot produce radiation with a brightness temperature exceeding $`10^{10}\mathrm{K}\gamma `$. However, in the case discussed here, this argument is not applicable, because there is no absorption. The emitted power of a radiation process is not limited any more if the radiation produced further inside is not absorbed on its way out of the magnetosphere.

Consider a cylindrical plasma column of radius $`R`$, length $`l`$, particle density $`n`$ and single particle power output $`P_{\mathrm{single}}`$ in forward direction. Then the total power adds up to

$$P_{\mathrm{ges}}=lR^2\pi nP_{\mathrm{single}}.$$
(5)
Thus the brightness temperature can reach arbitrarily high values on the condition that $`l`$ is large enough no matter, whether the emission process is coherent or not.

Consequently we use a different estimate for the brightness temperature. The average flux is received from an area that corresponds to the polar cap moved and expanded along the field lines up to the emission height:

$$A=A_{\mathrm{cap}}x_{\mathrm{em}}^3.$$
(6)
The total radio luminosity over a frequency band of $`\mathrm{\Delta }\nu `$ around the frequency $`\nu `$ can therefore be equated to the thermal flux of a blackbody with the appropriate brightness temperature. For broadband emission ($`\nu =\mathrm{\Delta }\nu `$) and using the definition of the brightness temperature we find

$$L=F_\nu \mathrm{\Delta }\nu A_{\mathrm{cap}}x^3\frac{2\pi \nu ^3}{c^2}k_\mathrm{B}T_\mathrm{B}A_{\mathrm{cap}}x_{\mathrm{em}}^3$$
(7)
which yields a brightness temperature of

$`T_\mathrm{B}`$ $`={\displaystyle \frac{Lc^2}{2\pi \nu ^3k_\mathrm{B}A_{\mathrm{cap}}x_{\mathrm{em}}^3}}=`$
$`=`$ $`1.8710^{23}\mathrm{K}\left({\displaystyle \frac{L}{10^{20}\mathrm{W}}}\right)\left({\displaystyle \frac{\nu }{400\mathrm{MHz}}}\right)^3\left({\displaystyle \frac{P}{\mathrm{s}}}\right)\left({\displaystyle \frac{x_{\mathrm{em}}}{50}}\right)^3.`$
Note that there is no additional geometric factor as the estimate of the total radio luminosity already averages over some typical pulse profile. The values for the brightness temperature found in (3.1) are consistent with observationally deduced results (Sutherland, 1979, Kramer, 1995) as the same principal method has been used in both cases.

High coherence is required for a different reason, namely that the radiating volume is small enough. In the next step we show that our model is able to reproduce these fluxes under the assumption of a thin layer emitting highly coherent radiation, but where the coherence factor is kept below the limit imposed by the model. As already mentioned, we propose a mechanism of the relativistic plasma emission type.

For the discussed process, particles are scattered by strong, nonlinear plasma waves (solitons) excited by the two-stream instability. Such a radiation mechanism can be described as an inverse Compton scattering (ICS) of solitons by relativistic electrons.

For each particle coherent ICS radiates a power of

$$P_{\mathrm{rad}}^{\mathrm{coh}}=N\sigma _\mathrm{T}c\frac{E_{\mathrm{wave}}^2\epsilon _0}{2}\gamma ^2$$
(9)
with the wave electric field $`E_{\mathrm{wave}}`$ if $`\mathrm{}\omega _{\mathrm{wave}}m_\mathrm{e}c^2`$, which is the case for radio emission<sup>1</sup><sup>1</sup>1For higher wave energies the cross section is no longer described by $`\sigma _\mathrm{T}`$, so that the Klein-Nishina cross section ($`\sigma _{\mathrm{KN}}\sigma _\mathrm{T}(\mathrm{ln}\gamma )/\gamma `$) must be used.. $`N`$ is the coherence number, and $`\gamma ^2`$ comes from the Lorentz transfomation of $`E_{\mathrm{wave}}`$. As long as the wave electric field is small, the wave can be treated linearly as a small density fluctuation. But in an unstable situation the wave grows exponentially. The growth saturates quickly because the wave will also accelerate particles as a back reaction. Thus, the wave growth has to be calculated using non-linear methods (Benford & Weatherall, 1992). As a result, the strongest possible wave has an electrostatic field energy comparable to the kinetic energy of the plasma. In the extreme, these two energy densities are equal

$$n_{\mathrm{GJ}}\gamma m_\mathrm{e}c^2=\frac{1}{2}\epsilon _0E_{\mathrm{wave}}^2$$
(10)
which means (9) can be rewritten as

$$P_{\mathrm{rad}}^{\mathrm{coh}}=N\sigma _\mathrm{T}cn\gamma ^3m_\mathrm{e}c^2.$$
(11)
As particles are not assumed to change their longitudinal momentum in the radio emission region, we can interpret the radiated power as an electric dissipation field (which is the field necessary to balance radiation losses). Its strength is

$$E_{\mathrm{diss}}=\frac{P_{\mathrm{rad}}^{\mathrm{coh}}}{ec}=\frac{2N\sigma _\mathrm{T}\mathrm{\Omega }B_0\epsilon _0m_\mathrm{e}c^2}{e^2x_{\mathrm{em}}^3}\gamma ^3.$$
(12)
### 3.2 Extension of the radiation region and the maximum luminosity

Now we calculate the radial width of the radiation zone. This is an important test of the model, as the very small RFM taken from the observations requires a fairly narrow emission zone. Apart from that, the highly modulated radiation observed in the radio band would be hard to explain for large widths.

In (12) the only the free parameter left is $`N`$. However, this number cannot directly be fixed from the brightness temperature as mentioned before. Nevertheless there are upper limits for the coherence factor which serve to derive a minimum radial width of the emission region.

As we consider a plasma wave being the cause of coherence, the maximum volume of coherently radiating particles is a sphere of around one Debye volume. Taking into account relativistic effects, Melrose (1992) finds a maximum coherence volume of

$$V_{\mathrm{coh}}=\left(\frac{c}{\omega _\mathrm{p}}\right)^3\frac{\gamma ^2}{\pi }$$
(13)
and therefore

$`N_{\mathrm{max}}`$ $`=`$ $`n_{\mathrm{GJ}}\left({\displaystyle \frac{c}{\omega _\mathrm{p}}}\right)^3{\displaystyle \frac{\gamma ^2}{\pi }}=\sqrt{{\displaystyle \frac{m_\mathrm{e}^3c^6\epsilon _0\gamma ^7x_{\mathrm{em}}^3}{2\mathrm{\Omega }}}B_0e^5}=`$ (14)
$`=`$ $`6.4110^{13}\gamma ^{7/2}\left({\displaystyle \frac{P}{\mathrm{s}}}\right)^{1/2}B_8^{1/2}\left({\displaystyle \frac{x_{\mathrm{em}}}{50}}\right)^{3/2}.`$
Another upper limit is found by equating the particle and the field energy density, since in the optimal case those quantities are comparable. More detailed calculations (Lesch & Schlickeiser, 1987) show that the field energy density is less than half of the particle energy density (as already used in eq. (10)).

Inserting this result into (12) we find a maximum coherence number of

$`N`$ $`=`$ $`\sqrt{{\displaystyle \frac{e^3x_{\mathrm{em}}^3}{2\gamma ^5\mathrm{\Omega }B_0m_\mathrm{e}c^2\epsilon _0^2\sigma _\mathrm{T}^2}}}=`$ (15)
$`=`$ $`3.7910^{15}\gamma ^{5/2}\left({\displaystyle \frac{P}{\mathrm{s}}}\right)^{1/2}B_8^{1/2}\left({\displaystyle \frac{x_{\mathrm{em}}}{50}}\right)^{3/2}.`$
Comparing the maximum coherence numbers from (14) and (15) we find that the latter is smaller whenever $`\gamma 1.97`$. Therefore it is a good approximation to use (15) for determining the maximum coherence number.

For estimates of the radial extension of the radio emission region we take the total radio power and calculate the radiating volume necessary to reproduce the observed luminosities.

The total energy loss in a resistive current reads

$$L=\underset{V_{\mathrm{rad}}}{}E_{\mathrm{diss}}jdV$$
(16)
where $`V_{\mathrm{rad}}`$ denotes the dissipation volume.

$`E_{\mathrm{diss}}`$ is given by (12) and (15), $`j(x)=n_{\mathrm{GJ}}(x)ec`$, and by assuming that the non-ideality extends over the entire cross-section of open field lines; using a dipolar geometry and neglecting angular corrections to the GJ- density all quantities in (16) only depend on the radial coordinate. The volume element d$`V`$ can be transformed to

$$\mathrm{d}V=r_{\mathrm{NS}}A_{\mathrm{cap}}x^3\mathrm{d}x.$$
(17)
Here we have used (6), as the cross-section is a simple function of $`x`$. Therefore (16) reduces to a one-dimensional integral.

Inserting all the previous expressions the integral reads

$`L`$ $`=`$ $`\left({\displaystyle \frac{2\pi ^2\nu ^2m_\mathrm{e}}{B_0\mathrm{\Omega }e}}\right)^{1/(2\alpha )}\sqrt{{\displaystyle \frac{2\mathrm{\Omega }B_0m_\mathrm{e}c^2}{e}}}{\displaystyle \frac{2\mathrm{\Omega }B_0\epsilon _0}{e}}ec`$ (18)
$`r_{\mathrm{NS}}{\displaystyle \frac{\pi \mathrm{\Omega }r_{\mathrm{NS}}^3}{c}}{\displaystyle \underset{x_{\mathrm{em}}}{\overset{x_{\mathrm{em}}+\stackrel{~}{l}_0}{}}}x^{3/(2\alpha )}x^{3/2}dx`$
where $`\stackrel{~}{l}_0`$ denotes the radial extension of the radiation region in units of $`r_{\mathrm{NS}}`$, whereas $`x_{\mathrm{em}}`$ stands for the height of its inner edge in pulsar radii.

Evaluating the integral and inserting numbers for a frequency of $`\nu =400`$ MHz we find

$$L=5.5610^{22}\mathrm{W}\left(\frac{P}{\mathrm{s}}\right)^{1/\alpha 5/2}B_8^{3/21/\alpha }F_\alpha $$
(19)
where

$$F_\alpha =0.0286^{1/\alpha }\{\begin{array}{cc}\stackrel{~}{l}_0& \mathrm{for}\alpha =1\\ 4\left[\left(x_{\mathrm{em}}+\stackrel{~}{l}_0\right)^{1/4}x_{\mathrm{em}}^{1/4}\right]& \mathrm{for}\alpha =2\\ \mathrm{ln}\frac{x_{\mathrm{em}}+\stackrel{~}{l}_0}{x_{\mathrm{em}}}& \mathrm{for}\alpha =3\end{array}$$
(20)
Solving this equation for $`\stackrel{~}{l_0}`$ in the case $`\alpha =3`$ renders

$`\stackrel{~}{l_0}=x_{\mathrm{em}}`$ (21)
$`\left\{\mathrm{exp}\left[5.8810^3\left({\displaystyle \frac{P}{s}}\right)^{13/6}B_8^{7/6}\left({\displaystyle \frac{L}{10^{20}\mathrm{W}}}\right)\right]1\right\}`$
which, for typical parameters ($`L=10^{20}`$ W, $`P=1\mathrm{s}`$, $`B_8=1`$, $`x_{\mathrm{em}}=50`$) yields $`\stackrel{~}{l_0}=0.295`$ which justifies the assumption of a narrow emission region. The dependence of $`\stackrel{~}{l_0}`$ and $`L`$ is shown in fig. 3.

From the result above one can obtain sufficient luminosity even for slow pulsars. However, for small $`\mathrm{\Omega }`$ the radial extension of the emission region increases. Especially, the corresponding electric field (12) can even exceed the vacuum field ($`E_{\mathrm{vac}}=\mathrm{\Omega }B_0r_{\mathrm{NS}}x^4`$) locally. In a shielded magnetosphere the kinetic energy density and by equipartition the dissipated electric field, may exceed the energy density of the vacuum field.

## 4 Application to PSR 0329+54

The pulsar PSR 0329 is an extensively researched object and one of the few where the spatial orientation of magnetic and spin axis is well determined (Rankin, 1983). Assuming that the profile peaks occur on the same field lines for all frequencies enabled Mitra and Rankin (2002) to obtain revised effective opening angles and emission heights for PSR 0329+54 amongst others. A good fit to the surface opening angle (including relativistic corrections) for the field lines that provide the two outer peaks of the profile is $`\varphi _0=0.66^\mathrm{o}`$. For any angle $`\varphi `$ between those peaks we can find an emission height of a profile at a specific frequency using the parametric equation of dipole field lines:

$$r(\varphi )=\frac{r_{ns}}{f_1}\frac{\mathrm{sin}(\varphi )}{\mathrm{sin}(\varphi _0)}$$
(22)
Mitra and Rankin (2002) found that the emission heights range from about 500 km for frequencies above 1 GHz up to 1100 km for the lowest frequency of 100 MHz. Together with the known radio fluxes (Kramer et al. 1997) we were able to match the described model to the data and determine the dissipation length scales for frequencies between 100 MHz and 30 GHz. To enhance the realism of the fit, the relativistic (Muslimov & Tsygan, 1992) corrections to the Goldreich-Julian densities and effective areas were also included. In this case we calculated the local Goldreich-Julian density through

$$n(r)=2ϵ_0B_0\left(\frac{r_{\mathrm{NS}}}{r}\right)^3\frac{f(r)}{f_1\sqrt{1\frac{r_\mathrm{g}}{r}}}$$
(23)
Here

$$f(r)=3(\frac{r}{r_\mathrm{g}})^3\left(\mathrm{ln}(1\frac{r_\mathrm{g}}{r})+\frac{r_\mathrm{g}}{r}(1+\frac{1}{2}\frac{r_\mathrm{g}}{r})\right)$$

and $`f_1=f(r_{\mathrm{NS}})`$ are the corrections due to general relativity (Muslimov & Tsygan). In this context $`r_\mathrm{g}=\frac{2\mathrm{\Gamma }m_{ns}}{c^2}`$ denotes the Schwarzschild radius of the neutron star.

Assuming that the observed emission is indeed caused by coherent inverse Compton scattering, we can calculate the Lorentz factor of the emitting electron current as a function of distance from the pulsar.

$$\gamma (\omega _{\mathrm{obs}})=\left(\frac{ϵ_0m_{\mathrm{NS}}}{n(r_{\mathrm{em}}(\omega _{\mathrm{obs}}))e^2}\right)^{\frac{1}{3}}\omega _{\mathrm{obs}}^{\frac{2}{3}}$$
(24)
It varies from 13 (100 MHz) to 267 (30 GHz). Using the relativistically corrected values for distance, density and Lorentz factor we obtain an improved estimate of the dissipation scale $`\lambda `$ which is given by

$$\lambda =\frac{L}{P_{\mathrm{coh}}(r,\gamma )n(r)2\pi r^2\left(1\mathrm{cos}(\varphi )\right)}$$
(25)
Here $`L`$ and $`\varphi `$ are the observed luminosity and opening angles of the pulsar beam at a given frequency and $`P_{\mathrm{coh}}`$ means the emitted power per coherently radiating particle. We estimated the luminosity as an upper limit from the radio flux by assuming that the smaller cap area given by the peak flux opening angle is filled with uniform emission can be a good representation of the conal beam structure.

Fig. 4 shows the development of the the dissipation scale with emission height. After a steep rise at 500 km ($`\nu >1000`$ MHz), $`\lambda `$ stays around 20 km up to a height 1100 km (100 MHz). Thus we have evidence for a quick onset of a strong and stable coherence at a height of 600 km, lasting for another 600 km. This is also borne out by the fact that the brightness temperatures of the source vary only a little in this regime where most of the pulsar radio luminosity is created. A typical fluctuation time scale of $`\gamma \lambda /c=1\mathrm{m}\mathrm{s}`$ is associated with the dissipation scale below 1 GHz. Above that frequency it is expected to decrease monotonically to $`20\mu \mathrm{s}`$ at 20 GHz. These timescales are in good agreement with observational evidence (Lange et al. 1998) where typical microstructure timescales were found to range within $`6001200\mu s`$ for that pulsar at 4.85 GHz. There is also some observational evidence for a decrease in microstructure timescales with increasing frequency for the pulsars PSR B1133+16 and PSR B2016+28 (Lange et al. 1998).

## 5 Micropulse flux enhancement by beaming effects

In this section we discuss anisotropy effects for pulse substructures which may enhance the observed brightness temperatures by up to several orders of magnitude. It is very unlikely that the flux enhancement only arises from a growth of the dissipation region. For the nanosecond giant pulses of the Crab pulsar this is even impossible (if the width of the radio emission region jumped from only 0.001 $``$ which is extremely narrow $``$ to around 1 pulsar radius, the peak would occur on a timescale of $`\tau _{\mathrm{mic}}r_{\mathrm{NS}}/c310^5`$ s). But as such small time scales also mean highly localized emission, the anisotropy of the elementary radiation process has to be considered. For an isotropic process there is a simple connection between the radiated power and the observed flux:

$$S_{\mathrm{obs}}(\nu )=\frac{1}{4\pi d^2}\frac{}{\nu }P_{\mathrm{rad}}(\nu )$$
(26)
where $`P_{\mathrm{rad}}(\nu )`$ is the total power emitted in frequencies below $`\nu `$, whereas $`d`$ denotes the transverse extension of the radiation region. For a rough quantitative estimate, we can set $`P_{\mathrm{rad}}(\nu _{\mathrm{crit}})=L`$, and replace the derivative by a quotient. Thus we obtain the simpler equation

$$S_{\mathrm{obs}}(\nu _{\mathrm{crit}})=\frac{L}{4\pi \nu _{\mathrm{crit}}d^2}.$$
(27)
This flux is drastically enhanced by two effects, producing a strongly anisotropic emission pattern. The first is the well-known relativistic lighthouse effect (Rybicki & Lightman, 1979). A relativistically moving source of isotropic radiation emits most of the power in a narrow cone with an angular opening of about $`1/\gamma `$ for $`\gamma 1`$, as can easily be verified by performing the Lorentz transformation. The exact result is

$$\theta ^{}=\mathrm{arctan}\left(\frac{1}{\gamma }\frac{\mathrm{sin}\theta }{\left(\mathrm{cos}\theta +\beta \right)}\right)$$
(28)
with $`\theta `$ being the angle between the direction of the flux and source velocity in the co-moving frame and $`\theta ^{}`$ denoting this angle in the observer’s system. For $`\theta =90^\mathrm{o}`$ and $`\gamma 1`$ the approximation $`\theta ^{}1/\gamma `$ is a good one.

Therefore we find a flux enhancement of

$$\xi _1=\frac{4\pi }{\pi /\gamma ^2}=4\gamma ^2$$
(29)
by the lighthouse effect.

Another anisotropy factor is the coherent radiation process itself, as has been shown by Kunzl et al. (1998b). To understand this effect, one has to recall that coherent emission requires a phase coupling of the emitting particles in one direction. For relativistic plasma emission, the preferred direction obviously is the direction of the streaming velocity (forward direction).

As long as the spatial dimension of the coherently radiating volume is small compared to the emitted wavelength, the interference is still constructive even under large angles to the forward direction. However, as soon as the extension of the coherence region becomes comparable to the wavelength, the coherent emitter can be seen as a phased antenna field producing intrinsically beamed emission.

The quantitative analysis of this effect has been performed by Kunzl et al. (1998b). The first minimum appears under the angle

$$\alpha _{\mathrm{min}}=\sqrt{1+\frac{2\pi }{T}}1.$$
(30)
Here $`T:=2\pi d/\lambda `$ with $`d`$ being the transverse extension of the emission region. The angle of the minimum is much smaller than unity for $`T1`$.

Similar to the lighthouse effect, we can use this angular extension to estimate the enhancement factor by coherence beaming. Thus the flux grows by a factor of

$$\xi _2=\frac{4\pi }{\pi \alpha _{\mathrm{min}}^2}=\frac{4}{\left(\sqrt{1+\frac{2\pi }{T}}1\right)^2}.$$
(31)
A more realistic treatment of the geometry uses a cylindrically symmetric soliton with a density profile proportional to $`\mathrm{cosh}^2(R_0/R)`$ where $`R_0\lambda _\mathrm{D}`$ is the typical length scale. Numerical calculations for this structure show a much stronger anisotropy (Kunzl et al. 1998b)

Combining the two anisotropy effects we find the total beaming factor

$$\xi :=\xi _1\xi _2=16\gamma ^2\left(\sqrt{1+\frac{2\pi }{T}}1\right)^2$$
(32)
which, even with quite moderate parameters (like $`\gamma =10,x_{\mathrm{em}}=50,N=10^{13}`$ for a typical pulsar), causes an enhancement of some $`10^3`$. Therefore it can easily explain strong substructures, and micropulse brightness temperatures that are some $`10^3`$ times above the mean value (Boriakoff, 1992).

As these beaming effects apply only to the elementary emission process, but do not enhance the mean flux significantly, it is an important check on whether the model can reproduce the observed brightness temperatures of up to $`10^{31}`$ K (Hankins, 1996) observed in giant pulses of the Crab pulsar.

Therefore we take the maximum possible number of coherently radiating particles, and the corresponding beaming effects. Let $`A`$ be an arbitrary cross section. Then starting with eq. (3.1) one can express the integrated flux $`I:=I_\nu \mathrm{\Delta }\nu =L_{\mathrm{radio},\mathrm{A}}/A`$ by the brightness temperature:

$$I:=\frac{L_{\mathrm{radio},\mathrm{A}}}{A}=\frac{2\pi \nu ^3k_\mathrm{B}T_\mathrm{B}}{c^2}.$$
(33)
Inserting the critical frequency of the Crab pulsar $`\nu =\nu _{\mathrm{crit}}=160`$ MHz and $`T_\mathrm{B}=10^{31}`$ K, the integrated flux is

$$I=3.9510^{16}\frac{\mathrm{W}}{\mathrm{m}^2}.$$
(34)
As this value means the beamed flux, the actual power per area is only

$$I_{\mathrm{real}}:=\frac{I}{\xi }=2.4710^{15}\frac{\mathrm{W}}{\mathrm{m}^2}\gamma ^2\left(\sqrt{1+\frac{2\pi }{T}}1\right)^2.$$
(35)
Now we can compute the minimum thickness of a layer producing the power per cross-section. To obtain the lower limit we take the strongest dissipation field possible (see eq. (12)) and assume a coherence cell with a lateral extension of $`c/\nu _\mathrm{p}`$, so that $`T=2\pi \gamma `$.

The total dissipated power from a volume with the cross-section $`A`$ and the radial extension $`d`$ is

$$P_{\mathrm{diss}}\stackrel{!}{=}I_{\mathrm{real}}A=E_{\mathrm{diss}}djA.$$
(36)
For a relativistic Goldreich-Julian current $`j=j_{\mathrm{GJ}}=n_{\mathrm{GJ}}(x)ec`$ eq. (36) can be solved for $`d`$. Inserting the Crab pulsar parameters ($`P=33.4`$ ms, $`B=3.810^8`$ T, $`\nu =160`$ MHz, $`x_{\mathrm{em}}=80`$) and additionally using eq. (4) we obtain

$$d=3.85\mathrm{m}\left(\frac{T_\mathrm{B}}{10^{31}\mathrm{K}}\right)$$
(37)
which corresponds to a time scale of

$$\tau :=\frac{d}{c}=1.2810^8\mathrm{s}\left(\frac{T_\mathrm{B}}{10^{31}\mathrm{K}}\right)=12.8\mathrm{ns}\left(\frac{T_\mathrm{B}}{10^{31}\mathrm{K}}\right).$$
(38)
Thus we can expect to see giant pulses on nanosecond timescales, with brightness temperatures above $`10^{30}`$ K, from far inside the light cylinder in the Crab pulsar, although the relative emission height will be considerably larger than for average pulsars.

## 6 Densities in the radiation region

Here we briefly discuss if the model described above could also work for particle densities significantly above or below the GJ value. We will show that even if there were a mechanism producing a current density and particle energies matching the requirements to produce the low radio frequencies, such a scenario would be very unlikely.

Assuming that the particle density was much larger than the value predicted by Goldreich & Julian (1969), the Lorentz factor corresponding to a certain frequency would drop according to (4). Especially for young pulsars this would quickly require $`\gamma <1`$ which is, of course, not possible (even $`\beta \gamma 1`$ should be fulfilled as the particles are definitely relativistic). This important fact has already been pointed out by Kunzl et al. (1998a).

A relativistic current significantly below the GJ current (meaning that the particle density is smaller than $`n_{\mathrm{GJ}}`$ faces even more difficulties. On the one hand, accelerating fields are incompletely shielded (the equilibrium Lorentz factor where acceleration and radiation losses $``$ now including coherent radio emission $``$ balance is extremely high, so coherence is highly doubtful). Furthermore the one-dimensional approach would no longer account as a permanent strong deviation from ideality allows efficient $`\stackrel{}{E}\times \stackrel{}{B}`$ drift. But also the existence of strong $`\gamma `$ -pulsars provides an argument against this solution:

As we observe pulsed $`\gamma `$-emission of pulsars whose high energetic luminosity reaches up to 10% of the spindown luminosity, some mechanism must efficiently convert particle energy into (hard) radiation (Kanbach, 2001). But the potential difference a particle can suffer is limited by the polar cap potential. Since the power dissipated by a relativistic GJ current flowing over the entire polar cap area with the polar cap voltage is comparable to the spindown power of the neutron star, a large fraction (well above 10%) of the generator’s power must be pumped into particles. This, of course, is only possible for currents of about the GJ current and thus for particle densities comparable to the Goldreich-Julian value.

## 7 Discussion

We have presented the outline of a model which may explain the creation of radio emission within the energy constraints of the observationally deduced emission heights of 50-100 pulsar radii. (e.g. Blaskiewicz et al. 1991; Taylor et al. 1993; Kijak and Gil 1997,1998; Kramer et al. 1997). Several other pulsar models fail to explain this feature, e.g. the model of Lyutikov et al. (1999) which uses Cherenkov drift resonance, but predicts the radiation zone to be in some 1000 pulsar radii.

One of the main differences to the standard model is that we do not assume an inner gap which accelerates particles to ultrarelativistic energies and is saturated by a pair avalanche. Instead we suggest a conducting wire model with a quasi- stationary, low relativistic current that suffers fluctuations on short time scales. Such a model is supported by calculations on particle densities available from the polar cap by thermal and field emission (Jessner et al. 2001) which show that inner gap models have severe difficulties in producing particles energetic enough to start pair cascades for parallel or nearly parallel rotators.

We showed that the brightness temperatures for the mean pulse can be reproduced by our model and also micropulse fluxes can be explained by appropriate beaming assumptions.

The radiation process in the framework of such a model could be some relativistic plasma emission like the ALAE process described by Melrose (1978). This mechanism requires mildly relativistic particles to be efficient. As shown also in Kunzl et al. (1998a) any plasma process responsible for the low frequency emission in the Crab pulsar must work with mildly relativistic particles. Von Hoensbroech et al. (1998) showed that in a Goldreich-Julian magnetosphere with low relativistic particles propagation effects can also qualitatively explain the observed polarization properties.

Finally we would like to note that our model does not rule out pair production in general. As can be seen from several observations of the Crab pulsar a dense particle wind is needed which transports energy to the nebula. The calculations of Hirotani and Shibata (1999) show that a significant fraction of the vacuum potential drop can occur in an outer gap if the X-ray photon density is not too high. For a detailled geometric examination of the outer gap structure see Cheng et al. (2000).

It has been shown that for reasonable parameters an outer gap model could easily explain the observed energy and particle fluxes in the pulsar wind. Recently Crusius et al. (2000) have shown that an outer gap model can explain the observed spectrum and energy output of the Crab pulsar from the infrared to X-ray frequencies which provides another hint that the main dissipation region is in the outer magnetosphere.

We gratefully acknowledge constructive discussions with J. Arons, D.B. Melrose M. Kramer, D. Mitra and A. von Hoensbroech. In addition D. Mitra and M. Kramer provided us with access to their radio flux- and emission height measurements for PSR 0329+54 for which we are particularly grateful.


format: text

text:
# 1 Introduction

## 1 Introduction

In recent years, the determination of quark masses has become one of the main research topics of lattice QCD simulations. An accurate determination of these masses is in fact of great importance for both phenomenological and theoretical studies. The masses of the charm and bottom quarks, for instance, enter the theoretical predictions of beauty hadron decay rates which, in turn, are relevant for the phenomenological analysis of the CKM unitarity triangle and thus of CP violation in the Standard Model. On the more theoretical side, a precise knowledge of quark masses may give insight on the physics of flavour, by revealing relations between masses and mixing angles or specific textures of the quark mass matrices, eventually due to still uncovered flavour symmetries.

The values of quark masses cannot be directly measured in the experiments, since quarks are confined inside the hadrons. On the other hand, being fundamental parameters of the theory, quark masses cannot even be computed on the basis of purely theoretical considerations. Their evaluation is based on the comparison of the theoretical estimate of a physical quantity, which depends on quark masses, and its experimental value. This is typically realized on the lattice by using, as experimental input, the values of pseudoscalar or vector meson masses.

In this paper we present the results of an extensive lattice calculation of the strange and the average up/down quark masses in the quenched approximation. Particular attention has been dedicated to the reduction and control of the systematic uncertainties, particularly in the case of the strange quark mass. Leading $`𝒪(a)`$ discretization effects (where $`a`$ is the lattice spacing) have been removed by using the non-perturbatively $`𝒪(a)`$-improved Wilson action and operators . The systematic uncertainty related to the evaluation of the quark mass renormalization constant in perturbation theory has been significantly reduced, by implementing the non-perturbative renormalization technique of ref.  in the $`\mathrm{RI}/\mathrm{MOM}`$ scheme. Conversion of the quark masses from the $`\mathrm{RI}/\mathrm{MOM}`$ scheme to the most popular $`\overline{\mathrm{MS}}`$ scheme has been performed by using continuum perturbation theory at the $`\mathrm{N}^3\mathrm{LO}`$ . Finite volume effects have also been studied. In order to estimate residual systematic uncertainties, we have compared the results obtained by using two alternative definitions of the lattice bare quark mass, related to the vector and axial-vector Ward identities respectively. Finally, with respect to previous determinations of light quark masses by our collaboration , in this study we have performed a calculation at four different values of the lattice scale, corresponding to an inverse lattice spacing in the range between approximately 2 and 4 GeV. In this way, we have been able to extrapolate our results to the continuum limit. The main relevant source of uncertainty, which is left in our calculation of the strange quark mass, is therefore the quenched approximation.

An additional uncertainty is present in the determination of the average up/down quark mass. Typical values of the lightest quark masses, used in the present and most of current lattice calculations, are approximately of the order of $`m_s/2`$, where $`m_s`$ is the strange quark mass. Therefore, a large chiral extrapolation is required to reach the physical values of the up and down quark masses. Chiral perturbation theory may be used as a guidance in this extrapolation, but the inclusion of higher order terms in the chiral expansion, necessary to increase the accuracy of this determination, requires simulations with many more (and preferably lighter) quark masses . In the region of masses considered in this paper, the pseudoscalar meson mass squared shows a good linear dependence on quark masses, and a linear or quadratic fit has been considered in performing the chiral extrapolation. The systematic uncertainty introduced by this extrapolation is difficult to be reliably quantified. It will be assessed only when simulations on larger lattice volumes and smaller values of quark masses become feasible. We stress, however, that this uncertainty does not affect the determination of the strange quark mass.

We conclude this section by summarizing the main results of this paper. For the strange quark mass and the average value of the up/down quark masses, $`m_{\mathrm{}}=(m_u+m_d)/2`$, quoted in the $`\overline{\mathrm{MS}}`$ scheme at the renormalization scale $`\mu =2\mathrm{GeV}`$, we obtain

$$m_s^{\overline{\mathrm{MS}}}(2\mathrm{GeV})=(106\pm 2\pm 8)\mathrm{MeV}$$
(1)
and

$$m_{\mathrm{}}^{\overline{\mathrm{MS}}}(2\mathrm{GeV})=(4.4\pm 0.1\pm 0.4)\mathrm{MeV},$$
(2)
in good agreement with the current lattice world averages . For the ratio of the strange to the average light quark mass we find

$$\frac{m_s}{m_{\mathrm{}}}=24.3\pm 0.2\pm 0.6.$$
(3)
This result is in good agreement with the prediction $`m_s/m_{\mathrm{}}=24.4\pm 1.5`$ based on NLO chiral perturbation theory .

## 2 Details of the lattice calculation

In determining the values of quark masses we used the standard procedures based on the vector and axial-vector Ward identities .

Vector Ward Identity (VWI): in the renormalized continuum theory the VWI reads

$$_\mu V_\mu (x)𝒪^{}(0)=(m_1(\mu )m_2(\mu ))S(\mu ;x)𝒪^{}(0),$$
(4)
where $`m_{1,2}`$ are the quark masses, $`V_\mu =\overline{q}_1\gamma _\mu q_2`$ is the vector current, and $`S=\overline{q}_1q_2`$ is the scalar density. The renormalized quark mass, $`m_q(\mu )=Z_m(\mu a)m_q(a)`$, is obtained from the bare mass $`m_q(a)`$ which, on the lattice with Wilson fermions, is equal to $`m_q(a)=(1/2a)(1/\kappa _q1/\kappa _{cr})`$. $`\kappa _q`$ is the Wilson hopping parameter and $`\kappa _{cr}`$ is its critical value, corresponding to the chiral limit. The VWI relates the quark mass renormalization constant to the one of the scalar density, i.e. $`Z_m=Z_S^1`$. The computation of the quark mass, using the VWI, can be summarized by the following formula:

$$m_q^{(\mathrm{VWI})}(\mu )=Z_m(\mu a)m_q^{(\mathrm{VWI})}(a)=Z_S^1(\mu a)\left(1+b_mam_q\right)m_q(a).$$
(5)
Notice that the parameter $`b_m=b_S/2`$ provides the elimination of $`𝒪(am_q)`$ effects .

Axial-Vector Ward Identity (AWI): the renormalized continuum AWI reads

$$_\mu A_\mu (x)𝒪^{}(0)=2m(\mu )P(\mu ;x)𝒪^{}(0)$$
(6)
where $`A_\mu =\overline{q}\gamma _\mu \gamma _5q`$ is the axial current, $`P=\overline{q}\gamma _5q`$ is the pseudoscalar density, $`𝒪`$ is a generic operator and we have considered quark fields with degenerate masses. The quark mass renormalization constant, in this case, is related by the AWI to the renormalization constant of the axial and of the pseudoscalar operators, i.e. $`\overline{Z}_m=Z_A/Z_P`$. By choosing $`𝒪=P`$ we have:

$$m_q^{(\mathrm{AWI})}(\mu a)=\overline{Z}_m(\mu a)m_q^{(\mathrm{AWI})}(a)=\frac{Z_A}{Z_P(\mu a)}\left(1+(b_Ab_P)am_q\right)\frac{{\displaystyle \underset{\stackrel{}{x}}{}}_0A_0^I(x)P^{}(0)}{2{\displaystyle \underset{\stackrel{}{x}}{}}P(x)P^{}(0)},$$
(7)
where the bare axial current is improved at $`𝒪(a)`$ as $`A_\mu ^I(x)=A_\mu (x)+ac_A_\mu P(x)`$. The coefficient ($`b_Ab_P`$) cancels the terms of $`𝒪(am_q)`$. For the time derivative, we consider the symmetric ($`𝒪(a)`$-improved) form, i.e. $`_0f(t)=(f(t+a)f(ta))/2a`$. The improvement coefficients $`b_m`$, $`b_Ab_P`$ and $`c_A`$, in eqs. (5) and (7), are only functions of the bare lattice coupling $`g_0^2`$.

Complete information about the lattice calculation performed in this study is provided in table 1. We generated $`𝒪(1000)`$ gauge field configurations in the quenched approximation at four values of the gauge coupling constant, corresponding to an inverse lattice spacing in the range between approximately 2 and 4 GeV. For each value of the lattice spacing, quark propagators have been computed at four light values of the bare quark mass, by using the non-perturbatively $`𝒪(a)`$-improved Wilson action .

In view of the final extrapolation of the lattice results to the continuum limit, an important requirement in this calculation is a precise determination of the lattice spacing which corresponds to the different values of the coupling used in this study. While the absolute value of the physical scale is affected by a rather large uncertainty in the quenched approximation (of the order of 10%), the ratio between two scales can be determined with better accuracy. To this purpose, we use the precise determination based on the study of the static quark anti-quark potential , which in the range $`5.7\beta 6.92`$ can be expressed in the form

$$\mathrm{ln}\left(\frac{a^1(\beta )}{a^1(\beta =6)}\right)=1.7331(\beta 6)0.7849(\beta 6)^2+0.4428(\beta 6)^3.$$
(8)
By using this formula, and varying the inverse lattice spacing at the reference point $`\beta =6.0`$ in the conservative range between 1.9 and 2.1 GeV, we obtain the estimates of the scale given in table 1.

In the same table, we also give the results for the relevant renormalization constants, $`Z_A`$, $`Z_S`$ and $`Z_P`$, in the chiral limit, which have been determined by using the non-perturbative renormalization method of ref. , in the $`\mathrm{RI}/\mathrm{MOM}`$ scheme. The scale dependent constants, $`Z_S`$ and $`Z_P`$, have been computed at the scale $`\mu =3\mathrm{GeV}`$, which lies in the allowed range $`\mathrm{\Lambda }_{QCD}<\mu <\pi /a`$ for all the values of the coupling considered in this study. For this reason, our non-perturbative results for the quark masses in the $`\mathrm{RI}/\mathrm{MOM}`$ scheme will be given at the reference scale $`\mu =3\mathrm{GeV}`$. Details of the non-perturbative calculation of the renormalization constants have been presented at the “Lattice 2002” conference and will be discussed in a separate publication .

Concerning the values of the improvement coefficients, $`c_A`$, $`b_m`$ and ($`b_Ab_P`$), we use the non-perturbative determinations of refs. , whose results are collected in table 1. Notice that at $`\beta =6.0`$ we opted for the value of $`c_A`$ obtained in  (and recently confirmed in ), whose absolute value is significantly smaller than the one obtained in ref. . Had we used the value of $`c_A`$ obtained in ref. , we would have found values of the AWI quark masses, at $`\beta =6.0`$, smaller by approximately 7%.

In table 2, we show the results for the pseudoscalar and vector meson masses, in lattice units, obtained by fitting the corresponding correlation functions at zero spatial momentum in the time intervals indicated in table 1. For each value of the hopping parameter, we also present in table 2 the corresponding values of the VWI and AWI quark masses, defined in eqs. (5) and (7), renormalized in the $`\mathrm{RI}/\mathrm{MOM}`$ scheme at the scale $`\mu =3\mathrm{GeV}`$.

In order to get the physical values of quark masses, we follow the usual procedure  and fit the pseudoscalar meson masses to the following form

$$(am_P)^2=Q_1\left(am_1^{(\mathrm{VWI})}+am_2^{(\mathrm{VWI})}\right)+Q_2\left(am_1^{(\mathrm{VWI})}+am_2^{(\mathrm{VWI})}\right)^2,$$
(9)
and similarly for the AWI quark masses. The subscripts $`1`$ and $`2`$ denote the flavour of the two valence quarks in the meson. Since in this study we only considered mesons consisting of two degenerate quarks, we always have $`m_1=m_2`$ in the fits. For the same reason, we did not include a term proportional to $`(m_1m_2)^2`$ in eq. (9). The physical values of the average up/down and of the strange quark masses are then obtained by substituting the experimental pion and kaon masses on the l.h.s. of eq. (9) and the values of the lattice spacing listed in table 1. Notice that we do not distinguish the up from the down quark mass and, by using eq. (9), we can only determine the average of the two, i.e. $`m_{\mathrm{}}=(m_u+m_d)/2`$. The fit of the pseudoscalar meson masses to eq. (9), at $`\beta =6.2`$, and the resulting extrapolation to the physical values, is shown in figure 1, for both the VWI and AWI quark masses. As can be seen from the figure, the effect of including a quadratic term in the chiral extrapolation of quark masses is rather negligible. This is true for all values of the lattice spacing considered in this study.

The results for the strange and the average up/down quark masses, in the $`\mathrm{RI}/\mathrm{MOM}`$ scheme, at the renormalization scale $`\mu =3\mathrm{GeV}`$, are collected in table 3. Also shown in the table are the values of the ratio $`m_s/m_{\mathrm{}}`$ (which is a both scheme and scale independent quantity). Note that this quantity can be determined with much better accuracy than the mass $`m_{\mathrm{}}`$ itself, since statistical and systematic errors largely cancel in the ratio. We emphasize that the results for quark masses presented in table 3 are obtained in a completely non-perturbative way.

## 3 Conversion to the $`\overline{\mathrm{MS}}`$ scheme and extrapolation to the continuum limit

We now convert the $`\mathrm{RI}/\mathrm{MOM}`$ quark masses, obtained in the previous section, to the popular $`\overline{\mathrm{MS}}`$ scheme, in which the light quark masses are conventionally expressed at the scale $`\mu =2`$ GeV. That allows to compare our results to the results obtained by other lattice groups and to the ones obtained by using QCD sum rules. It is only at this stage of the calculation that we are forced to introduce perturbation theory. This is because the $`\overline{\mathrm{MS}}`$ scheme, being related to dimensional regularization, is defined in perturbation theory only.

The conversion factor from the $`\mathrm{RI}/\mathrm{MOM}`$ to the $`\overline{\mathrm{MS}}`$ scheme is conveniently calculated by introducing the renormalization group invariant mass, $`m^{\mathrm{RGI}}`$, defined by dividing out from the renormalized quark masses the perturbative scale dependence

$$m^{\mathrm{RGI}}=\frac{m^{\mathrm{RI}}(\mu )}{c^{\mathrm{RI}}(\mu )}=\frac{m^{\overline{\mathrm{MS}}}(\mu )}{c^{\overline{\mathrm{MS}}}(\mu )}.$$
(10)
The mass $`m^{\mathrm{RGI}}`$ is, by definition, both renormalization scale and renormalization scheme independent. The beta function of QCD and the quark mass anomalous dimension, entering the functions $`c(\mu )`$ in eq. (10), are known to 4-loop accuracy, in both the $`\mathrm{RI}/\mathrm{MOM}`$  and the $`\overline{\mathrm{MS}}`$ schemes -. From these papers, we extract

$$c^{\mathrm{RI}}(\mu )=\left(\frac{\alpha _s(\mu )}{\pi }\right)^{\frac{12}{25}}\left[1+2.34747\frac{\alpha _s(\mu )}{\pi }+12.0599\left(\frac{\alpha _s(\mu )}{\pi }\right)^2+84.4076\left(\frac{\alpha _s(\mu )}{\pi }\right)^3\right],$$
(11)
and

$$c^{\overline{\mathrm{MS}}}(\mu )=\left(\frac{\alpha _s(\mu )}{\pi }\right)^{\frac{12}{25}}\left[1+1.01413\frac{\alpha _s(\mu )}{\pi }+1.38921\left(\frac{\alpha _s(\mu )}{\pi }\right)^2+1.09054\left(\frac{\alpha _s(\mu )}{\pi }\right)^3\right],$$
(12)
by using $`n_F=4`$ as the number of active flavours in the range of scales between 2 and 3 GeV. Then, by using $`\alpha _s(m_Z)=0.118`$, we obtain the conversion factor

$$R^{(4)}\frac{m_q^{\overline{\mathrm{MS}}}(2\mathrm{GeV})}{m^{\mathrm{RI}}(3\mathrm{GeV})}=\frac{c^{\overline{\mathrm{MS}}}(2\mathrm{GeV})}{c^{\mathrm{RI}}(3\mathrm{GeV})}=0.918.$$
(13)
This result has a $`\mathrm{N}^3\mathrm{LO}`$ accuracy in continuum perturbation theory. Therefore, we expect the perturbative error in our determination of the $`\overline{\mathrm{MS}}`$ quark masses to be completely negligible.

We use eq. (13) to convert the results for the $`\mathrm{RI}/\mathrm{MOM}`$ quark masses, presented in table 3, to the $`\overline{\mathrm{MS}}`$ scheme. The resulting values of $`m_{\mathrm{}}^{\overline{\mathrm{MS}}}(2\mathrm{GeV})`$ and $`m_s^{\overline{\mathrm{MS}}}(2\mathrm{GeV})`$, as obtained from the VWI and AWI respectively, are shown in table 4. We do not report in the table the values of the ratio $`m_s/m_{\mathrm{}}`$, presented in table 3, since this quantity is scheme independent.

The last step of our calculation is the extrapolation to the continuum limit. As discussed in the previous section, we find it convenient for that purpose to fix the relative values of the lattice spacing by using eq. (8), and taking as input the central value $`a^1(\beta =6)=2.0\mathrm{GeV}`$. Since our results for quark masses are free of leading $`𝒪(a)`$-effects, the first term in the extrapolation is of $`𝒪(a^2)`$. For this reason, we extrapolate our data to the continuum limit linearly in $`a^2`$, and obtain the results

$`m_s^{\mathrm{VWI}}[\mathrm{MeV}]=(102\pm 2)(2.01\pm 0.25)a^2,`$
$`m_{\mathrm{}}^{\mathrm{VWI}}[\mathrm{MeV}]=(4.20\pm 0.10)(0.100\pm 0.013)a^2,`$ (14)
$`(m_s/m_{\mathrm{}})^{\mathrm{VWI}}=(24.25\pm 0.17)+(140\pm 25)a^2,`$
from the VWI, and

$`m_s^{\mathrm{AWI}}[\mathrm{MeV}]=(106\pm 2)(1.39\pm 0.28)a^2,`$
$`m_{\mathrm{}}^{\mathrm{AWI}}[\mathrm{MeV}]=(4.35\pm 0.11)(0.076\pm 0.014)a^2,`$ (15)
$`(m_s/m_{\mathrm{}})^{\mathrm{AWI}}=(24.32\pm 0.17)+(127\pm 25)a^2,`$
from the AWI, where the lattice spacing must be expressed in units of fm. The illustration of the extrapolation of the $`\overline{\mathrm{MS}}`$ strange quark mass to the continuum limit is shown in fig. 2.

A pleasant feature of our results is that the quark masses extrapolated to the continuum as obtained by using either the AWI or the VWI lead to fully consistent determinations. The $`𝒪(a)`$-improved quark masses, with the lattice spacings used in our simulations, are very close to the continuum limit. Although the central values for the quark masses obtained at $`\beta =6.45`$ are slightly smaller than the ones obtained at $`\beta =6.4`$, they are completely consistent within the statistical errors.

## 4 Systematic uncertainties and final results

The errors quoted with the continuum determination of quark masses, in eqs. (3) and (3), are statistical only. In this section, we discuss the systematic uncertainties and present our final results.

The main sources of systematic errors, present in our calculation, are discussed below.

$``$Additive renormalization of the VWI quark mass: when using the VWI method, the renormalized quark mass is obtained from the bare one by implementing both a multiplicative and an additive renormalization. The latter is defined by the critical value of the Wilson hopping parameter, $`\kappa _{cr}`$, which, in this study, has been determined from the vanishing of the two-point correlation function of the divergence of the axial current. An equivalent possibility to fix $`\kappa _{cr}`$ is to require the pseudoscalar meson mass squared to vanish in the chiral limit. The values of light quark masses are rather sensitive to the precise choice of $`\kappa _{cr}`$, and, when using this alternative determination, we obtain the results

$`m_s^{\mathrm{VWI}}[\mathrm{MeV}]=(110\pm 2)(1.17\pm 0.28)a^2,`$
$`m_{\mathrm{}}^{\mathrm{VWI}}[\mathrm{MeV}]=(4.57\pm 0.11)(0.052\pm 0.016)a^2,`$ (16)
$`(m_s/m_{\mathrm{}})^{\mathrm{VWI}}=(23.94\pm 0.18)+(20\pm 27)a^2,`$
to be compared with those given in eq. (3). By combining the two sets of determinations in eqs. (3) and (4), we get the following estimates of the VWI quark masses

$`m_s^{\mathrm{VWI}}=(106\pm 2\pm 4)\mathrm{MeV},`$
$`m_{\mathrm{}}^{\mathrm{VWI}}=(4.38\pm 0.10\pm 0.18)\mathrm{MeV},`$ (17)
$`(m_s/m_{\mathrm{}})^{\mathrm{VWI}}=(24.10\pm 0.17\pm 0.15),`$
where the first error is statistical and the second represents the systematic uncertainty due to the spread of the two determinations of $`\kappa _{cr}`$.

The VWI results in eq. (4) are in perfect agreement with those obtained in eq. (3) by using the AWI method, the differences being smaller than 1%. The systematic uncertainty, however, is larger in the VWI case. For this reason, we will quote as our final central values of the results obtained from the AWI method, whereas the difference between the two methods will be included in the systematic error.

$``$Determination of the lattice spacing: our estimate of the lattice spacing has been performed by using eq. (8), in which the main source of uncertainty comes from the input value of the lattice scale at the reference point $`\beta =6`$. Our choice $`a^1(\beta =6)=2.0(1)\mathrm{GeV}`$ covers, in a rather conservative way, determinations of the scale based on different physical quantities, like $`f_\pi `$, $`m_\rho `$, $`m_K^{}`$, $`r_0`$, etc., which are not expected to produce the same estimate of the scale in the quenched approximation. In order to evaluate the effect of this uncertainty on the determination of the light quark masses, we have repeated the analysis by using for $`a^1(\beta =6)`$ the values 1.9 and 2.1 GeV respectively, and compared the results with those given in eqs. (3) and (3) obtained by using the central value $`a^1(\beta =6)=2.0\mathrm{GeV}`$. In this way, we find that the quark masses vary by approximately 5% (the mass increases as the lattice spacing increases).

$``$Renormalization constants: a reasonable estimate of the systematic uncertainty involved in the non-perturbative $`\mathrm{RI}/\mathrm{MOM}`$ calculation of the renormalization constants, $`Z_A`$, $`Z_S`$ and $`Z_P`$, can be obtained by comparing the results with those obtained by using the chiral Ward identity method . In the latter case, only scale independent quantities can be computed. In particular, the values of the renormalization constant $`Z_A`$ obtained from the two methods are in perfect agreement within the statistical errors, while a systematic difference of the order of 5% is observed in the case of the ratio $`Z_S/Z_P`$. We include this difference in the systematic uncertainty.

$``$Finite volume effects: the spatial extension of the lattices considered in this study is of the order of 1.6-1.7 fm, which is expected to be large enough for finite volume effects to be well under control. In order to verify this statement, we have performed an independent simulation, at $`\beta =6.0`$, on the volume $`24^3\times 64`$, which corresponds (in physical units) to a spatial extension of 2.3 fm. The results for the VWI and AWI strange quark masses, in the $`\overline{\mathrm{MS}}`$ scheme at the scale $`\mu =2\mathrm{GeV}`$, which are obtained from the simulation on the larger lattice, are $`m_s^{\mathrm{VWI}}=85(1)\mathrm{MeV}`$ and $`m_s^{\mathrm{AWI}}=93(3)\mathrm{MeV}`$ respectively. These results should be compared with the values $`m_s^{\mathrm{VWI}}=82(1)\mathrm{MeV}`$ and $`m_s^{\mathrm{AWI}}=92(1)\mathrm{MeV}`$ quoted in table 4. From this comparison, we get an estimate of finite volume effects which is of the order of 2%, which we also account for in our final systematic error.

$``$Continuum extrapolation: the extrapolation of quark masses to the continuum limit has been performed by considering only the effect of a linear term in $`a^2`$. One may wonder, however, whether higher order discretization effects are indeed negligible, particularly for the results obtained on the coarsest lattice. For this reason, we have also performed the continuum extrapolation without including the point at $`\beta =6.0`$ in the fit. In this way, we find that the results for the strange and the average up/down quark masses decrease by 1% and 2% respectively.

$``$Perturbative matching: the conversion factor $`R`$ of eq. (13), which translates the $`\mathrm{RI}/\mathrm{MOM}`$ quark masses to the masses in the $`\overline{\mathrm{MS}}`$ scheme, has been computed by using $`\alpha _s(m_Z)=0.118`$ and the number $`n_F=4`$ of active flavours in the range of scales between 2 and 3 GeV. By working in the quenched approximation, however, we could have equally computed this factor by using $`n_F=0`$ and $`\alpha _s`$ from $`\mathrm{\Lambda }_{\mathrm{QCD}}^{n_F=0}0.250`$ GeV. By proceeding in this way, we would have obtained $`R^{(0)}=0.944`$, instead of the result $`R^{(4)}=0.918`$ given in eq. (13). The difference between the two determinations, which is of the order of 3%, represents an intrinsic ambiguity in the quenched approximation. We include this difference in the final evaluation of the systematic uncertainty.

From the continuum results shown in eq. (3), and by adding in quadrature the systematic uncertainties discussed above, we finally obtain our best estimates of light quark masses,

$`m_s^{\overline{\mathrm{MS}}}(2\mathrm{GeV})=\left(106\pm 2\pm 8\right)\mathrm{MeV},`$ (18)
(19)
$`m_{\mathrm{}}^{\overline{\mathrm{MS}}}(2\mathrm{GeV})=\left(4.4\pm 0.1\pm 0.4\right)\mathrm{MeV}`$ (20)
and

$$m_s/m_{\mathrm{}}=24.3\pm 0.2\pm 0.6,$$
(21)
which have been also quoted in the abstract of this paper. Notice that most of statistical and systematic uncertainties cancel in the ratio $`m_s/m_{\mathrm{}}`$.

Our results in eqs. (18) and (21) are in agreement with the recent extensive lattice QCD calculations of the quark masses performed by using Wilson fermions in the quenched approximation -. They are also in good agreement with the current lattice world averages, presented in the reviews  and in the 2002 Review of Particle Physics . The main feature of the present study is, in our opinion, the special attention dedicated to the reduction and control of the systematic uncertainties within the quenched approximation.

## Acknowledgments

We thank Philippe Boucaud and Guido Martinelli for interesting discussions on the subject of this paper. Work partially supported by the European Community’s Human Potential Programme under HPRN-CT-2000-00145 Hadrons/Lattice QCD.