text:
# Lattice calculation of 1⁻⁺ hybrid mesons with improved Kogut-Susskind fermions

## I INTRODUCTION

The fact that gluons carry color charge suggests that they, like quarks, could be “valence” constituents of hadrons. In other words, we expect that the spectrum of QCD should contain glueballs and hybrids, or particles with both quarks and gluons as valence constituents. Hybrid mesons can have exotic quantum numbers, or $`J^{PC}`$ combinations not possible with a quark-antiquark state. However, a state with exotic quantum numbers is not necessarily a hybrid — it could be a $`\overline{q}\overline{q}qq`$ state, realized either as a single “bag” containing four quarks or as a “molecule” made of two $`\overline{q}q`$ mesons. Experimental evidence suggests the existence of one or more mesons with exotic quantum numbers $`J^{PC}=1^+`$, namely the $`\pi _1(1400)`$ pi1400  and the $`\pi _1(1600)`$ pi1600 . Analytic and numerical methods to predict the mass of light hybrid meson states include flux tube models Barnes:1995hc , the bag model Flensburg ; Barnes\_bag ; Chanowitz ; Barnes\_bag2 ; Barnes\_bag3 , QCD spectral sum rules Chetyrkin:2000tj ; Latorre:1985tg ; Narison:1999hg , relativistic WKB calculations Cornwall , and lattice QCD. Several lattice studies Lacock:1996ny ; Bernard:1997ib ; ZSU  have used quenched Wilson or quenched Wilson-clover fermions to calculate the masses of exotic hybrid states, although with quark masses much larger than the physical $`u`$ and $`d`$ quark masses. Lacock and Schilling have done a calculation in two flavor QCD, again with fairly heavy quarks Lacock\_dyn .

Here we report results of a lattice calculation of the mass of a $`1^+`$ hybrid meson using improved Kogut-Susskind quarks. The use of Kogut-Susskind quarks allows us to work at valence quark masses much smaller than were used in previous lattice calculations. In addition, the “$`a_{\mathrm{tad}}^2`$” action that we use has leading lattice spacing errors of order $`a^2g^2`$, while the clover-Wilson action has errors of order $`a^2`$. Our mass estimates in the quenched approximation are consistent with earlier Wilson quark results, but extrapolation to the physical light valence quark masses is under much better control. Preliminary results of this calculation were reported in Ref. LAT02\_hybrids .

We have also calculated hybrid meson propagators including the effects of three flavors of dynamical quarks, with light sea quark masses down to $`0.4`$ times the strange quark mass. We find that extracting mass estimates from the propagators in full QCD is difficult, and we argue that this difficulty is due to mixing of the hybrid meson with two meson states — the states into which it might decay.

## II $`1^+`$ HYBRID MESON OPERATOR

We can construct a $`1^+`$ hybrid meson operator as the cross product of a color octet $`1^{}`$ quark-antiquark ($`\rho `$ meson) operator and the chromomagnetic field, which has $`J^{PC}=1^+`$: $`\rho \times B`$ Bernard:1997ib . With staggered quarks we have several choices of rho meson operators, but it is convenient to choose the taste 111 We use the term “taste” to refer to the four types of quarks that are naturally present in the Kogut-Susskind formulation, while “flavor” can also distinguish quarks with an additional externally imposed label. For example, a meson with a source operator $`\overline{\psi }\gamma _5\mathrm{𝟏}\psi `$ but with disconnected diagrams not included would be a taste singlet but flavor non-singlet, and would be a pion in the continuum limit. singlet $`\rho _s`$, with the spin $``$ taste structure $`\gamma _i\mathrm{𝟏}`$.

$`1_i^+`$ $`=`$ $`ϵ_{ijk}\overline{\psi }^a\gamma _j\mathrm{𝟏}\psi ^bB_k^{ab}`$ (1)
$`=`$ $`2\overline{\psi }^a\gamma _j\mathrm{𝟏}\psi ^bF_{ij}^{ab},`$
where $`i`$, $`j`$ and $`k`$ are spatial indices and $`a`$ and $`b`$ are color indices. Each spin component of the $`1^+`$ includes two terms, for example:

$$1_x^+=\rho _yB_z\rho _zB_y,$$
(2)
so if we had chosen a spin $``$ taste structure like $`\gamma _i\gamma _i`$ the two components of $`1_x^+`$ would have different tastes.

The Kogut-Susskind $`\rho _s`$ meson operator, with spin aligned in the $`k`$ direction is $`\overline{\chi }\eta _kD_k\chi ,`$ Golterman:1985dz  where $`\chi `$ and $`\overline{\chi }`$ are the quark and antiquark fields respectively. The covariant symmetric shift operator is given by

$$D_\mu q(x)=\frac{1}{2}\left[U_\mu ^{}(x\widehat{\mu })q(x\widehat{\mu })+U_\mu (x)q(x+\widehat{\mu })\right].$$
(3)
We compute the field strength at each lattice point using the four plaquettes in each plane that have corners at this point, as described in Ref. Bernard:1997ib . In computing the field strength, we use links that have been smoothed with 32 iterations of APE smearing in the spatial directions only with relative weight of the staples set to 0.25 APE . This smearing removes short wavelength fluctuations in the gluon field, and reduces the noise in the hybrid propagator. (The smeared links are only used in constructing $`F_{\mu \nu }`$; the propagators are computed using the original links.)

Our zero momentum hybrid source and sink wave functions are constructed in Coulomb gauge and consist of a product of quark and antiquark fields with phases and offsets appropriate to a color octet $`\rho _s`$, as described above, and multiplied by the smeared field strength symmetrized with respect to the positions of the quark or antiquark to form the required $`C`$ even combination, as illustrated in Fig. 1:

$$\overline{\chi }ϵ_{ijk}\left(\eta _iD_iB_j+B_j\eta _iD_i\right)\chi $$
(4)
The operator is summed over all spatial sites and a trace is taken over the color indices.

The algorithm for constructing the meson propagator starts in Coulomb gauge with a quark “wall source”, consisting of a unit color vector field in a spatially constant direction, and applies the hybrid meson operator to form a source for the antiquark propagator. The calculation of the meson propagator is completed by acting upon the resulting antiquark propagator at an arbitrary time slice by the same hybrid operator and joining the resulting color vector field with the quark field propagated from the same wall source, summing over all sink spatial sites and color indices. The whole process is repeated, summing over the three wall source colors.

## III SIMULATION AND MEASUREMENT

We measured the connected correlator of the $`1^+`$ hybrid state on three sets of $`28^3\times 96`$ lattices generated with the “$`a_{\mathrm{tad}}^2`$” action Asqtad . To isolate the effects of dynamical quarks, we used matched quenched and full QCD lattices with $`10/g^2=8.40`$, $`m_{\mathrm{val}}a=0.016,0.04`$, for the quenched quarks, $`10/g^2=7.18`$ for lattices with three degenerate flavors of dynamical sea quarks at the strange quark mass ($`ma=0.031`$) and $`10/g^2=7.11`$ for lattices with $`m_{u,d}=0.4m_s`$ ($`ma=0.0124`$). These choices of $`10/g^2`$ give approximately the same lattice spacing ($`0.09`$ fm) in the three cases. The corresponding choices of quark mass allow simulation at roughly equivalent values of $`(m_{\mathrm{PS}}/m_V)^2`$, the square of the ratio of the pseudoscalar to vector meson masses. Table 1 summarizes the simulation parameters and fit results for the $`1^+`$ states, while Table 2 contains estimates for conventional hadron masses at these parameters.

The size of the datasets is comparable for quenched and full QCD runs. Successive full QCD lattices are separated by six molecular dynamics trajectories, with each trajectory one simulation time unit long. The full QCD lattices are not completely decorrelated but this autocorrelation has negligible effects on the hybrid mass fittings, since hybrid propagators have much larger statistical errors than, e.g. pion propagators. In particular, for the lightest sea and valence quark mass, $`am_q=0.0124`$, we calculated the normalized autocorrelations of the $`1^+`$ propagators separated by six simulation time units at each Euclidean time separation, or distance between the wall source and sink. For propagation distances zero through eight with the sample of 532 lattices we find 0.01, 0.13, -0.05, -0.00, 0.08, 0.01, -0.04, -0.16 and -0.08 respectively, instead of the uniformly positive autocorrelations that we would see if the propagators were systematically correlated from one stored lattice to the next. Although the statistical errors we quote come from the covariance matrix of the propagator, we have also performed a jackknife error analysis of each fitted mass and found jackknife error estimates to be consistent with errors from the covariance matrix. Varying the block size from 1 to 10 had no significant effect on the jackknife error.

In a separate study we have measured propagators of the pion, rho and nucleon. Statistical errors on these propagators are much smaller than for the $`1^+`$ propagator, so some effects of autocorrelations can be seen. For the nucleon at mass $`am_1=0.0124`$, which we use for comparison with the hybrid propagators, the data was grouped in blocks of four lattices, or 24 trajectories, before the covariance matrix was computed. Further blocking does not significantly increase the error bars. The fact that the nucleon mass fits have good $`\chi ^2`$ (in fact, better than the quenched nucleon fits) is also evidence that this blocking has removed most of the effects of the autocorrelations.

## IV RESULTS

We fit the measured correlators to the sum of oscillating and normal exponentials:

$$C(t)=A_1e^{M_{1^+}t}+A_2(1)^te^{m_2t}+A_3(1)^te^{m_3t},$$
(5)
where $`M_{1^+}`$ is the hybrid meson mass of interest and $`m_2`$ and $`m_3`$ are masses of non-exotic parity partner states which have oscillating correlators in the Kogut-Susskind formulation. In our case the oscillating parity partner is a $`1^{++}`$ ($`a_1`$) state, which is lighter than the $`1^+`$ hybrid, and the oscillating component dominates the correlator at large times. It is therefore essential to include the oscillating state(s) in our fits. We performed both four and five parameter fits. For the four parameter fits, we fix $`A_3=m_3=0`$, meaning that we include one state of each parity. For the five parameter fits we fix $`m_2`$ to an $`a_1`$ meson mass determined from propagators with a standard $`\overline{q}q`$ source operator, and fit for $`A_2`$, $`m_3`$ and $`A_3`$. We varied the range of the fit and tried to choose values for $`M_{1^+}`$ corresponding to high-confidence fits that were insensitive to $`D_{\mathrm{max}}`$ and $`D_{\mathrm{min}}`$, the limits of the fit range.

For the quenched lattices we were able to fit the propagators with reasonable confidence levels (25-50%) for valence quark masses $`ma=0.016`$ and $`ma=0.040`$. Figure 2 shows the measured propagator for $`ma=0.016`$. Note the oscillating component due to parity partner states. As expected, the oscillating component dominates at large distance, since the parity partner has lower mass than the $`1^+`$. Figure 3 shows mass fits for the quenched lattices for $`ma=0.040`$ and $`D_{\mathrm{max}}=15`$, with both the two particle (four parameter) and three particle (five parameter) fits. In the mass fit plots, we have included the small confidence level fits to illustrate how adjusting the fit range produces more optimal fits. Figure 4 shows the same plot for $`ma=0.016`$. In both plots the three particle fits exhibit a plateau with relatively small error bars ($`<1\%`$), demonstrating the stability of the result with respect to variations in the fit range. For the four-parameter fits, there is a slight oscillation of fitted values about the same plateau. Furthermore, the range of fits with high confidence level and relatively small errors is reduced. From plots like these, we picked a “best fit”, a value that met some balance of the following criteria: insensitivity to fit range, high confidence level, reasonable statistical errors. We can see that one might reasonably choose any one of several points as a “best fit”, and the range of resulting $`M_{1^+}`$ values is the basis of our estimate of the systematic error coming from the presence of higher mass states in the propagators. In all of these fit summary figures we include unused fits, that do not meet these criteria, say, because of low confidence level, to help illustrate how we selected the optimal fits.

For lattices with three degenerate sea quarks at $`m_s`$, we were also able to extract a value for $`M_{1^+}`$ in reasonable agreement with the quenched result. Four and five parameter fits are shown in Fig. 5. The fits exhibit larger statistical errors than the quenched lattice fits, and a slight dependence on range. The mass estimate in Table 1 reflects this with significantly larger statistical and systematic error bars than in the quenched case.

The lattices with $`m_{u,d}=0.4m_s`$ proved more interesting and difficult. The $`1^+1^{++}`$ propagator for valence mass $`am_q=0.0124`$ for this ensemble is shown in Fig. 6. Fits to the $`1^+`$ mass for both valence masses are illustrated in Figs. 7 and 8. The fitted mass agrees with those of the quenched and three-flavor results within two standard deviations, but with larger systematic errors, estimated from the dependence on fit range.

In the case of the light valence quark ($`ma=0.0124`$), we were unable to say much about the $`1^+`$ hybrid mass with any confidence. It is apparent from visual examination of the propagator (Fig. 6) that there is a lessening of the overall slope, suggesting that the non-oscillating piece may not be consistent with a single exponential. Indeed, the fits were very range dependent. Together these factors indicate the presence of lighter $`1^+`$ states, likely to be the states of two mesons into which the hybrid can decay. However, with the statistics available to us, we are unable to get convincing plateaus in the fits with more than one exponential in the $`1^+`$ channel.

We performed a linear extrapolation in quark mass of the quenched results to the physical value of $`(m_{\mathrm{PS}}/m_V)^2`$. Because the calculations at the two quark masses were done on the same set of quenched configurations, they are highly correlated, and a single elimination jackknife method was used to estimate the statistical error of the extrapolation.

## V DISCUSSION AND CONCLUSIONS

There are several sources of systematic error to be estimated. The largest of these, namely use of the quenched approximation, is inextricably mixed with the problem of determining the overall scale, or lattice spacing, so we will discuss these issues together.

The first source of systematic error is due to the possibility of mixing of higher mass states in the $`1^+`$ propagators. As described above we estimate this by looking at the mass range one might get by a reasonable variation of the fitting parameters.

We also have effects of finite lattice spacing. We obtained these results on lattices with $`a0.09`$ fm. For the conventional hadrons, we have masses at both $`a0.13`$ fm and $`a0.09`$ fm (Figs. 9 and 10). Since errors with this action are expected to be order $`a^2g^2`$, and the finer lattice spacing is about $`1/\sqrt{2}`$ times the coarser lattice spacing, we expect that the difference between $`a0.09`$ fm and $`a=0`$ masses is comparable to or slightly smaller than the difference between $`a0.13`$ and $`0.09`$ fm. For the quenched $`m_\rho /\sqrt{\sigma }`$ and $`m_N/\sqrt{\sigma }`$ we see differences as large as 3% between the two lattice spacings, and a difference of about 2% in the ratio $`m_N/m_\rho `$ at the light quark mass. Differences are smaller at the heavier mass — less than 1% in the nucleon to rho mass ratio. Therefore we expect effects of finite lattice spacing on our results based on hadron mass ratios to be around 1% for strange quarks, and we will use an estimate of 3% for light quarks.

The finite size of the entire lattice also introduces systematic error. The $`28^3\times 96`$ lattice corresponds to a box $`(2.5\mathrm{fm})^3\times 8.6`$ fm. In one case, three flavor QCD with light quark mass about 0.2 times the strange quark mass with a lattice spacing of 0.13 fm, we have calculated light hadron masses both in a $`2.5`$ fm box and on a larger $`3.6`$ fm spatial lattice. The $`\rho `$, $`\varphi `$ and nucleon masses decrease by a barely significant 0.9(7)%, 0.25(25)% and 0.9(6)% respectively as the lattice size goes from 2.5 to 3.6 fm. Since these effects are expected to fall exponentially with lattice size, we can simply take these numbers as an estimate of the effect of the 2.5 fm box size on the light hadron masses. However, hybrids are expected to be rather extended objects and may feel the influence of a finite lattice more than smaller particles, so we will use an estimate of 2% for this systematic error.

The largest systematic errors come from use of the quenched approximation, from the choice of quantity used to set the lattice scale, and the necessity for an extrapolation to the physical value of the valence quark mass. These effects are interrelated and so must be discussed together.

The hybrid mass estimates obtained above are in units of the inverse lattice spacing $`a^1`$, so to convert these to physical units we need to know $`a`$. The lattice spacing is determined by calculating some quantity that is known from experiment. In other words, the simulation actually produces the ratio of the hybrid mass to some other dimensionful quantity. In a simulation with sea quark masses at their physical values, the choice of quantity to fix the lattice spacing would be just a question of convenience. However, in the quenched approximation, we will not get the real world values for ratios of masses, so there is an important choice to be made. Because it is easily measured, and because it does not require an extrapolation in valence quark masses, the static quark potential is often used to determine the lattice spacing. In particular, we may use the string tension, $`\sqrt{\sigma }440\mathrm{MeV}`$, the coefficient of the linear term in $`V(r)`$. We might also use $`r_00.50`$ fm or $`r_10.34`$ fm, which are defined by $`r_x^2F(r_x)=1.65`$ or $`1.00`$ respectively. However, the shape of the static quark potential in quenched QCD differs from the shape with three dynamical flavors MILC\_spectrum . Hybrid mesons are expected to be large hadrons where the quarks are more likely to be in the linear part of the static quark potential, where $`\sigma `$ is defined, rather than the region of crossover between Coulombic and linear behavior, where $`r_0`$ and $`r_1`$ are defined. This suggests that plotting results in units of the string tension might minimize (although by no means eliminate!) effects of quenching. This expectation is borne out by calculations of the conventional hadron spectrum with this same improved action, where using $`\sigma `$ to define the lattice spacing produces better agreement of the quenched and three flavor results than using $`r_1`$ MILC\_spectrum . Figures 9 and 10 illustrate this with rho and nucleon masses plotted in units of $`r_1^1`$ and $`\sqrt{\sigma }`$ respectively. Since one of our important goals is to compare quenched and three flavor results, we therefore plot our results in units of the string tension. We also wish to compare our results with earlier results, and for this purpose the string tension in other published simulations is either available or can be reasonably estimated. In Fig. 11 we summarize our results along with the results of previous Wilson quark studies by the MILC collaboration Bernard:1997ib , the UKQCD collaboration Lacock:1996ny , the SESAM collaboration Lacock\_dyn , as well as recent results from the Zhongshan University group ZSU  using Wilson quarks on an anisotropic lattice. We use the string tension $`\sigma `$ to establish the lattice length scale and plot $`M_{1^+}/\sqrt{\sigma }`$. Our results are consistent with the earlier results at heavier quark masses.

To compare with experiment, we need to convert $`M_H/\sqrt{\sigma }`$ to physical units. Unfortunately, although phenomenological estimates are available, the string tension is not a parameter that is well known from experiment. The obvious workaround is to determine the string tension from the lattice results for $`m_\rho /\sqrt{\sigma }`$ etc., which in the end means that we are using the light hadron spectrum to set the length scale. Since ratios of quenched hadron masses are not quite those of the real world, we will get different estimates of the length scale depending on which hadron we choose. For the $`s\overline{s}`$ hybrid, the most reasonable choice for setting the length scale is a hadron with valence quark masses at the same value — the $`\varphi `$ meson or $`\mathrm{\Omega }^{}`$ baryon, which means that we are essentially quoting $`M_H/M_\varphi `$ or $`M_H/M_\mathrm{\Omega }^{}`$ with the quenched $`\varphi `$ mass and $`\mathrm{\Omega }^{}`$ masses defined to be 1020 MeV and 1672 MeV. Estimating the masses of the conventional hadrons on our quenched lattices from a linear extrapolation of results at $`am_q=0.015`$ and $`0.030`$, and setting the quenched string tension from the $`\varphi `$ or $`\mathrm{\Omega }^{}`$ gives $`\sqrt{\sigma }=436(4)`$ or $`437(9)`$ MeV respectively. (This remarkable agreement is surely coincidence, since other hadron mass ratios on these lattices differ by much larger amounts from the real world.) To estimate the light quark hybrid mass in MeV, we might use these estimates of $`\sqrt{\sigma }`$, or equally well argue that we should use light quark hadrons for comparison. Using the linearly extrapolated or interpolated $`\rho `$, $`K^{}`$, $`N`$ or $`\mathrm{\Delta }`$ masses to set the scale gives quenched $`\sqrt{\sigma }`$ of 389(5), 410(4), 380(5) or 400(21) MeV respectively, showing statistical errors only. These estimates are in reasonable agreement with phenomenological estimates from potential models on charmonium and bottomonium spectroscopy; for example $`\sqrt{\sigma }=384`$ MeV or $`427`$ MeV in Refs. RICHARDSON  and EICHTEN  respectively. Thus in estimating light quark hybrid masses in MeV we might consider a range of possible values for the quenched $`\sqrt{\sigma }`$ from around 380 to 440 MeV.

We begin with estimates for the $`s\overline{s}`$ hybrid masses. As mentioned above, it seems most consistent to use masses of hadrons made from strange quarks to set the lattice spacing in this case. If we use the $`\varphi `$ meson to set the length scale, using the results in Tables 1 and 2 we find, with statistical error only

$$M_{H,s\overline{s}}=1020\mathrm{MeV}\left(\frac{1.062(12)}{0.523(3)}\right)=2071(26)\mathrm{MeV}.$$
(6)
Systematic errors include fit choice, nonzero lattice spacing, finite spatial size, and effects of quenching. The first three have been discussed above. Effects of quenching can be estimated in part from the variation of our mass estimates among different ways of fixing the lattice scale, and in part from differences of other hadronic ratios between full and quenched QCD, as for example in Fig. 10. Here we include what we expect is a fairly conservative 5% error for this effect, thus estimating

$`M_{H,s\overline{s}}`$ $`=`$ $`2071(26)(39)(1\%)(2\%)(5\%)`$ (7)
$`=`$ $`2071(120)\mathrm{MeV},`$ (8)
where the errors are statistical, fit choice, lattice spacing, box size and quenching respectively. A similar calculation using the $`\mathrm{\Omega }^{}`$ mass to set the scale gives

$`M_{H,s\overline{s}}`$ $`=`$ $`1672\left({\displaystyle \frac{1.062(12)}{0.855(17)}}\right)`$ (9)
$`=`$ $`2077(48)(39)(1\%)(2\%)(5\%)`$
$`=`$ $`2077(129)\mathrm{MeV}.`$
We might also use the mass of a fictional octet baryon made from three quarks with the mass of the strange quark, assigning it a mass of $`M_{sss}=m_N+\frac{3}{2}\left(m_\mathrm{\Xi }m_N\right)=1507\mathrm{MeV}`$:

$`M_{H,s\overline{s}}`$ $`=`$ $`1507\left({\displaystyle \frac{1.062(12)}{0.771(2)}}\right)`$ (10)
$`=`$ $`2075(24)(39)(1\%)(2\%)(5\%)`$
$`=`$ $`2075(119)\mathrm{MeV}.`$
These three estimates are in remarkably close, and doubtless partly fortuitous, agreement.

Repeating this calculation with the three flavor lattices with $`m_{u,d}=m_s`$ with the $`\varphi `$, $`\mathrm{\Omega }^{}`$ and $`sss`$ baryon setting the scale produces

$`M_{H,s\overline{s}}`$ $`=`$ $`1020\left({\displaystyle \frac{0.986(30)}{0.4778(9)}}\right)`$ (11)
$`=`$ $`2105(64)(64)(1\%)(2\%)(3\%)`$
$`=`$ $`2105(120)\mathrm{MeV}`$
$`M_{H,s\overline{s}}`$ $`=`$ $`1672\left({\displaystyle \frac{0.986(30)}{0.7659(24)}}\right)`$ (12)
$`=`$ $`2152(66)(66)(1\%)(2\%)(3\%)`$
$`=`$ $`2152(123)\mathrm{MeV}`$
$`M_{H,s\overline{s}}`$ $`=`$ $`1507\left({\displaystyle \frac{0.986(30)}{0.6991(10)}}\right)`$ (13)
$`=`$ $`2125(65)(65)(1\%)(2\%)(3\%)`$
$`=`$ $`2125(121)\mathrm{MeV}`$
respectively. Here we have assigned an error of 3% for the partial quenching, or the remaining extrapolation of the sea quark masses to their physical values. Finally, we made an estimate of the $`s\overline{s}`$ hybrid mass from the run with $`m_{u,d}=0.4m_s`$. Although the error on this estimate, mostly coming from the choice of fit range, is too large for it to be very useful, we include it for completeness.

$`M_{H,s\overline{s}}`$ $`=`$ $`1020\left({\displaystyle \frac{0.911(34)}{0.4792(16)}}\right)`$ (14)
$`=`$ $`1939(73)(213)(1\%)(2\%)(2\%)`$
$`=`$ $`1939(233)\mathrm{MeV}.`$
Since the sea quarks here are much lighter, we used 2% as our estimate of the systematic error from partial quenching in this number. We can summarize this with an estimate of $`2100\pm 120`$ MeV for the mass of the $`s\overline{s}`$ $`1^+`$ hybrid meson.

To estimate the mass of a light quark $`1^+`$ hybrid meson we use the jackknife extrapolation of the quenched results to $`\left(m_{PS}/m_V\right)^2=0.033`$, $`am_H=0.919(39)`$. If we use the $`\varphi `$ to set the scale, this would correspond to a mass of

$`M_{H,u\overline{u}}`$ $`=`$ $`1020\left({\displaystyle \frac{0.919(39)}{0.523(3)}}\right)`$ (15)
$`=`$ $`1792(77)(36)(3\%)(2\%)(5\%)`$
$`=`$ $`1792(139)\mathrm{MeV},`$
with similar results using the $`\mathrm{\Omega }^{}`$ or $`sss`$ baryon. However, if we were to use the smaller estimates of the string tension obtained from linear extrapolations of light quark hadron masses to the physical light quark mass, we would obtain smaller values around 1600 MeV. As discussed above we have assigned a larger 3% systematic error for the effect of nonzero lattice spacing. We have also assigned a larger 5% error from quenching and chiral extrapolation. One reason that a larger systematic error is required here is that we are estimating the lattice spacing in large part from hadrons made up of strange quarks. Our strange quark mass was fixed by tuning the pseudo-scalar to vector meson mass ratio, and would have come out slightly different if we had used some other quantity. The effect of uncertainty in fixing the strange quark mass mostly cancels from mass ratios of hadrons made up of strange quarks, such as $`M_{H,s\overline{s}}/M_\varphi `$, but will be present when quantities such as $`M_\varphi `$ are used in estimating the mass of light quark hadrons. More evidence that this larger systematic error is required is seen in the extrapolations of conventional hadron masses to the physical light quark mass. If a naive linear extrapolation is made, and the resulting masses used to set the scale for the $`1^+`$ hybrid mass, the close agreements of the scales from various conventional hadrons that we found when using hadrons made from strange quarks is no longer present, as seen in the string tension estimates above.

Given the systematic errors from quenching and chiral extrapolation, our estimate for the mass of the light quark $`1^+`$ meson is not inconsistent with the experimental candidate $`\pi _1(1600)`$. In Fig. 11 we include the $`1^+`$ experimental candidates $`\pi _1(1400)`$ and $`\pi _1(1600)`$ at the physical value of $`(m_{\mathrm{PS}}/m_V)^2=(m_\pi /m_\rho )^2=0.033`$. These particles are represented by vertical bars, offset slightly to the left or right for clarity, representing the range of values for the quenched string tension from 380 to 440 MeV.

The $`m_{u,d}=0.4m_s`$ data illustrates that dynamical quarks introduce new and significant processes that contribute to the $`1^+`$ propagator. On this same set of lattices, mass fits for stable hadrons, even with $`ma=0.0124`$ valence quarks, display plateaus as functions of minimum included distance, $`D_{\mathrm{min}}`$ with fixed maximum distance, $`D_{\mathrm{max}}`$. The plateaus are similar for quenched and full QCD. In contrast, for the $`1^+`$, the full QCD fits do not show even the shorter plateau found in the quenched fits. We illustrate this by comparing fit plots for quenched and full QCD hybrids and nucleons in Figure 12. Fit plots for nucleon and quenched hybrids show a plateau, indicating the propagator has a single exponential form in the region $`D_{\mathrm{min}}`$ to $`D_{\mathrm{max}}`$. The full QCD hybrid fit plot deviates from a plateau in a significant manner — at minimum distance five, in the range which we have generally used for our quoted mass estimate, the low mass full QCD fits drop to a smaller value. Though quenching often introduces a systematic effect in the mass, this propagator is different in a way that suggests mixing of more than one exponential, representing propagators of different states with $`J^{PC}=1^+`$. Our hybrid propagators with light, dynamical quarks show features that are not evident either in hybrid propagators with heavier or quenched quarks, or in stable hadron propagators even with light dynamical quarks.

Four-quark states, molecular states of two mesons, or two independent mesons can have $`J^{PC}=1^+`$ without the gluonic excitations. For example the combination of $`b_1+\pi `$ can give $`1^+`$ with $`I=1`$, and as the sum of these masses is less than the predicted mass of the lowest $`1^+`$ hybrid, we expect that dynamical quarks introduce the possibility of the hybrid decaying into this two-meson state. In fact, at the values of the quark masses that we used the $`1^+`$ energies found in our dynamical simulations, while similar to the quenched hybrid masses, are also very close to the expected decay channel masses. For the run with three degenerate sea quarks at $`m_s`$, our $`1^+`$ mass is $`am_H=0.97(3)(3)`$, very close to the sum of the “$`\pi `$”and “$`b_1`$” masses: $`am_\pi +am_{b1}=0.32+0.68=1.00`$. For the run with $`m_{u,d}=0.4m_s`$, we would expect decays into a pseudo-scalar $`K`$ and a P-wave strange meson - a $`K_1`$. Again, our estimated mass for the $`s\overline{s}`$ $`1^+`$, $`am_H=0.90(4)(10)`$, is close to the sum $`am_K+am_{K1}=0.27+0.63=0.90`$.

We now have ahead of us the task of understanding these contributions so that we can make useful predictions of the $`1^+`$ hybrid mass in the presence of dynamical quarks. It is clear from our results with dynamical quarks that it will not be sufficient to simply do the same analysis that was done on the quenched gauge configurations, simply replacing them with full QCD configurations. One obvious avenue that may shed some light is to measure cross-correlators between the $`\rho \times B`$ operator and the two-meson state, as was explored with Wilson quarks in Ref. Bernard:1997ib . A more detailed study along these lines in the static quark (heavy quark) limit has been done by the UKQCD collaboration UKQCDDECAY . It may also be useful to study the dependence of the exotic energy as a function of valence quark mass (possibly with fixed sea quark mass) to look for an avoided level crossing as the decay threshold is crossed, as was done for the non-exotic $`0^{++}`$ meson in Ref. MILC\_spectrum .

## ACKNOWLEDGEMENTS

Computations for this work were performed at the San Diego Supercomputer Center (SDSC), the Pittsburgh Supercomputer Center (PSC), Oak Ridge National Laboratory (ORNL) and the National Energy Resources Supercomputer Center (NERSC). This work was supported by the U.S. Department of Energy under contracts DOE – DE-FG02-91ER-40628, DOE – DE-FG02-91ER-40661, DOE – DE-FG02-97ER-41022 and DOE – DE-FG03-95ER-40906 and National Science Foundation grants NSF – PHY99-70701 and NSF – PHY00–98395.


format: text

text:
# Quantitative spectral analysis of early B-type supergiants in the Sculptor galaxy NGC 3001footnote 11footnote 1Based on observations obtained at the ESO Very Large Telescope

## 1 Introduction

The 8-10 meter class telescopes and their new generation instruments make it possible to extend the quantitative stellar spectroscopy beyond the Local Group. Early B-type supergiant stars are ideal targets for detailed spectroscopy even at low resolution (R$``$1000). Their blue spectra are rich in metal features which allows us the analysis of chemical species like C, N, O, Si and Mg. Although our knowledge of the evolution of massive stars still has open questions, most of the recent works indicate that the blue luminous supergiants do not show any contamination of their oxygen surface abundances during the early stages of their evolution, neither the O-types (Villamariz et al., 2002), nor the B-types (Smartt, Dufton, & Lennon, 1997; Monteverde, Herrero, & Lennon, 2000; Smartt et al., 2002), nor the A-types (Venn, 1995; Takeda & Takada-Hidei, 1998; Przybilla, 2002), which enables a direct comparison between the stellar oxygen abundances and the ones derived from H II regions. This has become extremely important, especially in the extragalactic field where oxygen is used as the primary metallicity indicator, due to the fact that at high metallicity (larger than approx. 0.5 solar) strong line methods must be used, for which the choice of the calibration strongly influences the derived abundances (Kewley & Dopita, 2002; Pilyugin, 2002). In addition to chemical abundance studies, blue luminous stars have strong radiatively driven mass outflows which can provide us with information on extragalatic distances by means of the Wind Momentum - Luminosity Relationship, WLR (Kudritzki & Puls, 2000, and references therein).

Recently, and within a wide program aimed at the spectroscopy study of luminous blue stars beyond the Local Group, first steps have been done for A-type supergiants in NGC 3621 (6.7 Mpc away, Bresolin et al., 2001). Quantitative spectroscopy has been shown to be possible for A-type supergiants (Bresolin et al., 2002a) and Wolf-Rayet stars (Bresolin et al., 2002b) in NGC 300, 2.02 Mpc away in the Sculptor group. Here we report the first quantitative analysis of B-type supergiants (hereafter B-Sg) out of the Local Group, presenting the detailed chemical pattern along with the stellar parameters and the wind properties. The technique will be applied in a forthcoming paper to a large set of early B-Sg located at several galactocentric distances in order to derive radial abundance gradients of the $`\alpha `$-elements. Combined with the results of a similar study of A-type supergiants it will provide a wealth of information on the chemical evolution of the host galaxy NGC 300.

## 2 Observations

The stars are part of a spectroscopic survey of photometrically selected blue luminous supergiants in the Sculptor galaxy NGC 300, obtained at the VLT with the FORS multiobject spectrograph, and described in detail by Bresolin et al. (2002a), which presents a spectral catalog of 70 luminous blue supergiants in the blue region ($``$ 4000 - 5000 Å). The selected stars are identified as B-12 and A-9 in that spectral catalog (see their Table 2 and finding charts). In September 2001 the spectra of the H$`\alpha `$ region were obtained in order to measure the mass-loss rates, which provide us with a complete coverage of the 3800 - 7200 Å wavelength range at R$``$1000 resolution. The reader is referred to Bresolin et al. (2002a) for a detailed description of the observations and reduction process, as well as for the photometry and the spectral classification of the stars.

## 3 Spectral analysis

The spectra of early B-Sg are dominated by the O II lines, followed by N II/N III, Si III/Si IV, C II/C III and Mg II, in addition to H and He I lines. At high resolution it is possible to detect some other metal lines of Al III, S II/S III and Fe III but, due to their intrinsic weakness, these lines do not have any influence in the analysis at low resolution and could hardly be used to fix the abundance of such elements. Fig. 1 shows the high resolution - high S/N ratio (R$``$15000, SNR$``$350) blue spectrum of the Galactic supergiant HD14956 (B1.5Ia), and the same spectrum degraded to the resolution of the NGC 300 data, R$``$1000 (labeled as #d in the figure). We have also included the identification of the more important lines. As can be seen, only a few strong lines remain isolated at that low resolution, therefore the analysis must be based on the comparison of the observed spectra to a set of model atmospheres that include a vast number of lines in the calculation of the emergent fluxes. We have taken into account more than two hundreds metal lines in the 3800 - 6000 Å wavelength range. It is important to include extense metal line lists because of the fact that some spectral features are formed by the contribution of several chemical species (e.g. the strong blend of O, N and C at $``$ 4650 Å). We have excluded some strong isolated lines because our atomic models do not consider the levels involved in these transitions. Nevertheless, these lines are isolated and have no influence on the results.

Even considering the noise effects in the lower resolution FORS spectra (displayed also in Fig. 1), strong metal features can still be detected and used for a detailed chemical abundance analysis. In particular a wealth of information can be extracted from the selected regions at 4070, 4320, 4420 (O II), 4550 - 4570 (Si III), 4600 - 4660 (O II, N II, N III and C III) and 5010 (N II and He I).

### 3.1 Atmosphere models

We use the newest version of the FASTWIND code (first presented by Santolaya-Rey, Puls, & Herrero, 1997) which solves the radiation transfer in a moving media by means of suitable approximations which simplify the numerical treatment of the problem but without affecting the physical significance of the results. The atmospheric structure is treated in a consistent way, assuming a $`\beta `$-velocity law in the wind, ensuring a smooth transition between the ”photosphere” and the ”wind”; the temperature structure is approximated by means of non-LTE Hopf functions carefully chosen to ensure the flux conservation better than 2 % at any depth point; rate equations are solved in the co-moving frame scheme, with the coupling between the radiation field and the rate equations solved using local ALOs (following Puls, 1991). This new version includes the effects of the line blanketing. The reader is referred to Puls et al. (2003, in preparation) for a detailed description. We have analysed two Galactic stars, 10 Lac (O9V) and HD209975 (O9.5Ib) in order to compare our results with the ones obtained with other codes. In the case of 10 Lac, our results agree with the recent ones by Herrero, Puls & Najarro (2002, see their comparison to the results by Hubeny et al. 1998). The derived parameters for HD209975 are consistent with the results by Villamariz et al. (2002) which used plane-parallel model with line blocking.

A model is prescribed by the effective temperature $`T_{eff}`$, the surface gravity log g, the stellar radius $`R_{}`$ (all these three quantities are defined at $`\tau _{Ross.}=2/3`$), the mass-loss rate $`\dot{M}`$, the wind terminal velocity $`v_{\mathrm{}}`$, the $`\beta `$ exponent of the wind velocity law, the He abundance $`Y_{He}`$, the microturbulent velocity $`v_{turb}`$ and, in the case of B-type stars, the Si abundance. The $`T_{eff}`$ is well determined from the Si III triplet and the blends of Si IV (with O II at 4090 Å and with O II/He I at 4120 Å), and the surface gravity from the Balmer hydrogen lines, provided that the mass-loss rate information is extracted from the H$`\alpha `$ profile. An important issue concerns the wind terminal velocity, that must be adopted from a spectral type - v empirical calibration (Haser, 1995; Kudritzki & Puls, 2000). The assumed terminal velocity affects the derived $`\dot{M}`$ and the log g. But, with the joined information from H$`\alpha `$ and H$`\beta `$, the mass-loss rate and v can be constrained to yield reasonable uncertainties in log g. The stellar radius is derived interactively from the absolute magnitude, deduced from the apparent magnitude after adopting a distance modulus ($`\mu =26.53`$, Freedman et al., 2001), and the model emergent flux (Kudritzki et al., 1999), which also provides the reddening by the comparison of the synthetic colors with the observed ones.

### 3.2 Results

Best-fitting models are displayed in Fig. 1 and the results summarized in Tab. 1. The derived $`\beta `$ values are consistent with those obtained by Kudritzki et al. (1999) for Galactic B-Sg, with lower values excluded by the arise of emission wings in the synthetic H$`\alpha `$ profiles. We estimated an uncertainty of $`\pm `$0.25 in $`\beta `$. In the case of B-12 only the higher Balmer lines have been considered in the surface gravity determination, as the cores of H$`\gamma `$ and H$`\beta `$ are particularly affected by the sky substraction. As it has been quoted, the O and N abundances are very well constrained because of the large number of features from these species. The presence of a lot of weak metal lines in the 4600 - 4700 Å wavelength range makes the selection of the continuum level in this area difficult, good S/N ratio is also needed to disentangle between a real feature and the noise effects. Final abundance uncertainties are estimated to be $`\pm 0.2`$ dex from model comparisons (see Fig. 2).

We define the mean metallicity as the sum of the $`\alpha `$-elements abundances, $`X_{Si}+X_{Mg}+X_O`$ and refer it to the Sun abundances by Grevesse & Sauval (1998); at the early stages of massive star evolution, the O surface abundance is not affected by the CNO cycle, which means that the abundance of the $`\alpha `$-elements is a direct measurement of the ZAMS metallicity of the star. The results for B-12, located close to the galactic center, resembles the abundance patterns of the early B-Sg in the solar neighborhood, having a solar metallicity within the uncertainties of the analysis. On the other hand A-9, in the outskirts of the galaxy, has clearly a lower metallicity, around 0.3 $`Z_{}`$. This is in agreement with the results for A-8, a B9-A0 supergiant close to A-9, by Bresolin et al. (2002a, see the Fig. 2). These authors find a mean metallicity of 0.2 $`Z_{}`$ for A-8. We must emphasize that both the model atmospheres and the metallicity indicators are different, but the results agree extremely well. The metallicity and the spatial location of both stars in NGC 300 points to a M33-like radial metallicity gradient. The CNO abundances indicate a different degree of chemical evolution, while B-12 displays a normal CNO spectrum, A-9 shows indications of strong N enrichment.

Synthetic magnitudes and colors (see Tab. 2) are consistent with almost no reddening for both stars, except the observed $`(VI)`$ for B-12 that seems to be anomalous, probably reflecting the presence of the H II region. Fig. 3 shows the location of the stars on the Hertzprung-Russel diagram, along with theoretical stellar tracks without rotation at solar metallicity from Schaller et al. (1992). We have also added the location of the Galactic stars 10 Lac, HD209975 and HD14956 as a reference.

Comparing the wind momentum of both NGC 300 stars with the results for Galactic supergiants (Fig. 4), B-12 agrees well with the results by Herrero et al. (2002) for O-type supergiants in the Galactic association Cyg OB2, as does HD14956. Note, however that the Herrero et al. (2002) stars are considerably hotter than the ones considered here. The wind momentum of A-9 is also compatible with the WLR of Galactic early B-Sg as derived by Kudritzki et al. (1999). With respect to this relationship, however, B-12 (being an early B-type supergiant as well) shows an enhanced wind momentum rate, which might be related to clumping effects in the wind that would produce an overestimation of the mass-loss rate. The failure of our models to reproduce the blue absorption of H$`\alpha `$ for B-12, in parallel with an H$`\gamma `$ core which is too strongly refilled might then be explained by this effect, at least in part, and not only by the rather problematic sky substraction outlined above. The location of A-9, compared to HD14956, reflects the lower metal content of the NGC 300 supergiant. It must be considered here that we have adopted the same $`v_{\mathrm{}}`$ for both stars, HD14956 and A-9, while a lower value for A-9 could be expected due to its lower metallicity (Kudritzki & Puls, 2000). In any case the effect of the lower wind terminal velocity would reduce even more the wind momentum of A-9, reinforcing the difference with respect to the Galactic B1.5Ia.

Recently Kudritzki, Bresolin & Przybilla (2003) have proposed a new extragalactic distance indicator, the ”Flux-weighted - Luminosity Relationship (FGLR)”. The results for both NGC 300 B-type supergiants, B-12 and A-9, follow this relationship, within the observed scatter (see the Fig. 2 of the latter reference).

We are gratefull to L. J. Corral for making us available the spectrum of HD14956. MAU thanks F. Najarro for providing the routines for the computation of the synthetic magnitudes. AH and MAU thank the Spanish MCyT for a support under proyect PNAYA2001-0436, partially funded with FEDER funds from the EU. WG gratefully acknowledges financial support for this work from the Chilean Center for Astrophysics FONDAP 15010003.


text:
# Discovery of a new radio galaxy within the error box of the unidentified 𝛾-ray source 3EG J1735-1500

## 1 Introduction

The quest for the identification of the $`\gamma `$-ray sources detected by EGRET instrument of the Compton Gamma-Ray Observatory is one of the most important challenges of high-energy astrophysics in recent years. Among the 271 sources included in the Third EGRET (3EG) catalog (Hartman et al. 1999) there are a few confirmed pulsars, 66 blazars, and a bunch of miscellaneous identifications including one radio galaxy: the nearby Centaurus A, a Fanaroff-Riley (FR) Type I active galaxy (Sreekumar et al. 1999).

Most identified blazars are strong flat-spectrum radio sources (e.g. only 4 out 46 high-probability blazar detections in Mattox et al. (2001) analysis have fluxes below 1 Jy) with jets pointing close to the line of sight. On the contrary, in the case of Centaurus A the viewing angle is rather large ($`70\mathrm{°}`$, Bailey et al. 1986). Very recently, Mukherjee et al. (2002) have suggested that the radio galaxy NGC 6251 could also be associated with an EGRET source. If this is confirmed it could have important consequences because the spatial density of FR I radio galaxies is far higher than that of blazars, and hence, despite that these objects are expected to be weaker $`\gamma `$-ray emitters, there could be many other unidentified sources associated with them.

In this paper we report the discovery of a new radio galaxy within the location error box of the $`\gamma `$-ray source 3EG J1735$``$1500. The object was found during a re-analysis of the main sources in the radio field around the position of the EGRET detection. This research is part of a systematic program to study potential low-energy counterparts of unidentified $`\gamma `$-ray sources. Previous results regarding compact radio sources were published by Torres et al. (2001). Here we present the results of a more detailed study and new observations that led to the identification of the new radio galaxy among the potential counterparts of 3EG J1735$``$1500. We also provide information on all compact radio sources within the EGRET location error box, including spectral index determinations when possible. We have also found that the source PMN J1738$``$1502 is probably a weak flat-spectrum quasar with mild radio flux, which enhances its probability of being the counterpart of 3EG J1735$``$1500 respect to the a priori probability estimates by Mattox et al. (2001).

## 2 Radio data analysis and results

We have used the NRAO VLA Sky Survey (NVSS) (Condon et al. 1998) to study the small-scale radio emission within the inner 95 % location probability contours of the $`\gamma `$-ray source 3EG J1735$``$1500. Twenty-three radio sources with flux density $`>10`$ mJy at 1.4 GHz have been found, as can be seen in Fig. 1, where we have numbered the sources with increasing galactic longitude. The main characteristics of these sources are listed in Table 1. In particular, we provide, from left to right, the identification number corresponding to Fig. 1, the galactic coordinates, flux density at 1.4 GHz, spectral index (defined in such a way that $`S_\nu \nu ^{+\alpha }`$) when observations at other frequencies are available (365 MHz from Douglas et al. 1996 and 4.8 GHz from Griffith & Wright 1993) and the name in other catalogs if any.

With the exceptions of sources no. 8 and 10, all entries in Table 1 correspond to compact objects, at least at the present angular resolution (43″). Source no. 8, which appears as extended in our map, is actually a composite of three different (and probably compact) weak sources listed as no. 7, 8, and 9 in Torres et al. (2001). Source no. 10, instead, is really extended. We shall name this source hereafter as J1737$``$15, according to its radio coordinates. Its morphology at 1.4 GHz (Fig. 2, upper panel) is typical of double-sided FR II radio galaxies. There is a nuclear component as well as two jets ending in large radio lobes. The angular dimensions of J1737$``$15 are $`0.04\mathrm{°}\times 0.11\mathrm{°}`$, with an integrated flux of $`55.6\pm 1.5`$ mJy at 1.4 GHz. The flux density of the northern component is $`S_{\mathrm{north}}=21.6\pm 0.5`$ mJy, whereas the southern one is $`S_{\mathrm{south}}=30.7\pm 1.3`$ mJy. The central component has a flux density of $`S_{\mathrm{core}}=3.80\pm 0.05`$ mJy, and an estimated J2000.0 ICRS position of $`\alpha =17^\mathrm{h}37^\mathrm{m}12.9^\mathrm{s}\pm 0.3^\mathrm{s}`$, $`\delta =15\mathrm{°}11\mathrm{}02\mathrm{}\pm 15\mathrm{}`$.

The source is clearly non-thermal since it is not detected at 4.8 GHz (Griffith & Wright 1993). From the 4.8-GHz survey sensitivity we can infer an average steep spectral index $`\alpha 0.9`$, but it is not possible to say, within the present resolution, what is the index at the core because of contamination from the radio lobes, which have surely steeper spectra. The general morphology clearly supports the hypothesis of a radio galaxy with jets nearly perpendicular to the line of sight. Should the association between the radio galaxy and the EGRET source be confirmed, J1737$``$15 would be the second EGRET source (apart Centaurus A) to be detected with a large inclination angle.

## 3 Optical observations

With the aim of finding the optical counterpart of the radio galaxy we have made $`VRI`$ deep photometric observations of J1737$``$15 on 2002 June 9 at Calar Alto (Almería, Spain) with the 2.2-m telescope of the Centro Astronómico Hispano-Alemán (CAHA). Images were obtained using the Ritchey-Chrétien focus and CAFOS, with a scale factor of $`0.53\mathrm{}`$ per pixel and a $`8.8\mathrm{}\times 8.8\mathrm{}`$ field of view. A series of 600 s images of our target, as well as 200 s images of standard stars from the SA104 field of Landolt (1992), were obtained through the $`V`$, $`R`$ and $`I`$ Johnson filters. The observations were reduced using standard procedures (bias and dark subtraction and flat-field correction) within the IRAF software package. We performed a detailed astrometric reduction of the images using 20 field stars present in the USNO-A2.0 catalog (Monet et al. 1999), with an estimated astrometric accuracy of $`0.3\mathrm{}`$.

We show in Fig. 2 the obtained image through the $`I`$ Johnson filter, together with the contours of the radio source from the NVSS. It is clear at first sight the presence of a slightly elongated (east-west) optical object close to the center of the radio source. The fitted J2000.0 ICRS position of this object is $`\alpha =17^\mathrm{h}37^\mathrm{m}12.744^\mathrm{s}\pm 0.021^\mathrm{s}`$, $`\delta =15\mathrm{°}11\mathrm{}01.14\mathrm{}\pm 0.30\mathrm{}`$, well within the previously obtained error box in position for the core of the radio source. We also obtained absolute photometry of this object, which is believed to be accurate only to $`\pm `$0.2 magnitudes, because the crowded field around our target prevented a good estimate of the background. The obtained magnitudes are $`V=18.3\pm 0.2`$, $`R=16.9\pm 0.2`$, and $`I=15.1\pm 0.2`$. Using the galactic extinction estimates from Schlegel et al. (1998), implemented within the NASA/IPAC Extragalactic Database (NED), we obtain a color excess of $`E(BV)=0.495`$ in the direction of J1737$``$15. However, as pointed out by Arce & Goodman (1999) the values provided by the Schlegel’s model could be overestimated by a factor 1.3–1.5 in regions of smooth extinction with $`E(BV)>0.15`$ mag, as it happens in our case. As a resonable value, we have considered that the NED value is overestimated by a factor $`1.3\pm 0.2`$, and used an extinction of $`E(BV)=0.38\pm 0.05`$. This provides the following extinctions: $`A_V=3.315E(BV)=1.3\pm 0.2`$, $`A_R=2.673E(BV)=1.0\pm 0.1`$ and $`A_I=1.940E(BV)=0.7\pm 0.1`$ (Schlegel et al. 1998). Using these values we obtain the following dereddened magnitudes: $`V=17.0\pm 0.3`$, $`R=15.9\pm 0.2`$, and $`I=14.4\pm 0.2`$. These values imply color indices of $`VR=1.1\pm 0.3`$ and $`RI=1.5\pm 0.3`$ (where the errors have been computed directly from the uncertainties in measured magnitudes and assumed $`E(BV)`$).

If we suppose that this optical object is a star, the ranges of obtained dereddened colors only allow possibilities such as M3III or M0I, according to Ducati et al. (2001). However, comparison of the obtained apparent dereddened visual magnitude with the expected absolute visual magnitude (Wainscoat et al. 1992) implies distances of 33 and 465 kpc, respectively. The second value is unrealistic, while the first one implies a galactic height ($`b=9.2\mathrm{°}`$) of $`5`$ kpc, in clear disagreement of what is typically found for M3III stars. Therefore, the dereddened colors are clearly incompatible with the object being a single galactic star located by chance in the same direction of the radio galaxy core.

On the other hand, since the source is elongated, a 2-D Gaussian fit can provide information about the ellipticity and position angle. We have performed several Gaussian fits using different fitting box sizes. In all cases the ellipticity of the source seems to be around 0.14, while the position angle changes from $`85\mathrm{°}`$ (counterclockwise from north) to $`75\mathrm{°}`$ when reducing the size of the fitting box. In any case, the object is extended practically in the east-west direction, i.e., perpendicular to the radio jets. All these facts suggest that our identified optical counterpart is a galaxy, at the center of which the jets visible at radio wavelengths are produced.

In addition to the optical counterpart we have looked for infrared and X-ray sources at the position of J1737$``$15. Filtered IRAS images (Wheelock et al. 1991) at 12 and 100 microns of the region show no particular infrared enhancement at the central position of the radio galaxy. The ROSAT All-Sky Survey (0.1–2.4 keV) shows three X-ray sources inside of the $`\gamma `$-ray 95% contour but none of them is close to J1737$``$15.

## 4 Discussion

Depending on the jet and ambient medium parameters, most double-sided radio sources have sizes below $`300`$ kpc (Begelman et al. 1984). In the case of the radio galaxy reported here, and using standard Friedmann-Robertson-Walker formulae, this size translates into a possible distance smaller than 350 Mpc. If 3EG J1735$``$1500 is indeed the result of $`\gamma `$-ray emission in J1737$``$15, the intrinsic luminosity at $`E>100`$ MeV, at the distance quoted, should be smaller than $`2\times 10^{44}`$ erg s-1, although mild beaming could reduce this figure. A correct determination of the distance requires a knowledge of the redshift that is still unknown. The information we provide here for the optical counterpart will be useful to guide future spectroscopic observations. As in the case of Centaurus A and 3EG J1621+8203, the luminosity needed to generate the $`\gamma `$-ray source observed is several orders of magnitude less than those required for typical EGRET blazars ($`10^{45}`$$`10^{48}`$ erg s-1).

The $`\gamma `$-ray variability history of 3EG J1735$``$1500 is shown in Fig. 3. The source was detected only in three EGRET viewing periods (VPs) –in the remaining VPs only upper limits to the flux were established– and appears to be variable. Torres et al.’s (2001) index $`I`$ gives $`I=8.86`$, which place this source as one of the most variable gamma-ray sources in the EGRET catalog. Tompkins’ (1999) index $`\tau `$ is $`\tau =1.09_{0.00}^{10.1}`$, which, although inconclusive due to its lower limit, the central value and upper limit of the index also place this source as one of the most variable ones in the Catalog. (For a detailed discussion on variability indices see Torres, Pessah and Romero 2001). The likely variability status of this source might imply that the $`\gamma `$-ray emission is produced in the central region and not in the extended lobes, if related with the discovered radio galaxy. The photon spectral index seems to be unusually steep with a value $`\mathrm{\Gamma }=3.24\pm 0.47`$, steeper than the average of the $`\gamma `$-ray blazars spectra, a property also shared by Centaurus A and 3EG J1621+8203.

Since 3EG J1735$``$1500 is at a low galactic latitude it could be possible that some galactic object not visible at radio wavelengths be responsible of the high-energy emission. However, the source is not coincident with early-type stars, OB associations, or galactic X-ray binaries (Romero et al. 1999).

The only other known potential counterpart in addition to the galaxy discussed here is the source PMN J1738$``$1502. It is a compact radio source with a total flux at 1.4 MHz of $`0.3`$ Jy (source no. 16 in Table 1). We have calculated for this object a flat spectral index $`\alpha =0.17`$, which makes of it a clear blazar candidate. Mattox et al. (2001) give an a priori probability of only 0.07 % for a physical association. The a posteriori probability (see Mattox et al.’s paper for details) is slightly higher but still very small: 0.35 %. The radio source is weaker than most high-probability $`\gamma `$-ray blazars, but similar to most plausible (i.e. with lower a priori probability) identifications. Consequently, this source should be considered a serious alternative counterpart of the EGRET source despite the low probabilities given by Mattox et al. (2001). The object is similar to the plausible identification proposed recently by Wallace et al. (2002) in the case of 3EG J2006$``$2321, although we notice that in the case of 3EG J1735$``$1500 the high-energy photon index is far softer ($`\mathrm{\Gamma }=3.24\pm 0.47`$ vs $`\mathrm{\Gamma }=2.47\pm 0.44`$). Since $`\gamma `$-ray blazars tend to be highly polarized and variable, we are planning optical polarization observations of this source in order to confirm its blazar nature.

One of the main problems to explain the synchrotron emission in radio galaxies is that this radiation is distributed rather uniformly along distances $`100`$ kpc. To explain how the continuous acceleration of electrons take place in these huge sources remains as one of the most difficult topics in jet theory. Recently, it has been suggested by Neronov et al. (2002) that in some cases the relativistic leptonic population can be locally created by very-high energy $`\gamma `$-rays produced at the central source and injected into the jet, where they produce pairs through interactions with the 2.7 K cosmic background radiation. In this scenario lower energy (GeV) photons would escape from the radio galaxy and could even traverse the diffuse infrared background without further interaction reaching the Earth. The $`\gamma `$-rays could be produced in the central engine by photon-pion process involving disk photons (Neronov et al. 2002). Future multifrequency observations of J1737$``$15 could help to test the proposal of an association with 3EG J1735$``$1500. If the high-energy emission is confirmed, this radio galaxy could become an important natural laboratory to test the theories of non-thermal emission in extragalactic jet sources.

## 5 Conclusions

We have discovered a new radio galaxy inside the error box of the $`\gamma `$-ray source 3EG J1735$``$1500. It is a double-sided FR II type of object. Its linear angular size implies that it is relatively nearby. We have also identified the host optical galaxy, which appears elongated in a direction perpendicular to the radio structure. Future spectroscopic observations will help to fix the distance to this object. If it turns out to be closer that 100 Mpc, the radio galaxy could be a potential acceleration site of ultra-high energy cosmic rays, avoiding the so-called Greisen-Zatsepin-Kuz’min (GZK) cut-off (see Rachen & Biermann 1993).

We have also study the other important likely potential counterpart of 3EG J1735$``$1500. We have found that it is a flat-spectrum compact radio source, possibly a blazar, and hence perhaps a $`\gamma `$-ray emitting object. The remaining sources within the location confidence contours of the EGRET detection seem to be uninteresting from the point of view of the high-energy emission. Forthcoming observational studies of these two objects will shed light on the real nature of 3EG J1735$``$1500.

This work was partially supported by the Argentinian agencies CONICET (PIP 0430/98) and ANPCT (PICT 03-04881), as well as by the DGI of the Ministerio de Ciencia y Tecnología (Spain) under grant AYA2001-3092 and the European Regional Development Fund (ERDF/FEDER). Additional support was provided by Fundación Antorchas. The work of D.F.T. was performed under the auspices of the U.S. Department of Energy, National Nuclear Security Administration, by University of California Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48. D.F.T. is Lawrence Fellow in Astrophysics. M. R. is supported by a fellowship from CIRIT (Generalitat de Catalunya, ref. 1999 FI 00199). We thank the German- Spanish Astronomical Centre, Calar Alto, operated by the Max-Planck-Institut für Astronomie (Heidelberg) jointly with the Spanish National Commission for Astronomy, for carrying out the observations under the Director’s Discretionary Time Service.


text:
# Numerical simulation evidence of dynamical transverse Meissner effect and moving Bose glass phase

(March 15, 2024)

## Abstract

We present 3D numerical simulation results of moving vortex lattices in presence of 1D correlated disorder at zero temperature. Our results with field tilting confirm the theoritical predictions of a moving Bose glass phase, characterized by transverse pinning and dynamical transverse Meissner effect, the moving flux lines being localized along the correlated disorder direction. Beyond a critical transverse field, vortex lines exhibit along all their length a ”kink” structure resulting from an effective static ”tin roof” pinning potential in the transverse direction.

Depinning and dynamical ordering of periodic lattices, such as vortex lattices (VL) driven by an external force in presence of static point disorder has been widely studied these last years. A large velocity expansion Koshelev  showed that the motion reduces the effects of disorder. In this approximation the effects of disorder on moving VL results only in an additional ”shaking” temperature $`T_{sh}1/v`$ with expectation of moving crystalline VL at large velocities. It was however shown in Giamarchi  that, due to the components of the disorder and displacements transverse to the motion, a static disorder persists even at large velocity, and leads to a moving glass phase (MG). The motion of particles in the MG occurs along rough static channels. Depending on the velocity the channels can be coupled or decoupled due to dislocations, thus giving two generic moving glasses. For weak disorder or large velocity in $`d=3`$, channels are elastically coupled leading to a topologically ordered structure: the moving Bragg glass (MBG) Giamarchi  . For larger disorder or in $`d=2`$ channels are decoupled by dislocations giving a moving transverse glass (MTG) or smectic Giamarchi ; Balents  . The predicted channel structure and dynamical ordering was confirmed by numerical Simul ; Jensen ; Gronbech ; Fangohr  and further analytical Scheidl  studies and observed experimentally Marchevsky .

Within the MG model Giamarchi , the situation for correlated disorder such as heavy-ion columnar tracks has been recently investigated Chauve . For all velocities at $`T=0`$ or for weak velocities at $`T>0`$ it was predicted a moving Bose glass phase (MBoG) characterized by a transverse critical force and a diverging tilt modulus due to localization effect arising from the columns. This novel feature specific to correlated disorder results in a vanishing tilt response below a critical transverse field. This so-called dynamical transverse Meissner effect (DTME) is a crucial consequence of the MG theory; its observation should provide an unambiguous signature of the MBoG phase. However, up to now there is no experimental nor numerical simulation evidence for this phase.

This Letter reports 3D numerical simulation of field tilting in a moving layered vortex system at zero temperature in a random pinning landscape of columnar defects. The elastically moving vortex lines confirms the existence of MBoG displaying DTME. Above a critical angle, the vortex lines exhibit a ”kink” structure along their length. This results from the existence of an effective static pinning landscape which is almost periodic in the transverse direction to the vortex flow. Finite size effects are discussed within a simple model of a single flux line in a static periodic potential. In particular, it confirms the existence of a moving ”kink” structure observed in our simulation.

We model a stack of $`N_z`$ Josephson-coupled parallel superconducting planes of thickness $`d`$ with interlayer spacing $`s`$. Each layer in the $`(x,y)`$ plane contains $`N_v`$ pancake vortices interacting with a random pinning background of 1D correlated disorder, namely $`N_p`$ columnar pins parallel to the $`z`$ direction. At $`T=0`$ the overdamped equation of motion of a pancake $`i`$ in position $`𝕣_i(z)`$ reads

$$\eta \frac{d𝐫_𝐢}{dt}=\underset{ji}{}_iU^{vv}(\rho _{ij},z_{ij})\underset{p}{}_iU^{vp}(\rho _{ip})+𝐅^L+𝐅^{tilt}(z)$$

where $`\rho _{ij}`$ and $`z_{ij}`$ are the components of $`𝐫_{\mathrm{𝐢𝐣}}=𝐫_𝐢𝐫_𝐣`$ in cylindrical coordinates, $`\rho _{ip}`$ is the in-plane distance between the pancake $`i`$ and a pinning site in the same layer at $`𝐫_𝐩`$, and $`_i`$ is the 2D gradient operator acting in the $`(x,y)`$ plane. The viscosity coefficient is $`\eta `$, $`𝐅^L=F^L\widehat{𝐱}`$ is the Lorentz driving force due to an applied current and $`𝐅^{tilt}(z)`$ is the surface force due to the field tilting away from the $`z`$ axis in the $`y`$ direction. This force acts as a torque on each flux line, i.e. $`𝐅^{tilt}(z=0)=𝐅^{tilt}(z=N_zs)=F^{tilt}\widehat{𝐲}`$ and $`𝐅^{tilt}(z)=\mathrm{𝟎}`$ for pancakes in the bulk. The tilting force modulus is defined by $`F^{tilt}=ϵ^2\varphi _0H_y/4\pi =8\pi ϵ^2ϵ_0\lambda _{ab}^2H_y/\sqrt{3}a_0^2H_z`$, where $`ϵ_0=(\varphi _0/4\pi \lambda _{ab})^2`$, $`\lambda _{ab}`$ is the in-plane magnetic penetration depth, $`a_0`$ is the average vortex distance, $`H_y`$ is the transverse field component, and $`ϵ`$ is the anisotropy parameter. The intra-plane vortex-vortex repulsive interaction is given by a modified Bessel function $`U^{vv}(\rho _{ij})=2ϵ_0dK_0(\rho _{ij}/\lambda _{ab})`$. Following Ryu , the inter-plane attractive interaction between pancakes in adjacent layers of altitude $`z`$ and $`(z+s)`$ reads $`U^{vv}(\rho _{ij},z_{ij}=s)=(2sϵ_0/\pi )[1+ln(\lambda _{ab}/s)][(\rho _{ij}/2r_g)^21]`$ for $`\rho _{ij}2r_g`$ and $`U^{vv}(\rho _{ij},z_{ij}=s)=(2sϵ_0/\pi )[1+ln(\lambda _{ab}/s)][\rho _{ij}/r_g2]`$ otherwise; in this expression $`r_g=\xi _{ab}/ϵ`$, where $`\xi _{ab}`$ is the in-plane coherence length. Finally, the attractive pinning potential is given by $`U^{vp}(\rho _{ip})=\alpha A_pe^{(\rho _{ip}/R_p)^2}`$, where $`A_p=(ϵ_0d/2)ln[1+(R_p^2/2\xi _{ab}^2)]`$ as proposed in Mkrtchyan  and $`\alpha `$ is a tunable parameter. We consider periodic boundary conditions of $`(L_x,L_y)`$ sizes in the $`(x,y)`$ plane while open boundaries are taken in the $`z`$ direction. All details about our method for computing the Bessel potential with periodic conditions can be found in Olive . Molecular Dynamics simulation is used for $`N_v=30`$ vortex lines in a rectangular basic cell $`(L_x,L_y)=(5,6\sqrt{3}/2)\lambda _{ab}`$, and for a number of layers varying from $`N_z=19`$ up to $`N_z=149`$. The number of columnar pins is set to $`N_p=30`$. We consider the London limit $`\kappa =\lambda _{ab}/\xi _{ab}=90`$, with an average vortex distance $`a_0=\lambda _{ab}`$, and $`d=\mathrm{2.83\; 10}^3\lambda _{ab}`$, $`s=\mathrm{8.33\; 10}^3\lambda _{ab}`$, $`R_p=0.22\lambda _{ab}`$, $`ϵ=0.01`$, $`\eta =1`$. These values are coherent with strong anisotropy layered superconductors with $`B_z1000G`$. The choice of the second order Runge-Kutta algorithm time iteration step $`\delta t`$ is dictated by the dominant force $`F^L`$, and we take $`\delta t=\mathrm{7\; 10}^2t_0`$, where $`t_0=\eta \xi _{ab}/F^L`$. We choose the tunable pinning parameter $`\alpha =1/25`$ so that the maximum pinning force is $`F_{max}^{vp}F_0/5`$ where $`F_0=2ϵ_0d/\lambda _{ab}`$ is the unit force defined by the Bessel interaction. Finally, the driving force applied along a principal vortex lattice direction $`x`$ is set to $`F^L2.8F_0`$. This corresponds to the fully elastic flow limit since $`F^L25F_c^L`$ where $`F_c^L`$ is the critical Lorentz force along $`x`$.

The ”experimental” procedure is the folllowing: we start in 2D with the $`(x,y)`$ plane by randomly throwing $`N_v`$ pancakes and $`N_p`$ gaussian pins, and relaxation with zero Lorentz force yields a vortex structure with dislocations. The Lorentz force is then slowly increased far in the elastic phase up to $`F^L2.8F_0`$ with steps $`\delta F^L\mathrm{3.5\; 10}^3F_0`$ every $`\mathrm{3.6\; 10}^4t_0`$. The successive dynamical regimes observed in our 2D vortex system are the following: (i) pinned regime where all vortices have zero velocity; (ii) plastic channels flowing through pinned regions and where the motion is periodic (washboard frequency) as seen in Jensen ; Gronbech ; (iii) plastic turbulent flow with no stationnary vortices as seen in Gronbech ; Fangohr ; (iv) weakly decoupled channels named MTG since quasi long range longitudinal order exists in each channel; (v) elastic phase where all vortices have the same average velocity and named MBG since no dislocation appears and the motion occurs through rough static channels even in the high driving phase $`F^L25F_c^L`$ studied below. Our 2D numerical results are therefore in agreement with the moving glass theory developed in Giamarchi . Finally, note that they are consistent with the phase diagram region corresponding to the intermediate pinning of Ref Fangohr, .

Now starting from this situation, we extend the vortices and the pins in the $`z`$ direction in order to obtain 3D straight vortex lines moving elastically along the $`x`$ direction through a random columnar pin landscape. Then, starting from this situation, we now slowly tilt the magnetic field $`𝐇`$ away from $`z`$ in the $`y`$ direction by increasing the $`H_y`$ component transversally to the vortex flow. For several system sizes, Fig. 1 displays the vortex line response to low field tilting, i.e. the average vortex line inclination $`\mathrm{tan}\theta _B=B_y/B_z`$ versus the field inclination $`\mathrm{tan}\theta _H=H_y/H_z`$, compared with the linear response of the pinning free vortex system. It shows a partial screening of the transverse component $`H_y`$ of the field, which tends to be total for infinite vortex line length since the slopes tend towards zero when $`N_z\mathrm{}`$. This transverse screening is followed by a transition at a critical field $`H_y^c/H_z0.50`$, above which the average vortex line inclination recovers the one obtained without any pinning.

Therefore, Fig. 1 shows the existence of DTME which is the very signature of the MBoG phase as predicted in Chauve . As a good illustration of DTME, we show in the insert of Fig. 1, the vortex line configuration of the MBoG compared with the non-disordered lattice obtained for the same field tilting. One clearly sees that below the transition the transverse field enters at the surface without penetrating into the bulk, in opposition to the pinning free vortex system. So below the critical transverse field, the vortex lines remain parallel to the columnar defects even though they are moving at enough large velocity.

We shall now examine the situation above the critical transverse field. Far above $`H_y^c/H_z`$ and whatever the vortex line size is, the average vortex line inclination recovers the one obtained without any pinning (not shown). A more interesting case happens in the intermediate transverse field range above the transition and shown for $`N_z=49`$ layers in Fig. 2. Increasing the transverse field above its critical value $`H_y^c/H_z0.50`$ gives rise to plateaux in the $`\mathrm{tan}\theta _B`$ vs $`H_y/H_z`$ curve.

A consequence of these plateaux is that decreasing the transverse field generates hysteresis in the vortex inclination. A typical vortex line configuration above the transition is shown for $`N_z=149`$ layers in Fig. 3(a). A modulated structure $`z(y)`$ appears along their length.

This modulation is more obvious looking at the derivative $`dz/dy`$ shown for all the vortex lines in Fig. 3(b). By using the minima and maxima of the derivative we may roughly define two parts of small and big inclinations in the vortex lines with respect to $`z`$. The small ones are termed ”pinned” regions and are linked by ”unpinned” regions where the vortex inclination is bigger. The materialisation of the ”pinned” regions is shown in grey in Fig. 3(a). We visualize in this way a static pinning landscape in the transverse direction $`y`$ with effective columnar pin width corresponding to twice the input pinning ”radius” $`R_p`$. Each plateau of Fig. 2 can be thereby labelled with an odd number corresponding to the number of ”pinned” regions passed through by each vortex line. In this way, while increasing the transverse field above the transition, the vortex lines firstly pass through 3 pinned regions for a given range of $`H_y`$ defining the ”P3” plateau, then the system jumps to the following plateau ”P5” where the vortex lines pass through 5 pinned regions, and so on. The same process holds while decreasing the transverse field generating thereby hysteresis. Finally, the central point displayed in Fig. 3(a) is that the vortex lines which are moving in the $`x`$ direction perpendicularly to the figure plane see an effective static tin roof potential in the transverse direction $`y`$. This is a direct consequence of the existence of static pinned channels of motion Giamarchi ; Chauve , and clearly shows that the action of the disorder in the moving system cannot be simply reduced Koshelev  to a temperature. Note that this effective potential is not strictly periodic because of the channel roughness.

We now turn to the discussion of the finite size effects in our simulation by considering the tilt response of a single static elastic vortex line in a tin roof pinning potential at $`T=0`$. Such a simple model allows to get a valuable insight into the transverse properties of the moving vortex lattice in presence of correlated disorder, e.g. the TME and the kink vortex line structure. The energy of a line of length $`L`$ is given by:

$`E(u)={\displaystyle _0^L}𝑑z\left({\displaystyle \frac{c}{2}}\left({\displaystyle \frac{du}{dz}}\right)^2+V_{eff}(u)\right)+f\left(u(L)u(0)\right),`$ (1)
where $`u(z)`$ is the one-dimensional displacement field in the $`y`$ direction, $`c=ϵ^2ϵ_0`$ is the elastic constant, $`fH_y`$ is a surface force, and $`V_{eff}(u)=g\left[1\mathrm{cos}(4\pi u/\sqrt{3}a_0)\right]`$ with $`g`$ being a constant. There exists a ground state $`u_{L,f}`$ which minimizes $`E(u)`$. We calculate the tilt response of the vortex line defined by $`\mathrm{tan}\theta =(u_{L,f}(L)u_{L,f}(0))/L`$. In the limit $`L\mathrm{}`$, we find a critical force $`f_c=4\sqrt{gc}/\pi `$ below which $`u_{L=\mathrm{},f}0`$, giving $`\mathrm{tan}\theta =0`$ and hence TME. For $`ff_c`$ the vortex line is tilted and its ground state $`u_{L=\mathrm{},f}(z)`$ are expressed in terms of Jacobian elliptic functions which generate kink vortex line structures as observed in our simulation (see Fig. 3(a)). The tilt response is shown in the insert of Fig. 4 together with the pure elastic response, i.e. $`g=0`$. In the latter case, the ground states are linear functions and hence the vortex is a straight line of slope $`\mathrm{tan}\theta =f/c`$. The kink vortex line structure obviously rely only on the periodic potential. Note that $`f^{}=\sqrt{3}f/8\pi c`$ is the equivalent quantity of $`H_y/H_z`$. For a finite size line and $`f=0`$, the minimum energy solution is $`u_{L,f=0}0`$, i.e. the line is pinned. If $`f`$ is nonzero but sufficiently weak, both extremities of the line are pulled away from its minimum energy configuration in a continuous way, yielding a nonzero tilt response as observed in our simulation (see Fig. 1 and insert). Furthermore, for given values of $`f`$, the vortex line abruptly stretches to a more favourable position resulting in jumps in the tilt response. We show in Fig. 4 the tilt response computed for $`L=0.4a_0`$ which is the equivalent size of the 49 layers in our simulation results shown in Fig. 2; the straight line displays the pure elatic response $`g=0`$. We therefore find plateaux and hysteresis which are very close to the ones observed in our simulation and shown in Fig. 2.

In conclusion, our numerical simulation results give the first confirmation of the existence of the MBoG phase displaying a DTME developed in Chauve . Furthermore, above the critical transverse field, the simulation revealed a ”kinked” vortex line structure resulting from an almost periodic effective static pinning potential in the transverse direction.

###### Acknowledgements.

We are grateful to the Laboratoire de Mathématiques et Physique Théorique - CNRS UMR 6083 Université F. Rabelais - Tours (France) for our extensive use of their computers.


text:
# Presupernova Evolution of Rotating Massive Stars and the Rotation Rate of Pulsars

## 1. Introduction

Stars more massive than about $`810\mathrm{M}_{}`$ end their life as supernova, leaving a compact remnant - a neutron star or a black hole. Generally speaking, the evolution of such a star follows a well understood path of contraction to increasing central density and temperature. This path of contraction is interrupted by nuclear fusion – first hydrogen to helium, then helium to carbon and oxygen, followed by carbon, neon, oxygen, and silicon burning, until finally a core of iron is produced. Each fuel burns first in the center of the star, then in one or more shells. In Table 1 we summarize the burning stages and their durations for a $`20\mathrm{M}_{}`$ star. The time scale for helium burning is about ten times shorter than that of hydrogen burning, mostly because of the lower energy release per unit mass. The time scale of the burning stages and contraction beyond central helium-burning is greatly reduced by thermal neutrino losses that carry away energy in situ, instead of requiring that it be transported to the stellar surface by diffusion or convection. These losses increase with temperature as roughly $`T^9`$ (see Woosley, Heger & Weaver for a more extended review). When the star has built up a large enough iron core, exceeding its effective Chandrasekhar mass, it collapses to form a neutron star or a black hole. A supernova explosion may result (e.g., Colgate & White 1966), or, in rare cases, a gamma-ray burst (Woosley 1993; MacFadyen & Woosley 1999; MacFadyen, Woosley, & Heger 2001).

Unevolved massive stars, however, are known to rotate rapidly at several times 10 % of their break-up velocity at the surface. Characteristic values of the equatorial rotation velocity on the main sequence are around $`200\mathrm{km}\mathrm{s}^1`$ (Fukuda 1982; and other contributions in these proceedings). In this paper, we address what happens to the rotation in the interior of these stars as they evolve. In general, the ongoing contraction leads to faster rotation, but transport processes (circulation, viscosity, magnetic fields, etc.) lead to loss of angular momentum from the stellar core. We are particularly interested in the model predictions for the rotational periods of the collapsed remnants left by the supernova.

## 2. Rotating stars

### 2.1. Angular momentum transport

On the main sequence, massive stars quickly settle into rigid rotation (e.g., Zahn 1992; Talon et al. 1997, Meynet & Maeder 2000; Maeder & Meynet 2000). Contraction increases the rotational energy faster than gravitational binding energy and would eventually make the core unstable to triaxial deformations (Ostriker & Bodenheimer 1973). Without the transport of angular momentum, for typical initial rotation rates, this would happen during central helium burning, but clearly the transport of angular momentum in the stellar interior cannot be neglected. Important processes for angular momentum transport include hydrodynamic instabilities (convection, shear, circulation, etc.; see, e.g., Endal & Sofia 1978; Knobloch & Spruit 1983; Heger, Langer, Woosley 2000; Tassoul 2000) and magnetic fields (e.g., Spruit & Phinney 1998; Spruit 2002). The instabilities, however, are inhibited by thermal, and, most importantly, compositional stratification inside the star (e.g., Kippenhahn 1974, Knobloch & Spruit 1983; Maeder & Meynet 1997; Maeder 2003).

Early calculations by Kippenhahn & Thomas (1970); Kippenhahn, Meyer-Hofmeister, & Thomas (1970) and by Endal & Sofia (1976) that ignored angular momentum transport, or only included it in chemically homogeneous regions, did indeed find that the star became secularly unstable at the end of central helium burning or reached critical rotation before central carbon burning. Later models by Endal & Sofia (1978), that considered Eddington-Sweet circulation (Eddington 1924, 1929; Vogt 1925), shear instabilities (e.g., Zahn 1974, 1975; Maeder 1997), GSF instability (Goldreich & Schubert 1967; Fricke 1968), and Solberg-Høiland instability (Wasiutyński 1946) found significant transport of angular momentum that allowed their stars to go well beyond central carbon burning without reaching the limit of secular instability. Still the models rotated about three times faster, after central helium burning, than the more recent models of Heger, Langer & Woosley (2000). This more recent study reduced the efficiency of composition gradients in inhibiting instabilities and of compositional mixing relative to angular momentum transport (Pinsonneault et al. 1989; Chaboyer & Zahn 1992) in order to reproduce observed stellar surface abundance anomalies (e.g., Giess & Lambert 1992; Venn et al. 2002). Calculations by the Geneva Group also found, after central helium burning, rotation rates similar to those of Heger et al. (2000; Maeder, private communication) and the most recent calculations by Hirschi, Meynet, & Maeder (2003) find comparable rotation rates during oxygen shell burning, $`0.16\mathrm{rad}\mathrm{s}^1`$ in $`15\mathrm{M}_{}`$ star of initially $`300\mathrm{km}\mathrm{s}^1`$ compared to $`0.12\mathrm{rad}\mathrm{s}^1`$ for a star with initially $`200\mathrm{km}\mathrm{s}^1`$ in Heger et al. Others have also studied massive star rotation, but concentrated on the main sequence and helium burning (e.g., Urpin, Shalybkov, & Spruit 1996; Talon & Zahn 1997; Denissenkov, Ivanova, & Weiss 1999).

### 2.2. Mixing

Rotation also causes mixing even in regions that, without rotation, are stably stratified. For a review of mixing processes in stars, see Pinsonneault (1997). Observational evidence for this is the surface enrichment of elements produced deep inside the star. For example, an enrichment of nitrogen with respect to oxygen and carbon, or of helium, and the depletion of boron during central hydrogen burning (e.g., Giess & Lambert 1992; Venn et al. 2002; Fliegner, Langer, & Venn 1996). The enrichment of helium above the convective core during central hydrogen burning can also lead to an increase of the luminosity similar to core overshooting that was often assumed in non-rotating stellar models.

During central helium burning, rotationally-induced mixing can prevent the formation of layering in the convective core (Heger et al. 2002). If the mixing above the convective core is fast enough, carbon and oxygen can even be mixed outwards into hydrogen-rich regions and make primary nitrogen, though more typically in stars of solar metallicity, this is suppressed by a rise in entropy and the composition barrier at the base of the hydrogen shell. Additionally, the ratio of circulation time-scale to evolutionary time is smaller than during hydrogen burning, unless the core is spinning faster ($`\tau _{\mathrm{circ}}\tau _{\mathrm{KH}}(\mathrm{\Omega }/\mathrm{\Omega }_{\mathrm{crit}})^2`$). Note in particular, that hydrogen cannot be mixed down into the core, since it will rapidly burn at the higher temperatures further down in the star. However, for stars of very low or zero metallicity, this barrier is reduced and some primary nitrogen is made (e.g., Heger, Woosley, & Waters 2000a; Meynet & Maeder 2002; Marigo, Chiosi, & Kudritzki 2003). Unfortunately, the “mixing events” in such stars are currently not well modeled and the amount of primary N14 produced is not predicted very reliably.

Since the energy loss to neutrinos accelerates the later burning stages, most rotationally-induced instabilities become too slow to mix radiative regions (Heger et al. 2000), though shear instabilities may still function at convective boundaries (see also Hirschi et al. 2003 in these proceedings).

### 2.3. Evolution till core collapse and pulsar rotation

In Figure 1 we give the evolution of the internal rotation profile for a $`15\mathrm{M}_{}`$ star and, in Table 2, numerical values for the specific angular momentum, radius, and angular velocity at a mass coordinate of $`m=1.7\mathrm{M}_{}`$. This corresponds to the edge of the neutron star that forms later.

During hydrogen burning, angular momentum loss by stellar winds and the contraction of a core bounded by an expanding envelope leads to angular momentum loss from the core. As a gradient in mean molecular weight develops, however, angular momentum gets “trapped” in the core and differential rotation develops between the core and envelope. By the end of central helium burning, rotation in the core is reduced by another factor of two. This is, in part, because the growing core incorporates material of low angular momentum that had been slowed down during the first dredge-up, and, in part, due to angular momentum transport. After central helium burning, hydrodynamic rotationally-induced instabilities are not fast enough for significant angular momentum transport (see above) and angular momentum is essentially locally trapped in the cores of different composition. Only convection is fast enough to redistribute angular momentum over larger regions of the stars and the rigidly rotating convective shells leave an imprint in the rotation profile that are each characterized by a steep drop of the specific angular momentum at its base since, typically, the inner radius of such a shells is much smaller than its outer radius.

In the case of a $`15\mathrm{M}_{}`$ star (solar metallicity) with an initial rotation rate of $`200\mathrm{km}\mathrm{s}^1`$ the core spins fast enough to form a neutron star in excess of critical rotation (last line of Table 2) if no further loss of angular momentum occurs. Since rotation close to breakup clearly is in contradiction to the rotation rates of young pulsars (Table 3) some additional spin-down is necessary. However, this much rotation in collapsing massive star cores would be preferable for the “collapsar” model for gamma-ray bursts (Woosley 1993). For some time, the $`r`$-mode instability (e.g., Lindblom, Tohline, & Vallisneri 2001) seemed to be a promising explanation for spin-down, so that both slow pulsars and collapsars could result from rapidly rotating core, but recent evaluation (Arras et al. 2003) seem to indicate that it may be not efficient enough (see also Woosley & Heger 2003 in these proceedings). An alternative possibility for removing this discrepancy may be angular momentum transport by magnetic fields.

## 3. Magnetic stars

Spruit & Phinney (1998) assumed that differential rotation would torque an initially arbitrarily weak radial field component into a toroidal field, that, when strong enough, would become unstable, producing a new radial field that is wound up again. With such strong fields they obtained a star that was kept in rigid rotation until about central carbon burning and estimated pulsar rotation rates at birth $`100\mathrm{s}`$ result. If mass loss is taken into account, the resulting spin would have been even lower by orders of magnitude.

In a more recent approach, Spruit (2002) introduced a dynamo model that considers stabilization of the stratification by thermal and chemical stratification (see also Spruit 2003 in these proceedings). It is assumed that the dynamo reaches a steady state field on a time-scale short compared to that on which the structure and rotation profile of the star evolve. From this dynamo and its equilibrium field an effective diffusivity for chemical mixing and an effective viscosity for angular momentum transport can be derived (Spruit 2002). We have implemented both as time-dependent diffusive processes (Heger, Woosley, & Spruit 2003). The results for a $`15\mathrm{M}_{}`$ model with the same initial properties as in Figure 1 and Table 2, but now including the dynamo process, is given in Figure 2 and Table 4. At the end of central hydrogen burning the star is essentially still in solid body rotation – the small steps seen in Figure 1 between $`46\mathrm{M}_{}`$, marking the composition gradient left behind by central hydrogen burning, are not present – and it rotates significantly slower. By the end of central helium burning, the core is slowed down even further. At this point it rotates at less than a tenth the speed of the non-magnetic star. From here to core collapse the angular momentum in the inner $`2\mathrm{M}_{}`$ is reduced by an additional factor $`3`$. Despite significant differential rotation, essentially no angular momentum is transported across the edge of the helium core where entropy strongly increases. In Table 4 we also give the central angular velocity at late evolution stages until core collapse along with the equilibrium radial and toroidal magnetic fields. This compilation demonstrates an important property of of Spruit’s model: the radial field component is much weaker than the toroidal field component. This is also one of the significant differences with the model of Spruit & Phinney (1998). From these stellar models we obtain neutron star rotation rates of $`47\mathrm{ms}`$ for $`1525\mathrm{M}_{}`$ stars assuming all the angular momentum in the core were conserved to for a neutron star (Table 5) and the remnants of more massive stars are rotating faster. These rotation rates are still faster than that of observed young pulsars (Table 4) and additional spin-down mechanisms are still needed to reduce the discrepancy (see Woosley & Heger 2003 in these proceedings). In the derivation of the dynamo model by Spruit (2002) several quantities of “comparable size” have been set equal, introducing parameters “of the order unity” in the theory. Most significant for the angular momentum transport and the presence of effective barriers, like the edge of the helium core mentioned above, are the influence of stabilization of the stratification by thermal and compositional buoyancy. We have varied both, as well as the overall magnetic stress from the dynamo by a factor 10 up and down (Table 6). Even these large variation only result in variations of the neutron star rotation by about a factor $`2`$. The reason for the smallness of the effect is that the stress resulting from the dynamo scales as the 6th power of the shear rate.

## 4. Summary and conclusions

Without the action of magnetic fields, hydrodynamic rotationally-induced instabilities and convection alone (Figure 1) do not transport angular momentum efficiently enough to avoid forming pulsars rotating at break up (assuming that most of the angular momentum is conserved during the collapse; Table 2). When a dynamo process like that of Spruit (2002) is included in the models (Table 4 and Figure 2), one obtains pulsar periods more than a factor 10 slower, about $`47\mathrm{ms}`$ (Table 5). Though these periods may still increase by $`20\%`$ due to neutrino losses from the proto-neutron star, they are still too fast compared even to the fastest-rotating observed young pulsars (Table 3) and additional spin-down mechanisms may still be needed in single stars (see Woosley & Heger 2003 in these proceedings). However, the rotation rates found in the magnetic stellar models would be too slow for most current gamma-ray burst models which require rapidly rotating stellar cores. However, interacting close binary stars may lead to much faster or slower core rotation, and, e.g., to pulsars of up to a second spin period (Langer et al. 2003).

This research has been supported by the NSF (AST 02-06111), NASA (NAG5-12036), and the DOE Program for Scientific Discovery through Advanced Computing (SciDAC; DE-FC02-01ER41176). AH is supported, in part, by the Department of Energy under grant B341495 to the Center for Astrophysical Thermonuclear Flashes at the University of Chicago, and a Fermi Fellowship at the University of Chicago.

## References

Arras, P., Flanagan, E. E., Morsink, S. M., Schenk, A. K., Teukolsky, S. A., & Wasserman, I. 2003, submitted to ApJ, astro ph-0202345

Chaboyer, B., Zahn, J.-P. 1992, A&A, 253, 173

Colgate, S. A., & White, R. H. 1966, ApJ, 143, 626

Denissenkov, P. A. and Ivanova, N. S. and Weiss, A. 1999, A&A, 341, 181

Eddington, A. S. 1925, Observatory, 28, 73

Eddington, A. S. 1929, MNRAS, 90, 54

Endal, A. S. & Sofia, S. 1976, ApJ, 210, 184

Endal, A. S. & Sofia, S. 1978, ApJ, 220, 279

Fliegner, J., Langer, N., Venn, K. A. 1996, A&A, 308, L13

Fricke, K. 1968, ZfAp, 68, 317

Fukuda, I. 1992, PASP, 94, 271

Gies, D. R., Lambert, D. L. 1992, ApJ, 387, 318

Goldreich, P., Schubert, G. 1967, ApJ, 150, 571

Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368

Heger, A., Woosley, S. E., Waters, R. 2000a, in *The First Stars. Proceedings of the MPA/ESO Workshop*, eds. A. Weiss, T. G. Abel, V. Hill, Springer, p. 121

Heger, A., Woosley, S. E., & Spruit, H. 2003, ApJ, in preparation

Hirschi, R., Meynet, G., Maeder. A. 2003, in Proc. IAU Symp. 215, *Stellar Rotation*, eds. A. Maeder, P. Eenens, (this volume)

Kippenhahn, R. & Thomas, H. C. 1970, in Proc. IAU Coll 4, *Stellar Rotation*, ed. A Slettebak, Reidel, Dortecht, p. 20

Kippenhahn, R., Meyer-Hofmeister, E., Thomas, H. C. 1970, A&A, 155

Kippenhahn, R. 1974, in *Late Stages of Stellar Evolution*, eds. R. J. Taylor, J. E. Hesser, IAU Symp. 66

Knobloch, E. & Spruit, H. C. 1983, A&A, 125, 59

Langer, N., Yoon, S.-C., Petrovic, J., Heger, A. 2003, in Proc. IAU Symp. 215, *Stellar Rotation*, eds. A. Maeder, P. Eenens, (this volume)

Lattimer, J. M. & Prakash, M. 2001, ApJ, 550, 426

Lindblom, L., Tohline, J. E., & Vallisneri, M. 2001, Phys. Rev. Lett, 86, 1152

MacFadyen, A., Woosley, S. E. 1999 ApJ, 524, 262

MacFadyen, A., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 410

Maeder, A. 2003, A&A, 2003, in press

Maeder, A. 1997, A&A, 321, 134

Maeder, A. & Meynet, G. 1997, A&A, 321, 465

Maeder, A. & Meynet, G. 2000, ARA&A, 38, 143

Marigo, P., Chiosi, C., Kudritzki, R.-P. 2003, A&A, in press; astro-ph/0212057

Marshall, F. E., Gotthelf, E. V., Zhang, W., Middleditch, J., & Wang, Q. D. 1998, ApJ, 499, 179

Meynet, G. and Maeder, A. 2000, A&A, 361, 101

Meynet, G. and Maeder, A. 2002, A&A, 390, 561

Ostriker, J. P. & Bodenheimer, P. 1973, ApJ, 180, 171

Pinsonneault, M. H. 1997, ARA&A, 35, 557

Pinsonneault, M. H., Kawaler, S. D., Sofia, S., Demarque, P. 1989, ApJ, 338, 424

Spruit, H. C. 2002, A&A, 381, 923

Spruit, H. C. 2003, in Proc. IAU Symp. 215, *Stellar Rotation*, eds. A. Maeder, P. Eenens, (this volume)

Spruit, H. C. & Phinney, E. S. 1998, Nature, 393, 139

Talon, S., Zahn, J.-P. 1997, A&A, 317, 749

Talon, S., Zahn, J.-P., Maeder, A., Meynet, G. 1997, A&A, 322, 209

Tassoul J.-L. 2000, *Stellar Rotation*, Cambridge University Press, ISBN 0521772184

Urpin, V. A., Shalybkov, D. A., Spruit, H. C. 1996, A&A, 306, 455

Venn, K.A., Brooks, A.M., Lambert, D., Lemke, M., Langer, N., Lennon, D.J., Keenen, F.P. 2002, ApJ, 565, 571

Venn, K. A., Lambert, D. L., Lemke, M. 1996, A&A, 307, 894

Vogt, H. 1925, Astr. Nachr., 223, 229

Wasiutyński, J. 1946, Ap. Norveg., 4, 1

Woosley, S. E. 1993, ApJ, 405, 273.

Woosley, S. E., Heger, A. 2003, in Proc. IAU Symp. 215, *Stellar Rotation*, eds. A. Maeder, P. Eenens, (this volume)

Zahn, J.-P. 1974, in *Stellar instability and Evolution*, eds. P. Ledoux, A. Noels, A. W. Rogers, IAU Symp. 59, Redel, Dortrecht, p. 185

Zahn, J.-P. 1975, Mém. Soc. Roy. Sci. Liège, 6e Serie, VIII, p. 31

Zahn, J.-P. 1992, A&A, 256, 115