text:
# Monopolar Pulsar Spin-Down

## 1 Introduction

Pulsars are accepted to be rapidly rotating and highly magnetized neutron stars powered by the loss of rotational kinetic energy. Yet, the mechanisms through which pulsars lose their rotational energy are not very well understood. The discovery of pulsars was preluded by the rotating magnetic dipole model of Pacini (pacini67 (1967)), which predicted a spin-down relation $`\dot{\nu }\nu ^3`$, where $`\nu `$ is the rotational frequency. Pulsar spin-down has then been studied through the power-law differential equation,

$$\dot{\nu }=k\nu ^n,$$
(1)
where $`k`$ is a constant and $`n`$ is called the “braking index”, determined by the physical mechanisms spinning-down the star. For a pure magnetic dipole $`n`$ is equal to 3 (Gold gold (1968); Pacini pacini68 (1968)), while a pure gravitational or electromagnetic quadrupole one has $`n=5`$ (Ostriker & Gunn ostriker (1969); Ferrari & Ruffini ferrari (1969)). The expansion of the Larmor equation into electromagnetic moments higher than dipolar gives $`n>3`$. Departures from these basic models have been studied. For example, the deformation of magnetic field lines by corotation with the magnetosphere might produce $`1n3`$ (Manchester & Taylor manchester77 (1977)), pulsar winds $`n<3`$ (Blandford & Romani blandford (1988); Manchester et al. manchester85 (1985)), while magnetic field decay (Chanmugam & Sang chanmugam89 (1989)) and the alignment between rotation and magnetic axis tend to give $`n>3`$ (Goldreich goldreich (1970)). Physical processes such as phase transitions in the interior of the fast rotating neutron stars can produce changes in the moment of inertia $`I`$ and led also to deviations from $`n=3`$ even if the pulsar is spinning-down by pure dipolar radiation (Chubarian et al. chubarian2000 (2000); Glendenning et al. glendenning97 (1997)). However, these phase transitions are significant only for very fast pulsars, with periods $`P<1.5`$ ms. For pulsars with periods $`P\begin{array}{c}>\\ \end{array}3`$ ms, these changes in the moment of inertia can be neglected. But while physical mechanisms to produce $`n<3`$ have been proposed, the reconciliation of pulsar evolution with $`n<3`$ remains to be studied. This question, addressed to the case of young classical pulsars, is one of the main objectives of this work.

The braking index $`n`$ is frequently defined in terms of observational quantities, with a definition generalized to a second braking index, $`m`$:

$$n\frac{\ddot{\nu }\nu }{\dot{\nu }^2},m\frac{\stackrel{\mathrm{}}{\nu }\nu ^2}{\dot{\nu }^3}.$$
(2)
Accurate measurement of both braking indices requires long-term timing measurements, a difficult task due to the existence of glitches and timing noise in young pulsars, affecting the measurement of $`\ddot{\nu }`$ and $`\stackrel{\mathrm{}}{\nu }`$. Out of more than 1300 detected pulsars, only a handfull have reliable braking index measurements (Table 1). Still, measurements continue to be made and objects like PSR J1119–6127 and PSR J1846–0258 might probably be recent additions to the list (Camilo et al. camilo00 (2000); Mereghetti et al. mereghetti (2002)). The second braking index $`m`$ is even more difficult to measure and has only been estimated for PSR B1509–58 (Kaspi et al. kaspi (1994)) and the Crab pulsar (Lyne et al. lyne88 (1988)) (Table 1). It is well-known that measured values of $`n`$ and $`m`$ are below, but close, to those of a magnetic dipole in vacuum ($`n=3`$ and $`m=15`$).

Pulsar evolution is mostly studied through the $`P\dot{P}`$ diagram, which shows a sparse distribution of pulsars clustered roughly at $`\mathrm{log}P0.3`$ and $`\mathrm{log}\dot{P}15`$. In principle it should be possible to link the bulk of the pulsar population in the $`P\dot{P}`$ diagram with the younger pulsars through evolutionary tracks consistent with a general spin-down equation. However, neither the magnetic dipole model or a constant braking index equation can link these groups and it has been argued that pulsars might be born with small period derivatives, in a region of the $`P\dot{P}`$ diagram where no pulsar detection has been reported (Camilo camilo96 (1996)). Alternatively, magnetic dipole moments might decay during the pulsar life (Camilo camilo96 (1996); Tauris & Konar Tauris (2001)). On these lines, Colpi et al. (colpi2000 (2000)) studied the distribution and evolution in the $`P\dot{P}`$ diagram of anomalous X-ray pulsars (AXPs), considering dipolar spin-down of magnetars. They conclude that magnetic field decay is required to account for the observed distribution and X-ray luminosities.

We propose here a spin-down model derived from the very general spin-down law $`\dot{\nu }=f(\nu ,t)`$. Although the derivation is empirical, the model accomodates the usual spin-down mechanisms, namely magnetic dipole radiation and gravitational radiation. In addition, a monopolar term is introduced, aiming to reconcile braking indices with trajectories in the $`P\dot{P}`$ diagram. We calculate evolutionary trajectories for PSR B1509–58, PSR B0540–69, the Vela pulsar (PSR B0833–45) and the Crab pulsar (PSR B0531+21). The model is presented in Section 2, while in Section 3 we apply the model to the four pulsars just mentioned, assuming time-independent spin-down mechanisms. In Section 4 we study the model considering time-dependent spin-down mechanisms. In Section 5, a discussion of possible physical mechanisms related to the spin-down terms is given, to lead to the concluding Section 6.

## 2 The Model

The basic assumption of this work is that the frequency derivative $`\dot{\nu }`$ can be described as a function of frequency $`\nu `$ and time $`t`$ only:

$$\dot{\nu }=f(\nu ,t).$$
(3)
We then assume the following properties:

1. $`f>0`$ for all $`(\nu ,t)`$, as $`f`$ describes energy and angular momentum losses.
2. $`f`$ is a continuous function of time. This requires that the cumulative effect of glitches on evolutionary paths in the $`P\dot{P}`$ diagram can be neglected.
3. $`f`$ is an antisymmetric function, $`f(\nu ,t)=f(\nu ,t)`$. This is equivalent to ascribe a sign to the frequency $`\nu `$ and its derivative $`\dot{\nu }`$: $`\nu `$ and $`\dot{\nu }`$ have the same signs when the pulsar is spinning-up and opposite signs when the pulsar spins-down. This conditions will be further discussed below.

The model consists then of a simple Taylor expansion of $`f`$, restricted to the three lowest order terms:

$$\dot{\nu }=s(t)\nu r(t)\nu ^3g(t)\nu ^5.$$
(4)
Higher order terms can be neglected for the frequencies measured in known isolated pulsars111Binary and millisecond pulsars are not considered here.. A variable moment of inertia, $`\dot{I}0`$, would introduce an extra positive monopole term $`\dot{I}/2I`$ that can be incorporated in $`s(t)`$. This would affect only very fast pulsars and we will consider $`I`$ constant in here. Equation (4) represents a very general spin-down model for isolated pulsars. Beyond its simple empirical derivation, it can represent the loss of rotational energy through physical processes described by standard electromagnetic and gravitational radiation multipolar terms, together with the ad-hoc “monopolar” term. We will argue later that this term has to be associated with other type of processes, like particle acceleration or pulsar winds. We underline that if this model cannot be used to interpret the $`P\dot{P}`$ diagram, then there is no analytical equation of the form given by equation (3), consistent with the required conditions, able to reconcile braking indices with pulsar evolution.

## 3 The stationary multipole spin-down model

Leaving momentarily aside the discussion of the origin of the different terms, we consider now the model given by equation (4) when the coefficients $`g,r`$ and $`s`$ are constant, following Alvarez & Carramiñana (alvarez98 (1998)). These coefficients can be calculated if the frequency and both braking indices are known. Taking the first and second derivative of equation (4) and using equations (2), one arrives to a matrix equation which can be inverted to give:

$$\left(\begin{array}{c}g\nu ^5/(\dot{\nu })\\ r\nu ^3/(\dot{\nu })\\ s\nu /(\dot{\nu })\end{array}\right)=\frac{1}{8}\left(\begin{array}{ccc}\hfill 3& \hfill 3& \hfill 1\\ \hfill 10& \hfill 10& \hfill 2\\ \hfill 15& \hfill 7& \hfill 1\end{array}\right)\left(\begin{array}{c}1\\ n\\ mn^2\end{array}\right).$$
(5)
The factors $`g\nu ^5/(\dot{\nu })`$, $`r\nu ^3/(\dot{\nu })`$ and $`s\nu /(\dot{\nu })`$ represent the fractions of the total energy loss contained in the quadrupole, dipole and monopole terms respectively. Constraining each term to be between 0 and 1 leads to the bounds on their values as function of $`n`$ represented in Figure 1. The value of $`m`$ is restricted to the range

$$\mathrm{Max}(n^2+3n3;n^2+7n15)mn^2+5n5,$$
(6)
while $`1n5`$. The minimum and maximum allowable values for $`m`$ are 1 and 45, which happen for $`n=1`$ and $`n`$=5 respectively. Accordingly, pulsars must be in a restricted region of the $`(n,m)`$ plane defined by equation (6), and shown in Figure. 2. Spin-down evolution is initially dominated by the quadrupole term ($`n5`$), afterwards by the dipole term ($`n3`$) and finally by the monopole term ($`n1`$). Under this model the braking index itself would be a rough qualitative age indicator.

From the known timing parameters, $`\nu `$, $`\dot{\nu }`$, $`n`$ and $`m`$, one determines $`\{g,r,s\}`$ uniquely, from the matrix relation (5), with no degree of freedom. This could only be done for PSR B1509–58, while for the Crab, PSR B0540–69 and Vela pulsars we assumed $`m`$ within the constraints of the model. Evolutionary tracks were calculated for these pulsars, integrating equation (4) backwards one dynamical time and forward to a point where the monopole term dominates the spin-down evolution (Fig. 3).

The results for each of the four pulsars can be summarized as follows:

* the second braking index of the Crab pulsar, $`m=10.23`$, is below the minimum value allowed by the model ($`m10.82`$) and in order to compute evolutionary trajectories we used $`m=11`$. This value of $`m`$ is consistent with the initial period calculated in other models for this pulsar (Mereghetti et al. mereghetti (2002); Glendenning Glendenning (1996)). The Crab is the only pulsar studied here with the precise age known, allowing some specific comparisons. Integrating back to its actual birth (949 years ago) one finds an initial period $`\begin{array}{c}<\\ \end{array}18\mathrm{ms}`$ (Fig 4). If we integrate back one dynamical time (1258 years) from now, we get an initial period of 10 ms, lower but within a factor of two of the period at birth.
* for PSR B1509–58 the measured values of $`m`$ and $`n`$ are consistent with the model. However, the integration of equation (4) could not be carried out one entire dynamical time backwards, as $`𝑑t`$ converges at about 1350 years, i.e., before the present dynamical time (1553 years). Integrating back to a dynamical age of 1 year, we obtain a putative “initial period” of about 14 ms and a $`\dot{P}`$ two orders of magnitude higher than for any of the other three pulsars. This pulsar seems to have been formed in a particular region of the $`P\dot{P}`$ diagram. Its quadrupolar parameter $`g`$ is at least two orders of magnitude higher that for the other pulsars. An initial phase of strong gravitational spin-down would naturally lead to the high $`P`$ and $`\dot{P}`$ observed in this pulsar.
* PSR B0540–69: its second braking index is still unknown so we took $`m=8.0`$, in the middle of the range of allowed values. As the known timing parameters of PSR B0540–69 are similar to those of the Crab, albeit a somewhat longer period, their evolutionary paths are similar, with an initial $`P\dot{P}`$ position of PSR B0540–69 relatively close to that of the Crab.
* Vela pulsar: the frequent glitches of Vela prevent any estimate of the second braking index, and render difficult measuring its first braking index. We assumed $`m=3.5`$, consistent with the low first braking index, $`n=1.4\pm 0.2`$ (Lyne et al. lyne96 (1996)). As $`n<2`$, Vela has already passed its minimum $`\dot{P}`$ value and its predicted evolution is close to a $`n=1`$ line. We note that the track computed for the past evolution of Vela is close of that of PSR B0540–69. However, its motion in the $`P\dot{P}`$ diagram is very slow and when integrating back one dynamical time Vela only reaches $`P50\mathrm{ms}`$ with a braking index almost equal to 2.

Although exact tracks cannot be computed for older pulsars, their known timing parameters not constraining $`\{g,r,s\}`$ sufficiently, we evolved backwards 17 pulsars with $`t_{dyn}10^5\mathrm{years}`$, assuming $`n=1.5`$, in order to locate the region of their likely initial conditions in the $`P\dot{P}`$ diagram (Fig. 5). We stopped the computation when $`t_{dyn}=1\mathrm{year}`$. We note that most of these pulsars arrive to higher $`\dot{P}`$ than all of the youngest pulsars, with the exception of PSR B1509–58. In order to get to $`P\dot{P}`$ positions more consistent with that of the Crab, these pulsars would need to have a braking index closer to 1. A prediction of the stationary multipole model is that $`n\begin{array}{c}>\\ \end{array}1`$ for the majority of the pulsar population.

A more detailed discussion of the stationary model can be found in Alvarez (alvarez02 (1998)). Here we want to point out several inconsistencies of the model, which lead us to reject its non-time dependent version. As it can be seen in Fig. 3, the evolutionary tracks of the youngest pulsars take them to a region of the $`P\dot{P}`$ diagram where no pulsars have been found. Although a large $`\dot{P}`$ might difficult the detection of pulsars located there, the region corresponds to pulsars slower but more energetic than some known radio-pulsars, specifically those with same $`\dot{P}`$ but lower $`P`$. Pulsars in that region should be detectable. We believe these tracks are indicative of the inconsistency between the model and the observations. A second discrepancy is the eventual dominance of the monopole term, which implies $`\nu /\dot{\nu }\mathrm{constant}`$ (Fig 6). As $`n1`$ the dynamical time tends to a constant, $`t_{dyn}1/2s`$, which values between $`2.8\times 10^3`$ years for PSR B0540–69 and $`13\times 10^3`$ years for Vela. The dynamical ages of the rest of the pulsars in the $`P\dot{P}`$ diagram are much larger and cannot be reproduced by the model from the data of these young pulsars, a clear inconsistency between the predicted evolution of the young pulsars and the timing parameters of older pulsars.

Through the study of the stationary multipole model proposed here we reach the general conclusion that there is no time-independent evolutionary equation $`\dot{\nu }=f(\nu )`$ able to reconcile braking indices with the $`P\dot{P}`$ diagram. Note that including a $`\dot{\nu }\nu ^2`$ term222inconsistent with the antisymmetric condition. cannot prevent the eventual dominance of the monopolar term, needed to explain the braking index of Vela (and probably also of PSR B0540–69). The impossibility to have a time-independent expression means that the physical properties involved in the dynamics of pulsars must vary on timescales comparable to their life as active radio-pulsars. The most obvious form of time dependence is through the magnetic moment of the star, namely magnetic field decay or alignment between the magnetic and rotational axes. We consider in the next section the time-dependent multipole model using physically expected time-dependences.

## 4 The time-dependent multipole model

Modelling the distribution observed in the $`P\dot{P}`$ diagram requires considering plausible time dependences in physical parameters involved like the magnetic dipole moment of the neutron star. Observational evidence in favor of magnetic field decay is discussed by Tauris & Konar (Tauris (2001)), who analyzed the relation between characteristic time against rotational energy of the pulsars in the Princeton Catalog (Taylor et al. taylor (1993)). Additional evidence comes from the measurements of their proper motions and from the study of ultramagnetized neutron stars or magnetars. Sang & Chanmugam (sang (1990)) argued that dynamical and kinetic ages of pulsars can only be reconciliated with decaying magnetic fields of the form $`B(t)=B(0)/(1+t/t_c)`$. Further evidence of magnetic field decay is provided by anomalous X-ray pulsars (AXPs) interpreted as magnetars. Colpi et al. (colpi2000 (2000)) required the magnetic field in these neutron star related objects to decay as a power law on timescales of the order of $`10^4`$ years.

We introduce time dependence in our model through the coefficients $`\{r(t),s(t)\}B(0)^2\psi (t/t_c)`$, where $`t_c`$ is the characteristic time-scale for field decay and $`\psi `$ a dimensionless function satisfying $`\psi (0)=1`$. Frequently used functional forms of $`\psi `$ are exponential or inverse linear. For simplicity, we will only consider future evolution for the youngest pulsars, so that the quadrupolar coefficient $`g(t)`$ can be neglected hereafter333Neglecting $`g`$ also means that we do not require knowing $`m`$.. We can then study evolutionary tracks on the $`P\dot{P}`$ diagram writing equation (4) in terms of the period and its derivative as

$$\dot{P}=\left(\frac{r_0}{P}+s_0P\right)\psi (t/t_c).$$
(7)
This equation describes the period evolution of pulsars for a decaying magnetic field. It can be integrated analytically, leading to:

$$P^2(t)=\left(r_0/s_0+P_0^2\right)\mathrm{exp}\left[2s_0t_c\mathrm{\Psi }(t/t_c)\right]r_0/s_0,$$
(8)
where $`P_0P(0)`$ and $`\mathrm{\Psi }(x)=_0^x\psi (u)𝑑u`$. If $`\psi 0`$ for $`tt_c`$, $`\mathrm{\Psi }`$ will generally converge to a value $`\mathrm{\Psi }_{\mathrm{}}`$ and the period will tend to a constant, $`P_{\mathrm{}}`$. This occurs because the spin-down mechanisms disappear as $`\psi `$ approaches zero. If one knows the present values of $`\psi `$ and $`\dot{\psi }`$, the initial parameters $`r_0`$ and $`s_0`$ can be derived from equation (7) and its time derivative:

$`s_0`$ $`=`$ $`{\displaystyle \frac{1}{2\psi }}\left({\displaystyle \frac{\dot{P}}{P}}\right)\left(3n+{\displaystyle \frac{(\dot{\psi })/\psi }{\dot{P}/P}}\right),`$
$`r_0`$ $`=`$ $`{\displaystyle \frac{P\dot{P}}{2\psi }}\left(n1{\displaystyle \frac{(\dot{\psi })/\psi }{\dot{P}/P}}\right),`$ (9)
where $`n`$ is the present braking index and a decay-law means $`\dot{\psi }<0`$. The decay timescale $`t_c`$ is taken here as a free parameter to fit the data, to be determined by selecting physically meaningfull trajectories in the $`P\dot{P}`$ diagram and $`r_0>0`$ and $`s_0>0`$. Note that $`r_0`$ and $`s_0`$ are both positive if

$$n1>\frac{(\dot{\psi })/\psi }{\dot{P}/P}>n3,$$
(10)
which relates $`t_c`$ with the dynamical time and requires $`n1`$ always. On the other hand, old pulsars can now have arbitrarily large braking indices, as $`n`$ can go to $`\mathrm{}`$ when $`\psi /\dot{\psi }\mathrm{cte}`$, $`PP_{\mathrm{}}`$ and $`\dot{P}0`$. The largest loss of rotational energy occurs when $`\dot{\psi }0`$, $`\psi \mathrm{constant}`$ i.e. when $`n`$ is close to one.

In the next sections we present evolutionary tracks for the cases of exponential and linear magnetic field decay.

### 4.1 Exponential decay

Based on the observed deficit of pulsars with periods larger than 1 second, Ostriker & Gunn (ostriker (1969)) proposed a rapid decrease of pulsar radio luminosity as result of exponentially decaying magnetic fields,

$$B(t)=B_0\mathrm{exp}\left(t/t_c\right),$$
(11)
a natural result of ohmic dissipation. In this case the model becomes

$$\dot{P}=e^{2t/t_c}\left(r_0/P+s_0P\right),$$
(12)
and we have $`\mathrm{\Psi }_{\mathrm{}}=1/2`$. For the youngest pulsars we obtain the trajectories shown in Figure 7. Satisfactory trajectories are obtained for $`t_c10^4`$ yrs (Table 2), which correspond to extremely rapid evolution, as the magnetic field becomes negligible in just $`10t_c`$ and pulsars cross the death-line of the $`P\dot{P}`$ diagram in $`10^4`$ to $`10^5`$ years. This timescales are too short to be consistent with estimated birth rates of one pulsar every 90 years (Brazier & Johnston brazier99 (1999)) and supernova rates of one every 10-30 years (Van der Bergh & Tammann vanderbergh (1991)). We conclude that exponentially decaying magnetic fields cannot be used to properly fit the $`P\dot{P}`$ diagram without contradicting present known pulsar and supernova rates.

### 4.2 Inverse linear decay

Sang & Chanmugam (sang (1990)) proposed an inverse linear decay law for the magnetic field of pulsars,

$$B(t)=\frac{B(0)}{1+t/t_c},$$
(13)
to reconcile dynamical and kinetic ages. Using this decay law in the model, relatively short timescales, $`t_c10^4`$ years, are needed to fit the evolutionary trajectories of the four pulsars considered with the $`P\dot{P}`$ distribution, as shown in Figure 8 for the Crab, PSR B1509–58, PSR B0540–69 and Vela. Except for the Vela pulsar, to be discussed below, trajectories can be fitted through the bulk of the pulsar population. Although the timescales needed to fit the data are similar to those for the case of exponential decay, linear decay is far more convenient as pulsars remain active for $`10^310^4t_c`$, crossing the death line of radio emission in some $`10^7`$ years, in reasonable agreement with pulsar and supernova rates.

For the Vela pulsar, the minimum decay timescales consistent with the model, i.e. with positive $`r_0`$ and $`s_0`$, do not manage to turn its evolutionary tracks towards the bulk of pulsars. However, the uncertainty in the measurement of $`n`$ allows us to consider a range of values. Evolutionary trajectories with $`n1.4`$ take Vela to a region of large periods and high period derivatives where no radio pulsars have been found. For $`n=1.2`$, Vela takes about $`2.0\times 10^5`$ yrs to slow down to $`P10\mathrm{s}`$, and 3.6$`\times 10^5`$ years for $`n=1.3`$. Taking these numbers and considering one pulsar like Vela born every 11000 years, we should expect to see between 18 and 27 pulsars in the same region of the $`P\dot{P}`$ diagram occupied by Vela. Indices $`1.5n1.6`$ can produce reasonable tracks (Fig. 9), leading to longer periods than those of the other three young pulsars, but consistent with long-period pulsars present in the $`P\dot{P}`$ diagram. It does not seem feasible to reconcile Vela trajectories with those of the Crab, PSR B1509–58 and PSR B0540–69, which pass though the bulk of the pulsar population without impossing a braking index $`n>1.6`$, inconsistent (at least at the 1$`\sigma `$ level) with observations. So, even though the model can fit the data, the evidence points to the Vela pulsar as a particular type of pulsar. Pulsar evolution might distinguish between Crab-like pulsars, with magnetic-fields decaying at periods $`P_{\mathrm{}}\begin{array}{c}<\\ \end{array}1\mathrm{s}`$, and -less common- Vela-like pulsars, where the magnetic-field decay is slower and becomes spun-down to asymptotic periods between 1 and 10 seconds.

## 5 Physical Spin-down Mechanisms

The main feature of the model presented here is the inclusion of the monopolar term, required to derive reasonable evolutionary tracks considering the measured braking indices. The other two terms of the model have been widely used in the past and are usually interpreted in terms of electromagnetic or gravitational radiation losses. An inclined rotating magnetic dipole in vacuum produces a dipolar term given by

$$\dot{E}_{dip}=\frac{2\mu _{}^2}{3c^3}\mathrm{\Omega }^4,$$
(14)
where $`\mathrm{\Omega }=2\pi \nu `$ is the angular frequency and $`\mu _{}`$ the component of magnetic dipole moment perpendicular to $`\mathrm{\Omega }`$. Pacini (pacini67 (1967)) was the first to equate $`\dot{E}_{dip}=I\mathrm{\Omega }\dot{\mathrm{\Omega }}`$, introducing the basic equation of pulsar dynamics. A quadrupolar term can be associated either to the next electromagnetic moments or to the lowest gravitational radiation multipole moment:

$$\dot{E}_{quad}=\frac{32G}{5c^5}I^2e^2\mathrm{\Omega }^6,$$
(15)
where $`Ie`$ the product of the moment of inertia and equatorial ellipticity of the star. The fact that observed indices are close to 3 indicates that gravitational radiation is not the dominant form of energy loss in pulsars at present, although it might be important in the very early stages of pulsar evolution.

Directly equating the Larmor equation to rotational losses cannot lead to the monopolar term proposed in equation (4). The Larmor equation arises from considering the energy flux of an electromagnetic wave, $`S|E_{rad}|^2\widehat{k}`$, with the radiation field $`E_{rad}`$ proportional to the acceleration of the charged particles, scaling therefore as $`\mathrm{\Omega }^n`$ with $`n`$ even and $`4`$. Even in the framework of the models mentioned in Section 2, it is unlikely that radiative emission processes can give a different dependence on $`\mathrm{\Omega }`$.

On the other hand, there is clear evidence that particle acceleration and massive winds are important energy loss processes in pulsars. This evidence includes the short lifetimes of relativistic electrons in the Crab Nebula and the need of a source to account for its luminosity (as pointed out first by Oort and Walraven oort56 (1956)), observational changes directly observed in the vicinity of the Crab pulsar (Hester et al. hester (1995)) and the X-ray images of the Crab, PSR B1509–58 and the Vela pulsar (Weisskopf et al. weisskopf (2000), Gaensler et al. gaensler (2002), and Helfand et al. helfand (2001) respectively). The monopolar term cannot be of radiative origin, so we propose here its relation with particle and/or mass loss processes. We note we cannot discard a $`\dot{E}\mathrm{\Omega }^3`$ term but this term would be neither necessary nor sufficient to explain the data; therefore we avoid introducing it and assumed the antisymmetric condition for simplicity. Neglecting the $`\dot{E}\mathrm{\Omega }^3`$ term (equivalent to $`\dot{\nu }\nu ^2`$) we put forward the issue of whether particle acceleration and/or mass-loss follow an antisymmetric spin-down equation, like the Larmor equation does.

A monopolar spin-down term arises from any energy loss mechanism with a $`\dot{E}\mathrm{\Omega }^2`$ dependence, which can be expressed as

$$\dot{E}_{mon}=\beta \left(\frac{\mu _{}^2}{c^3}\mathrm{\Omega }^4\right)\left(\frac{c/\mathrm{\Omega }}{R_{}}\right)^2,$$
(16)
where $`\beta `$ is a dimensionless factor that can depend on time but not explicitly on $`\mathrm{\Omega }`$. It is an open question whether energy-loss processes like particle acceleration and pulsar winds can provide this sort of expression and we can only especulate about this point. Michel (michel (1969)) estimated the torque exerted on a neutron star by a pulsar wind as

$$T=(\pi \mathrm{\Phi }^2/4c)\mathrm{\Omega },$$
(17)
where $`\mathrm{\Phi }`$ is the magnetic flux carried away by the wind. Assuming $`\mathrm{\Phi }`$ can be constant, a pulsar wind would lead to a monopolar term in the energy loss equation. A detailed analysis of pulsar wind energetics has been given by Harding et al. (harding (1999)), who assumed an electric current flowing in the magnetosphere causing a magnetic torque on the star. The corresponding loss of rotational energy due to this torque is given by:

$$\dot{E}_{wind}=\frac{c}{3}B_{}^2R_{}^6\left(\frac{1}{R_{open}^2R_{LC}^2}\right),$$
(18)
where $`R_{LC}=c/\mathrm{\Omega }`$. $`R_{open}`$ is the distance from the center of the star to the magnetosphere region at which the magnetic field lines become open due to the pulsar wind. For $`R_{open}`$ independent of $`\mathrm{\Omega }`$, this implies $`\dot{E}_{wind}\mathrm{\Omega }^2`$, i.e., a monopolar-like term.

## 6 Conclusions

A multipole equation for pulsar evolution can be derived naturally as the Taylor expansion of the general equation $`\dot{\nu }=f(\nu ,t)`$. While it might give insight into short term evolution of pulsars, the time independent version of the multipole equation is inconsistent with observations, indicating that no equation of the form $`\dot{\nu }=f(\nu )`$ is able to fully model pulsar spin-down evolution. A time dependent multipole equation considering exponentially decaying magnetic-fields cannot reproduce the data without shortening the active life of pulsars below $`10^5`$ years. Introducing an inverse linear magnetic field decay with timescales of the order of $`10^4`$ years gives good agreement with observations and allows pulsars to remain active about $`10^7`$ years, in agreement with pulsar and supernova birth rates. The model is consistent with observations of the Vela pulsar only if its braking index is $`\begin{array}{c}>\\ \end{array}1.5`$. The different trajectories of this pulsar relative to those of the Crab, PSR B1509–58 and PSR B0540–69 suggest the existence two types of pulsars: those with relatively short decay timescales which evolve to periods $`P_{\mathrm{}}\begin{array}{c}<\\ \end{array}1\mathrm{s}`$ and Vela-like pulsars, with more persistent magnetic fields evolving to the larger period pulsars observed in the $`P\dot{P}`$ diagram. We believe the need to introduce a monopolar term in the spin-down equation is an indirect observational evidence that particle acceleration and/or mass-loss processes must follow $`\dot{E}\mathrm{\Omega }^2`$.

###### Acknowledgements.

César Alvarez thanks Dany Page for his valuable comments during the preparation of his Ph-D thesis. César Alvarez also thanks CONACyT for the Ph. D. support grant received (reference number 86472).


format: text

text:
# Lattice QCD meets experiment in hadron physics

## 1 Introduction

QCD is a key component of the Standard Model of particle physics. It gives us a rich spectrum of bound states of quarks and gluons whose properties are predictable from QCD if we can solve the theory. QCD is strongly coupled in this regime, however, and we need the non-perturbative techniques of lattice QCD to do this from first principles.

Most (but not all) questions which lattice QCD can address require calculations with a precision of a few percent to answer them. These include the spectrum of hadrons, their internal structure and decay rates. In particular, the hunt for internal inconsistencies in the Standard Model which could lead to new physics requires calculations of hadronic weak matrix elements to 2-3% to match the experimental errors that will become possible.

Figure 1 shows the recent status of the combined experimental and theoretical efforts to pin down the vertex of the Cabibbo-Kobayashi-Maskawa unitarity triangle ckmfitter . The different circular constraints on the vertex come from different decay rates of $`B`$ and $`K`$ mesons and are found by dividing experimental rates by theoretical results obtained in lattice QCD. The current constraints are strongly limited by current lattice QCD errors of around 20%.

Lattice QCD is hard and numerically very expensive. Recent progress us  has at last made precision calculations look possible and we will concentrate on that work and its implications in this review.

## 2 Lattice QCD calculations

Lattice QCD calculations proceed by the discretisation of a 4-d box of space-time into a lattice. The QCD Lagrangian is then discretised onto that lattice. The spacing between the points of the lattice, $`a`$, is $``$ 0.1fm in current calculations and the length of a side of the box is $`L`$ 3.0fm. Thus our simulations can cover energy scales from $``$ 2 GeV down to $``$ 100 MeV.

The Feyman Path Integral is evaluated numerically in a two-stage process. In the first stage sets of gluon fields (‘configurations’) are created which are representative ‘vacuum snapshots’. In the second stage, quarks are allowed to propagate on these background gluon field and hadron correlators are calculated. The dependence of the correlators on lattice time is exponential. From the exponent the masses of hadrons of a particular $`J^{PC}`$ can be extracted, and from the amplitude, simple matrix elements.

QCD as a theory has a number of unknown parameters, the overall dimensionful scale of QCD ($``$ the bare coupling constant) and the bare quark masses. To make predictions, these parameters must be fixed from experiment. In lattice QCD we do this by using one hadron mass for each parameter. The quantity which is equivalent to the overall scale of QCD on the lattice is the lattice spacing.

Lattice calculations are hard and time-consuming. Progress has occurred in the last thirty years through gains in computer power but also, often more importantly, through gains in calculational efficiency and physical understanding. One particular area which revolutionised the field from the mid-1980s was the understanding of the origin of discretisation errors and their removal by improving the lattice QCD Lagrangian. Discretisation errors appear whenever equations are discretised and solved numerically. They manifest themselves as a dependence of the physical result on the unphysical lattice spacing. In lattice QCD, as elsewhere, they are corrected by the adoption of a higher order discretisation scheme. The complication in a quantum field theory like QCD is the presence of radiative corrections to the coefficients in the higher order scheme which must be determined (using perturbation theory, for example howard ).

Physical understanding of heavy quark physics on the lattice has also made a huge difference to the feasibility of calculating matrix elements relevant to the $`B`$ factory programme on the lattice. The use of non-relativistic effective theories requires the lattice to handle only scales appropriate to the physics of the non-relativistic bound states and not the (large) scale associated with the $`b`$ quark mass. $`B`$ physics is now one of the areas where lattice QCD can make most impact.

One area which has remained problematic, but which this year’s results have addressed successfully, is the handling of light quarks on the lattice. In particular the problem is that of how to include the dynamical (sea) $`u/d/s`$ quark pairs that appear as a result of energy fluctuations in the vacuum. We can safely ignore $`b/c/t`$ quarks in the vacuum because they are so heavy, but we know that light quark pairs have significant effects, for example in screening the running of the coupling constant and in generating Zweig-allowed decay modes for unstable mesons.

Because quarks are fermions, they cannot be simulated directly on the computer, but must be ‘integrated out’ of the Feynman Path Integral. This leaves an QCD Lagrangian in terms of gluon fields which includes ln(det($`M`$)) where $`M`$ is an enormous ($`10^7\times 10^7`$) sparse matrix. The inclusion of dynamical quarks is then numerically very expensive, particularly as the quark mass is reduced towards the small values which we know the $`u`$ and $`d`$ quarks have.

Many calculations even today use the ‘quenched approximation’ in which the light quark pairs are ignored. Results then suffer from a systematic error of $`𝒪`$(20%). A serious problem with the quenched approximation is the lack of internal consistency which means that the results depend on the hadrons that were used to fix the parameters of QCD. This ambiguity plagues the lattice literature.

Other calculations have included 2 flavours of degenerate dynamical quarks, i.e. $`u`$ and $`d`$, but with masses 10-20$`\times `$ the physical ones. This approximation is better than the quenched approximation but large uncertainties remain because the $`s`$ quark is omitted. Results must also be extrapolated to the physical $`u/d`$ quark mass and chiral perturbation theory is a good tool for this. However, chiral perturbation theory only works well if the $`u/d`$ quark mass is light enough and, for errors at the few percent level, this means less than $`m_s/2`$. This has been impossible to achieve in most calculations.

New results this year us  have included $`u,d`$ and $`s`$ quarks in the vacuum, with light enough $`u/d`$ masses to perform accurate chiral extrapolations. The results use a new discretisation of the quark action - the numerically fast improved staggered formalism. This formalism is well-matched to the supercomputing power of a few Tflops that is currently achievable.

### 2.1 Improved staggered quarks

The starting point for the staggered quark formalism is the naive discretisation of the Dirac quark action onto a lattice. This action has good features: chiral symmetry and discretisation errors that appear only as the square and higher powers of the lattice spacing. The naive discretisation suffers from the notorious doubling problem, however. A single quark species on the lattice gives rise to $`16`$ quark species, or tastes, on a 4-d lattice. The additional tastes appear around the edges of the Brillioun zone, where $`p\pi /a`$, as copies of a $`p0`$ quark. This would not be a problem if there were no interaction between the different tastes since the quark action would then fall apart into 16 different pieces in an appropriate basis and we could take $`\mathrm{det}(M)^{(1/16)}`$ in simulations to give the effect of 1 quark flavour.

There is interaction between the different tastes, however. It is mediated by highly virtual gluons, with momenta around $`\pi /a`$. A quark of one taste can absorb or emit such a high momentum gluon and turn into a quark of another taste. The effects of this taste-changing interaction are quite severe for the naive action, giving rise to large discretisation errors (even though formally of $`𝒪`$$`(a^2)`$) and large perturbative renormalisation factors, e.g. for the quark mass, when translating from the lattice scheme to the continuum. The degeneracy in mass of mesons made from quarks of different taste is lost. This is most noticeable for the pions because there is a light Goldstone boson.

Because the taste-changing interaction is a high momentum one it can be understood in lattice perturbation theory. In particular, the effects can be significantly improved by suppressing the coupling of quarks to gluons of momenta $`\pi /a`$ in any direction. This is achieved by ‘smearing’ the gluon field in the action in a particular way gpl ; orginos , and can be thought of as part of the standard Symanzik programme for systematically removing discretisation errors from lattice actions.

It is simple to ‘stagger’ the naive action and its improved variant to remove an exact degeneracy of a factor of 4 in tastes which arises from the spin degree of freedom. This results in an action with 4 doublers which can be simulated on the lattice using $`\mathrm{det}(M)^{(1/4)}`$ per flavour. It is very fast numerically because there is only one spin degree of freedom per site and the eigenvalues of $`M`$ are well behaved. This is what has allowed the MILC collaboration to generate ensembles of configurations which include $`u`$, $`d`$, and $`s`$ quarks in the vacuum with much more realistic masses than before milc .

Some worries remain about potential non-locality in the action as the result of taking the fourth root. However, this causes no problem in perturbative QCD where a simple power series in $`x`$ is obtained for an action with $`\mathrm{det}(M)^x`$. Stringent non-perturbative tests are also then needed. Luckily these tests are possible in this formalism with present day computers because of its speed, and are exactly the calculations required to test (lattice) QCD. The results, shown in the next section, speak for themselves.

## 3 Recent results

The MILC collaboration have made sets of ensembles of gluon field configurations which include 2 degenerate light dynamical quarks ($`u,d`$) and 1 heavier one ($`s`$milc . Taking the $`u`$ and $`d`$ masses the same makes the lattice calculation much faster and leads to negligible errors in isospin-averaged quantities. The dynamical $`s`$ quark mass is chosen to be approximately correct based on earlier studies (in fact the subsequent analysis shows that it was slightly high and further ensembles are now being made with a lower value). The dynamical $`u`$ and $`d`$ quarks take a range of masses down as low as a sixth of the (real) $`m_s`$. The sets of ensembles divide into two different values of the lattice spacing, 0.13fm and 0.09fm, and the spatial lattice volume is $`(2.5\mathrm{fm})^3`$, reasonably large. Analysis of hadronic quantities on these ensembles has been done by the MILC and HPQCD collaborations us .

There are 5 bare parameters of QCD relevant to this analysis: $`\alpha _s,m_{u/d},m_s,m_c`$ and $`m_b`$. The lattice spacing takes the place of $`\alpha _s`$ in lattice QCD. It is important that these parameters are fixed using the masses of ‘gold-plated’ hadrons, i.e. hadrons which are well below their strong decay thresholds. Such hadrons are well-defined experimentally and theoretically and should be accurately calculable in lattice QCD. Using them to fix parameters will not then introduce unnecessary additional systematic errors into lattice results for other quantities. This has not always been done in past lattice calculations, particularly in the quenched approximation. It becomes an important issue when lattice QCD is to be used as a precision calculational tool. We use the radial excitation energy in the $`\mathrm{{\rm Y}}`$ system (i.e. the mass splitting between the $`\mathrm{{\rm Y}}^{}`$ and the $`\mathrm{{\rm Y}}`$) to fix the lattice spacing and $`m_\pi `$, $`m_K`$, $`m_{D_s}`$ and $`m_\mathrm{{\rm Y}}`$ to fix the quark masses.

We can then focus on the calculation of other gold-plated masses and decay constants. If QCD is correct and lattice QCD is to work it must reproduce the experimental results for these quantities precisely. Figure 2 shows that this indeed works for the unquenched calculations with $`u,d`$ and $`s`$ quarks in the vacuum. A range of gold-plated hadrons are chosen which range from decay constants for light hadrons through heavy-light masses to heavyonium. This tests QCD in different regimes in which the sources of systematic error are very different and stresses the point that QCD predicts a huge range of physics with a small set of parameters.

References gottlieb ; gray ; charm ; bernardlat03  give more details on the quantities shown in Figure 2. Here we will discuss some of these. Figure 3 shows the radial and orbital splittings in the $`b\overline{b}`$ ($`\mathrm{{\rm Y}}`$) system for the quenched approximation ($`n_f`$ = 0) and with the dynamical MILC configurations with 3 flavours of dynamical quarks. Our physical understanding of the $`\mathrm{{\rm Y}}`$ system is very good and there are a lot of gold-plated states well below decay thresholds, which makes it a valuable system for lattice QCD tests. We use the standard lattice NRQCD effective theory for the valence $`b`$ quarks, which takes advantage of the non-relativistic nature of the bound states. The lattice NRQCD action is accurate through $`v^4`$ where $`v`$ is the velocity of the $`b`$ quark in its bound state. This means that spin-independent splittings, such as radial and orbital excitations, are simulated through next-to-leading-order and should be accurate to $``$ 1%. Thus the test of QCD using these splittings is a very accurate one. The fine structure in the spectrum is only correct through leading-order at present and more work must be done to bring this to the same level and allow tests against, for example, the splittings between the different $`\chi _b`$ states gray .

The $`\mathrm{{\rm Y}}`$ system is a good one for looking at the effects of dynamical quarks because we do not expect it to be very sensitive to dynamical quark masses. The momentum transfer inside an $`\mathrm{{\rm Y}}`$ is larger than any of the $`u,d`$ or $`s`$ masses and so we expect the radial and orbital splittings to simply ‘count’ the number of dynamical quarks once we have reasonably light dynamical quark masses. The righthand plot of Figure 3 shows this to be true - the splittings are independent of the dynamical $`u/d`$ quark mass in the region we are working in (and therefore for the points plotted in the left hand figure of Figure 3 and in Figure 2).

The $`\pi `$ and $`K`$ decay constants are important light hadron matrix elements, related to the purely leptonic decay rate via a $`W`$, and experimentally well-known. These are very sensitive to light quark masses and require a well-controlled extrapolation in the $`u/d`$ quark mass and interpolation in the $`s`$ quark mass to get accurate results to compare to experiment. Chiral perturbation theory can be used to perform the $`u/d`$ quark mass extrapolation provided the masses used on the lattice are small enough for the expansion in powers of quark mass ($`m_\pi ^2/(1\mathrm{G}\mathrm{e}\mathrm{V}^2)`$) and its logarithms to work well. In practise this means that second order chiral perturbation theory should work at the 2% level for $`m_{u/d}<m_s/2`$. Note that the error is set by the largest quark mass used in the chiral fits, not the smallest.

Figure 4 shows the results and chiral extrapolation for the decay constants on the ensembles of MILC configurations with $`m_{u/d}^{sea}=m_s/2.3`$ and $`m_s/4.5`$ at a lattice spacing of 0.09 fm. The curves in the left plot show the chiral extrapolation using only results with $`m_{u/d}^{valence}<m_s/2`$. This extrapolation has to be corrected, using the lattice results, to interpolate to the physical $`s`$ quark mass for both sea and valence $`s`$ quarks. This then gives the results shown in Figure 2 which agree with experiment. The plot on the right shows what happens when the chiral extrapolation fit obtained in the left plot is evaluated for larger valence $`m_{u/d}`$. The $`f_\pi `$ results start to show clear disagreement for $`m_{u/d}>m_s/2`$, which makes the problem of performing accurate chiral extrapolations using results with $`m_{u/d}>m_s/2`$ obvious. Previous lattice calculations have been forced by computing cost to work only in this regime, with the added problem that the sea $`m_{u/d}`$ is also large aoki .

Another gold-plated hadron mass is that of the nucleon. A full chiral extrapolation of results for this on the MILC configurations has not yet been done. The left-hand plot of Figure 5 shows very encouraging signs that an answer in agreement with experiment will be found gottlieb . There is a clear sign of dependence on the lattice spacing, however, which will have to be taken into account. Combinations of baryon masses can be made which are relatively insensitive to $`u/d`$ quark masses and other effects and it is one of these, $`3m_\mathrm{\Xi }m_N`$, which is plotted in Figure 2.

It is important to realise that accurate lattice QCD results are not going to be obtainable in the near future for every hadronic quantity of interest. What these results show is that ‘gold-plated’ quantities should now work. Gold-plated hadrons are those well below decay threshold for strong decays. Unstable hadrons, or even those within 100 MeV or so of Zweig-allowed decay modes, have a strong coupling to their real or virtual decay channel which is not correctly simulated on the lattice. The problem is that, with the lattice volumes being used, the allowed non-zero momenta are typically greater than 400 MeV and this significantly distorts the decay channel contribution. Much larger simulations will be necessary to handle these hadrons.

Gold-plated hadrons include: $`\pi `$, $`K`$, $`D`$, $`D_s`$, $`J/\psi `$, $`\mathrm{{\rm Y}}`$, $`B`$, $`B_s`$, $`p`$, $`n`$, $`\mathrm{\Lambda }`$, $`\mathrm{\Omega }`$ etc. The following are not gold-plated: $`\rho `$, $`\varphi `$, $`D^{}`$, $`D_{sJ}`$, $`\mathrm{\Delta }`$, $`N^{}`$, pentaquarks, glueballs and hybrids in general. Lattice calculations will not get the masses right for non-gold-plated hadrons even when light dynamical quarks are included. This does not preclude lattice calculations giving useful qualitative results and insight but these points should be borne in mind for any quantitative comparison.

Figure 5 also shows the spectrum of $`D_s`$ states obtained on the dynamical MILC configurations fnalds . The valence $`c`$ quarks are simulated using an effective theory which, in a similar way to the $`\mathrm{{\rm Y}}`$ above, should be accurate for spin-independent splittings and not quite so accurate for fine structure in the spectrum. The hyperfine splitting between the $`D_s`$ and $`D_s^{}`$, for example is currently missing a radiative correction to the term in the action proportional to the spin coupling to the chromo-magnetic field. This is being calculated in lattice perturbation theory howard . Also shown are the scalar and axial vector orbital excitations compared to the recent experimental results for these mesons. The lattice calculation is giving a high result, albeit with large statistical errors at present. However, a high result is consistent with the fact that these mesons are not gold-plated and the lattice calculation does not currently include correctly the coupling to their decay modes.

Decay rates which can be accurately calculated for gold-plated hadrons are those in which there is at most one (gold-plated) hadron in the final state. This therefore includes leptonic and semi-leptonic decays and the mixing of neutral $`B`$ and $`K`$ mesons. Luckily there is a gold-plated decay mode available to extract each element (except $`V_{tb}`$) of the CKM matrix which mixes quark flavours under the weak interactions in the Standard Model:

$$\left(\begin{array}{ccc}𝐕_{\mathrm{𝐮𝐝}}& 𝐕_{\mathrm{𝐮𝐬}}& 𝐕_{\mathrm{𝐮𝐛}}\\ \pi l\nu & Kl\nu & B\pi l\nu \\ & K\pi l\nu & \\ 𝐕_{\mathrm{𝐜𝐝}}& 𝐕_{\mathrm{𝐜𝐬}}& 𝐕_{\mathrm{𝐜𝐛}}\\ Dl\nu & D_sl\nu & BDl\nu \\ D\pi l\nu & DKl\nu & \\ 𝐕_{\mathrm{𝐭𝐝}}& 𝐕_{\mathrm{𝐭𝐬}}& 𝐕_{\mathrm{𝐭𝐛}}\\ B_d|\overline{B}_d& B_s|\overline{B}_s& \end{array}\right)$$

As described earlier, the determination of the CKM elements and tests of the self-consistency of the CKM matrix are the current focus for the search for Beyond the Standard Model physics and lattice calculations of these decay rates will be a key factor in the precision with which this can be done.

First calculations on the dynamical MILC configurations have concentrated on the $`B`$ and $`B_s`$ leptonic decay rates wingatelat03 ; fnalds , because these are simplest. They are parameterised by the decay constants, $`f_B`$ and $`f_{B_s}`$, and these are an important component of the mixing rate for these mesons, which constrains $`V_{ts}`$ and $`V_{td}`$. Again one issue in extracting reliable lattice results for $`f_B`$ and $`f_{B_s}`$ is the chiral extrapolation in the $`u/d`$ quark mass. Figure 6 shows results on the MILC configurations for the ratio of $`f_{B_s}/f_{B_d}`$ plotted against the valence $`u/d`$ quark mass kronfeld . The data extend into the region $`m_{u/d}<m_s/2`$ which will allow an accurate chiral extrapolation for the first time. Although the statistical errors are currently rather large, it seems likely that the result for this ratio will be larger than previous estimates based on extrapolations from larger $`u/d`$ masses, and including only two flavours of dynamical quarks aoki . Further calculations of gold-plated matrix elements are in progress slstuff .

## 4 Conclusions

The impact of lattice QCD calculations has been hindered by the difficulty of including a realistic QCD vacuum. This has led to a level of systematic error far greater than the few percent needed to provide input to tests of the Standard Model, particularly those testing the CKM matrix at $`B`$ factories. New results this year look set to herald a brighter future in which accurate calculations, at least for gold-plated quantities, are available from the lattice at last.

We are grateful to PPARC, NSF, DoE and the EU for funding this work, and to all our collaborators on us  for many useful discussions. A version of this talk was also presented at the 2003 Lepton Photon Conference and is published in the proceedings of that meeting.


text:
# Minimalist’s Linux Cluster

## 1 Motivation

Constructing a Linux cluster using commodity PC’s and commodity networking hardware became quite easy and using such a cluster for a lattice QCD project is very popular. This increases the range of computing power available to those who have only moderate means. However, from our experience of using such a cluster, we found that there is a room for improvement in scaling up the current Linux cluster architecture : first, if many PC’s are just stacked on top of each other in rows, soon the cluster begins to occupy too large physical space. Secondly, not all the components in a standard PC is essential for a lattice simulation. By getting rid of unnecessary parts, one may reduce overall cost and electrical power requirement for each PC’s. Third, switched ethernet hub is usually used in a linux cluster and providing full bisection bandwidth using such switches is expensive. Building a cluster without a switch may be more scalable. On the other hand, for those who have limited resources like us, building everything (in particular, hardware components) from the scratch to alleviate the above problems is not sensible because it will probably take too long to develop such components. Thus we looked for a solution which does not require custom made hardware components and is re-usable in the future once developed so that the evolutionary upgrade does not introduce delays.

## 2 Architecture and Hardware

Each node is an extremely thin node which consists only of a Intel pentium IV 2.4GHz CPU, 512 Mbytes DDR SDRAM, a mother board, and 4 fast ethernet network interface cards(NIC). One of 4 NIC’s has a socket for EPROM or EEPROM for the bootcode. Table1 shows hardware components in the 36-node cluster and their costs. Except the chasis and the network cables, everything is off-the-shelf components and there is nothing special about them. The chasis is designed so that each crates accept any standard ATX-size mother board and an upgrade means just replacing mother boards with new one. The current chasis size is suitable for 64-node configuration and has room for additional 28 nodes. One PC with 360 GBytes hard disk serves as a front end server.

Thus, development effort for our thin node cluster is mostly involved with setting up necessary software environment : booting, OS, and MPI parallel programming. Since there is no permanent storage device on each nodes, booting is a little bit tricky and Linux operating system needs to be configured dynamically after the boot. Fortunately, there is a Linux solution, called “Linux Terminal Server Project”(LTSP), which is developed for the server-client situation similar to our case, a server booting up hosts of diskless client computers. In this scenario, instead of booting from the kernel image on a permanent media such as hard disk, floppy disk or flash memory device, an NIC which has a small size EPROM or EEPROM (for example 64 Kbytes) on the mother board does network booting. On power-up, this network card on the client node executes its bootcode and broadcast its IP request and its MAC address to the local network by use of Dynamic Host Configuration Protocol(DHCP). The server responds to this DHCP request and replies with the basic IP information such as client node IP address, netmask setting, root file directory and kernel image name depending on the client MAC addresses. With the reply from the server, the client node configures its TCP/IP and fetches kernel image from a host computer by Trivial File Transfer Protocol(TFTP). Once the kernel image is loaded on the node memory, the kernel starts executing and initializes the client node and set it up for normal operation. One may choose whether application softwares run on client nodes or runs on the server node.

The main difference between LTSP setup and ours lies on the assumed network topology. LTSP relies on the star network connection and our project adopts 2-dimensional mesh structure. In our case, each nodes once booted, must act as a DHCP server and a TFTP server to the next client node in constrast to the LTSP situation that central server controls the other client nodes. This booting process may progress in parallel to speed up : the front end server in our cluster starts booting processes on 6 nodes simultaneously and then these 6 nodes boot the next 6 nodes, etc (see Fig.1). After the booting procedure is completed, 2-D mesh routing is achieved by explicit ‘route’ command in a script called from Linux “init” script. ‘route’ assigns algorithmically one of four ethernet devices, eth0, eth1, eth2, and eth3 depending on the destination IP addresses and the logical node ID. Since the size of Linux routing table may grow upto 2048 elements by just changing kernel compile option, this kind of explicit routing work fine with a moderate size cluster. Ideally, one would like to have a distributed routing mechanism implemented on the kernel level but it is not part of the current Linux kernel. Linux distribution used on the cluster is ‘Wow Linux version 7.1’, which is equivalent to Red Hat Linux 7.1 and the kernel version is 2.4.9. The version of LTSP package which we modified for our need is 3.0.5. MPICH and LAM implementation of MPI parallel programming environment is available on the cluster.

## 3 Performance and Discussion

Fig. 2 shows the code performance of hybrid molecular dynamics simulation of two staggered quark flavor with $`m_qa=0.01`$ on a $`8^3\times 512`$ lattice (the single node benchmark is for $`8^3\times 32`$ lattice). One-dimensional ring ($`N_t=512`$ is distributed over the nodes) layout of lattice sites is used for the code and the code is not yet optimal for the 2-D mesh structure of the cluster. However, the code performance scales up nicely between 1 to 8 nodes. Sustained speed is about 2.25 GFLOPS on 8 node and is 11% of the theoretical peak speed. Thus, our cluster achieved $`0.5`$ MFLOPS/$`\$`$ with a straight FORTRAN code with no assembly language subroutine. We find that using more than 8 nodes with the current full QCD test code quickly degrades cluster performance due to non-optimal communication pattern of the test code.

A conventional PC with Intel Pentium IV CPU is ordinarily equiped with 350W power supply. Since we put 180W power supply for each node and the cluster operates fine with this condition, the overall power requirement is successfully reduced to a half of usual Linux cluster. Also, the footprint of our cluster is $`200\times 91\times 75`$ cm, which is considerably smaller than that of stacking 36 PC’s. The physical dimension of the full cluster (64 node) will be even more beneficial since the same chasis will be used. Saved node cost would be $`100\$(15\%)`$ from doing without a hard disk. The whole construction is reusable as we planned since the mother board size is the only factor which needs to considered in an upgrade and the standard ATX size of mother board will stay with us for a while.

Global MPI operations such as “MPI\_ALLREDUCE” involves multiple hops in our cluster. Since LAM or MPICH relies on TCP/IP and each hops contributes to software and hardware latency in message passing, transversing many nodes reduces the efficiency of a program in our cluster. However, since the software latency is larger than the hardware latency, multiple hop will be less severe problem when user space devices such as Infiniband becomes cheaply available.


text:
# Parsec - scale Herbig-Haro Outflows from Intermediate Mass Stars

## 1 Introduction

HH objects are the shock - excited nebulous tracers of outflows from pre - main sequence stars. In many cases these outflows are collimated in the form of highly supersonic jets the existence of which appear to be intrinsicly linked to accretion by the underlying young stellar object (YSO) (Hartigan, Edwards, & Ghandour, 1995). Most known optical jets have low mass ($``$ 1M ) sources: either the embedded (IRAS Class I) counterparts of classical T Tauri stars or classical T Tauri stars themselves (Reipurth, Bally, & Devine, 1997).

Turning to higher mass YSOs ($`\mathtt{>}`$ 10M ), such as those driving the Orion OMC 1 or Cepheus A outflows, one sees a very different picture (Allen & Burton, 1993; O’Dell et al., 1997; Hartigan, Morse, & Bally, 2000). Their outflows, although highly energetic, often appear poorly collimated and more chaotic (Reipurth & Bally, 2001) than their low mass counterparts. This transformation begs two obvious inter-related questions: at what point does the transition occur and is it a smooth function of mass? To answer these questions one must examine outflows from YSOs of mass greater than $``$2 M .

Optical outflows have been observed from a number of intermediate - mass (2 M $``$ M $``$ 10 M ) YSOs, for example R Mon, LkH$`\alpha `$ 234 , and AFGL 4029 (Ray et al., 1990 and references therein). Such stars, where optically visible, are known as Herbig Ae/Be stars (HAEBES) although their embedded counterparts have also been seen. However optical outflows from these YSOs are rare and there are a number of reasons for this. The initial mass function favours the production of low mass stars and therefore intermediate mass YSOs tend to be found at relatively large distances. More massive stars also have a faster evolutionary timescale – i.e. they evolve more quickly than low - mass YSOs and so their outflow phase is shorter. This would also make their outflows more difficult to detect. Another contributing factor could be that massive stars tend to be more obscured : intermediate and massive stars tend to be surrounded by large amounts of circumstellar gas and dust, making it harder to find an outflow at visual wavelengths. Finally there may well have been a historical bias towards studying outflows from low - mass stars. The situation however has changed in recent years as more and more studies focus on higher mass YSOs.

With these ideas in mind, we have investigated the occurrence of large-scale outflows from intermediate mass stars. By large scale we mean those stretching several parsecs e.g. the PV Cephei outflow at 2.6 pc (Gomez, Kenyon, & Whitney, 1997; Reipurth, Bally, & Devine, 1997), the HH 80/81 5.3 pc outflow (Marti, Rodriguez, & Reipurth, 1993) and the HH 354 outflow at 2.4 pc (Reipurth, Bally, & Devine, 1997). We emphasise that these outflows have vastly longer associated timescales than those of “traditionally” observed flows. A small - scale HH jet close to its source has a dynamical timescale of only a few hundred years, whereas the HH objects in these parsec - scale flows trace mass ejection over tens of thousands of years. They are, in effect, fossil records of the mass - loss histories of their parent star.

Newly detected parsec - scale outflows around five intermediate - mass young stars are discussed here (see also Table 1 and Table 2). Of these, LkH$`\alpha `$ 198, LkH$`\alpha `$ 233 and LkH$`\alpha `$ 234 are of spectral type A and 1548C27 is A7 - F0. All of these stars were known to possess small-scale optical outflows. The one optically invisible YSO in our sample, IRAS 19395+2313, was not previously known to drive any outflow.

Section 2 describes how we made our observations and in Section 3 we present our newly discovered large-scale flows. The implications of our findings are discussed in Section 4 and our conclusions are drawn in Section 5.

## 2 Observations

To carry out our survey we used the Wide Field Camera (WFC) at the prime focus of the 2.5m Isaac Newton Telescope at El Observatorio del Roque de los Muchachos (La Palma, Canary Islands). The WFC consists of four thin - coated EEV CCDs each with 2048 x 4100 15$`\mu `$m2 pixels. One pixel projects to 0$`\stackrel{}{.}`$33 on the sky. Three of the CCDs are positioned from north to south with their long axes adjoining. The fourth is attached to the west to form a square mosaic (34$`\stackrel{}{.}`$2 wide) with its northwestern corner missing.

Our images were taken on nights between the 13th and the 21st of July 1998. Seeing was moderate at 1$`\stackrel{}{.}`$15 – 1$`\stackrel{}{.}`$5 as measured from the images. HH objects were identified using a number of narrowband emission line filters: H$`\alpha `$($`\lambda _c`$ = 6568Å, $`\mathrm{\Delta }\lambda `$(FWHM) = 95Å), \[SII\]($`\lambda _c`$ = 6725Å, $`\mathrm{\Delta }\lambda `$(FWHM) = 80Å) and \[OIII\]($`\lambda _c`$ = 5008Å, $`\mathrm{\Delta }\lambda `$ (FWHM) = 100Å). To distinguish HH emission from reflection nebulosity, we also took broadband images in V and I. Exposure times for the narrowband and broadband images were typically 30 and 10 minutes respectively. The data was reduced using standard IRAF reduction procedures.

## 3 Results for Individual Regions

### 3.1 LkH$`\alpha `$198 & V376 Cas

LkH$`\alpha `$ 198 and its nearby companion, V376 Cas, are both Herbig Ae stars (Herbig, 1960) located in the small dark cloud L1265, at a distance of 600 pc (Chavarria-K., 1985). An asymmetrical, bipolar molecular outflow in this region was noted by Canto et al. (1984). Strom et al. (1986) subsequently found the first optical outflow tracer, HH 161, a bright HH object some 12″ from LkH$`\alpha `$ 198 at a position angle (P.A.) of 100°. Further observations by Goodrich (1993) yielded another HH knot 81″ away at a P.A. of 153°. This object was rediscovered by Aspin & Reipurth (2000) who refer to it as HH 461.

The discovery, however, of LkH$`\alpha `$ 198 B, a deeply embedded companion to LkH$`\alpha `$ 198, by Lagage et al. (1993) raised the question of which of these two stars is the primary outflow source in this region. In their study Corcoran, Ray, & Bastien (1995) (hereafter referred to as CRB) concluded there are two separate outflows with their origin in the vicinity of LkH$`\alpha `$ 198: one driven by LkH$`\alpha `$ 198 itself and the other by LkH$`\alpha `$ 198 B. Their observations of HH 161 revealed a tail pointing back towards LkH$`\alpha `$ 198 B and they also discovered a suspected bow shock 39″ southeast of this source. The bow shock (their knot B), HH 161 and its tail are all aligned and appear to constitute a one-sided outflow from LkH$`\alpha `$ 198 B. To date no counterflow has been seen. CRB also found a number of faint HH knots (HH 164 C, D and E – see inset in Fig. 1) with the same P.A. from LkH$`\alpha `$ 198 as HH 461. Thus HH 164 C, D and E and HH 461 appears to delineate an outflow from LkH$`\alpha `$ 198. To the north they also discovered a faint knot (HH 164 F) which may be a tracer of the counterflow from LkH$`\alpha `$ 198. Finally to the east of V376 Cas a number of bright HH emission knots, HH 162, were also found by CRB. These knots were seen again in our images, although their origin remains unclear.

Two of the three HH complexes found here are extensions of the outflow from LkH$`\alpha `$ 198 discovered by CRB. HH 801 has the same P.A. of 340° with respect to LkH$`\alpha `$ 198 as HH 164 F and is therefore almost certainly part of the same flow. It appears to be an asymmetrical bow shock of which we primarily see the western wing (B in Fig. 2) some 45″ in length. There are also a number of more easterly knots, C - F. Knot A may also be part of the western wing.

In the diametrically opposite direction from HH 801 through LkH$`\alpha `$ 198 we find HH 802. It consists of a number of features, A – I (Fig. 3) which at first sight look somewhat chaotic. Feature E, however, like its counterpart in HH 801, could be the western wing of an asymmetrical bow shock. The total length of HH 802 is some 2. It has a P.A. of 160° with respect to LkH$`\alpha `$ 198 and it is aligned with knots C, D and E (CRB) of HH 164 and HH 461 i.e. it is the counterflow of HH 801 and HH 164 Knot F. The furthest knot in HH 802 (Knot I) is at a distance of 8$`\stackrel{}{.}`$16 (1.4 pc) from LkH$`\alpha `$ 198, implying the total projected extent of the flow, from HH 801 to HH 802, is some 2.3 pc.

HH 800 to the northwest of HH 801 is unlikely to be part of the HH 801 – HH 802 outflow unless the outflow direction has changed abruptly. Although changes in flow direction have been observed in other parsec-scale outflows (Reipurth & Bally, 2001), they tend to be more gradual. Moreover HH 802 is even further from LkH$`\alpha `$ 198 than HH 800 (at least in projection) and the flow associated with the former appears to have maintained a constant outflow direction. We should also add that HH 800 is probably not part of the counterflow from LkH$`\alpha `$ 198 B as its outflow is at a P.A. of 132°, while HH 800 is at $``$330° with respect to this source. We cannot however exclude the possibility that the outflow axis may have swung through 18°. Finally either of the two optically invisible IRAS sources in the region (Fig. 1) could be its driving source. IRAS completeness in this region is of the order of 5L . Proper motion studies would clearly help to identify its origin.

### 3.2 1548C27

The cometary-shaped reflection nebula 1548C27, and its associated H$`\alpha `$ emission line jet (HH 165) were first noted by Craine, Byard, & Boeshaar (1981). The optical jet, at a P.A. of 54°, was later confirmed by Mundt et al. (1984). A low - mass, poorly collimated molecular outflow was also observed in the region by Dent & Aspin (1992).

Near-infrared photometry in the immediate vicinity of 1548C27 by Vilchez et al. (1989) yielded two sources. One of these appears to be a foreground star but the other, IRS 1, is located near the apex of the nebula, and they suggest this to be the driving source of HH 165 (see Fig. 4). The IRAS PSC (point source catalogue) shows IRAS 19407+2316 to be located close to, but not coincident with, IRS 1. Vilchez et al. (1989) however conclude that both near - and far - infrared sources are the same object, which, for convenience, we will refer to here as IRS 1.

The optical jet (HH 165) is very narrow, with a width of 2′′ \- 3′′ and its length is estimated to be $``$45′′. There is a gap of $``$13′′ between the source and the jet and two bright knots are visible at 23′′ and 29′′ from the source (Mundt et al., 1984). The kinematic distance of 1548C27 is 2.4 kpc Dent & Aspin (1992) which implies a luminosity for IRS 1 of approximately 580 L.

Another star S2, 10″ northeast of 1548C27, was found by Scarrott, Rolph, & Tadhunter (1991). This star illuminates the nebula along with IRS 1. Scarrott, Rolph, & Tadhunter (1991) suggest HH 165 curves towards this star, implying that it is the driving source, however we find no evidence in our images to support this idea.

HH 365 to the northeast of 1548C27 was briefly referred to by Alten et al. (1997) as being bow-shaped and possibly associated with HH 165. This object was mentioned in their paper but no image of it was included. Its morphology is clearly seen here in Figs. 4 and 5. HH 365 is 8$`\stackrel{}{.}`$13 (5.7 pc) from IRS 1 at a P.A. of 46° i.e. close to that of the HH 165 jet. From our images it appears to consist of two bright regions, Knot A and Feature B that extends to the northwest some 8′′ (Fig. 5). The overall shape of Feature B is suggestive of an asymmetrical bow, the axis of which points roughly back towards 1548C27.

Our survey also revealed a number of possible new HH objects (see Fig. 4) although several are very faint. Moreover, it is unclear whether Knots A, B and C (Fig. 6), for example, to the northwest of HH 165 are HH objects, as there is a lot of contaminating HII nebulosity in the region. However the fact that these lie on the path between HH 165 and HH 365 would suggest they might be. Further study is necessary, however, to determine their nature and for this reason we will desist from assigning them HH numbers. In any event it seems likely that HH 165 and HH 365 are part of the same outflow from IRS 1 and that Knots A, B and C may also be part of this flow.

A counterjet from IRS 1 was found recently in the near-infrared (Whelan, private communication). It lies along the same line as HH 165 at a P.A. of 234° and extends for at least 5″. A number of faint \[FeII\] emitting knots were also seen beyond the counterjet. Neither the counterjet nor any of these knots are observed in our optical images, presumably because of extinction.

HH 803 (Fig. 7), 2$`\stackrel{}{.}`$63 (1.85 pc) southeast of IRS 1 at a P.A. of 223°, has a very interesting morphology. It appears to be bow-shaped but convex towards IRS 1. It is 50″ in width and contains a 13″ “jet - like” feature bisecting the bow. The “ring” at the southern edge of the bow and the diffuse appearance of Knot A to the north add to the complexity of this object. Note that the jet - like feature does not quite point back towards IRS 1. We should also add that apart from IRS 1, no other IRAS sources were found in the region strengthening the possibility that HH 803 is driven by IRS 1.

If we include Knots A, B and C as part of the HH 165/HH 365/HH 803 outflow then its overall appearance suggests that it may be slowly precessing with shifts in the outflow direction of at least 10°. The sense of precession (i.e. point-like symmetry through IRS 1) to the northeast is also consistent with the position of HH 803 to the southwest.

### 3.3 The IRAS 19395+2313 Region

A number of HH objects and possible sources were found in this region (see Fig. 8), which lies approximately 18 west of 1548C27 and is also in the vicinity of the young open cluster NGC 6823. It is highly unlikely that any of these newly-discovered HH objects are part of the 1548C27 IRS 1 outflow although we will assume they are at the same distance, i.e. 2.4 kpc.

Two of the newly discovered HH objects are found close together - HH 806 is 30″ west of HH 807 (Fig. 8). The region between them coincides with a gap in the CCD mosaic although a cursory inspection of the Palomar Sky Survey Red (E) plate shows there is a star in the “gap”, ALS 10422 or IRAS 19399+2312 (at 19h42m05.5s +23°1859′′, J2000) which, at first sight, might be a possible HH driving source. This object is however an AGB star (Parthasarathy, Vijapurkar, & Drilling, 2000) so we can disregard it. We also see on the Palomar Sky Survey red (E) plate, a conical nebula $``$36″ southeast of HH 806 (at 19h42m02.12s +23°1930′′, J2000) which may be associated with the driving source of this object. In fact HH 806 lies along the major axis of this conical nebula. Although this nebula is not seen in Fig. 8, its position is marked. There is no obvious driving source for HH 807.

HH 805 has an interesting morphology (Fig. 9) and IRAS 19395+2313 seems an obvious candidate to be its driving source given its position. We have estimated the luminosity of IRAS 19395+2313 to be $``$ 320 L . The angular extent of HH 805 is approximately 45″ (0.5 pc) and it can be seen from Fig. 9 that this outflow is poorly collimated. Morphologically it appears to have a knotty ring-like structure and is reminiscent of the HH complex associated with V380 Ori (Corcoran & Ray, 1995). HH 804 (Fig. 8) may be part of the counterflow from IRAS 19395+2313 although we emphasise that this association is highly uncertain. It is at a P.A. of 80° with respect to the latter. HH 804 is at a distance of 6$`\stackrel{}{.}`$3 from IRAS 19395+2313, implying the total projected extent of the flow, assuming HH 804 is part of it, is $``$5 pc.

### 3.4 LkH$`\alpha `$ 233

LkH$`\alpha `$ 233 is an A5e - type pre - main sequence star (Corcoran & Ray, 1997) at a distance of 880 pc and is associated with a bipolar nebula that is approximately 0.1 pc in size (Calvet & Cohen, 1978; Staude & Elsasser, 1993). The nebula has a distinct X - like morphology with bright reflection limbs at 50° / 90°and 230° / 270°.

The discovery by Leinert, Haas, & Weitzel (1993) of a light scattering “halo” $``$1′′ in size around LkH$`\alpha `$ 233 led them to suggest that the star is highly embedded and optically visible largely through scattered light. Corcoran & Ray (1998) discovered a bipolar jet and counterjet (HH 398) spectroscopically at P.A.s of approximately 245° and 65° that bisect the X - shaped nebula. In their spectrograms, the redshifted counterjet is seen to begin $`0^{\prime \prime }.7`$ from the centre of the stellar continuum, whereas the blueshifted jet can be traced right back to the continuum peak. The jet, and counterjet, can be seen in our continuum subtracted image (Fig. 12) along with a number of HH objects along the outflow direction.

Evidence for the presence of a large “polarisation disk” centred on LkH$`\alpha `$ 233 with a radius of about 15000 AU was found by Aspin, McCaughrean, & McLean (1985). The position angle of this disk is about 155°, which places it perpendicular to the observed optical outflow. The existence of this “polarisation disk” combined with the fact that the counterjet is not seen spectroscopically close to the star leads to the conclusion that there is a circumstellar disk present that obscures the beginning of the receding flow.

Close to LkH$`\alpha `$ 233 the \[SII\] emission can be resolved into two velocity components (Corcoran & Ray, 1998). The high velocity component is identified with the optical jet while the low velocity component, which extends to less than 2′′ from the star, may be modelled as a rotationally broadened disk wind (see Kwan & Tademaru, 1995).

This survey revealed a number of previously unknown HH objects in the vicinity of LkH$`\alpha `$ 233 . It is possible that not all of these objects can be attributed to LkH$`\alpha `$ 233 and the positions of a nearby IRAS source, IRAS 22317+4024, complicates our analysis (see Fig. 10). Candidate driving sources for all the new HH objects are suggested here.

Continuum subtracted images of the nebula surrounding LkH$`\alpha `$ 233 reveal a number of HH features which are not seen, at least so clearly, in the \[SII\] image alone. Fig. 12 shows the first optical images of an $``$7′′ jet emerging from the LkH$`\alpha `$ 233 nebula at a P.A. of $``$248°, which is relatively close to the P.A. of the jet as spectroscopically determined by Corcoran & Ray (1998). A possible counterjet to the northeast of LkH$`\alpha `$ 233 is also seen in this image. But it is difficult to determine whether this is actually a counterjet or simply a residual of the continuum subtraction process. There are two other emission knots to the southwest of the source. The first of these, HH 808, is situated close to a diffraction spike from a bright star to the west of LkH$`\alpha `$ 233 . HH 808 is 1$`\stackrel{}{.}`$05 from the source at a P.A. of 250°. The second knot, HH 809, is 2$`\stackrel{}{.}`$2 from the source at a P.A. of 241°. Two other objects were seen to the northeast of LkH$`\alpha `$ 233 at a distance of $``$2, the first has a P.A. of 63° and the second is at 65°. These objects are outside the area shown in Fig. 12 and it is unclear at present whether these are HH objects.

Also to the northeast of LkH$`\alpha `$ 233 , we discovered HH 814 (Fig. 10) at a distance of $`5^{}.18`$ (1.3 pc). The morphology of this object may be studied more clearly in the continuum subtracted image inset in Fig. 10. HH 814 appears to be a broad bow shock facing back towards LkH$`\alpha `$ 233 at a P.A. of 47° with respect to LkH$`\alpha `$ 233 . Corcoran & Ray (1998) determined a P.A. of approximately 65° for the counterjet, suggesting that if HH 814 is part of the same flow, its direction has changed by $``$20°. Note however that the P.A. determined by Corcoran & Ray (1998) is very crude as it was deduced by slit sampling at various P.A.’s.

A number of other objects were discovered to the southwest of LkH$`\alpha `$ 233 (see Fig. 11). HH 810, HH 811 and HH 812, are at 4$`\stackrel{}{.}`$5 (1.2 pc), 6 (1.5 pc) and 6$`\stackrel{}{.}`$5 (1.7 pc) respectively from LkH$`\alpha `$ 233 , all at a P.A. of 236°. HH 813, at a distance of 7 (1.8 pc) has a P.A. of 247° with respect to LkH$`\alpha `$ 233 . Considering the possibility that outflows from higher mass stars may not be as collimated as those from low mass stars, HH 810 - HH 814 could be the optical tracers of the edges of a moderately collimated flow driven by LkH$`\alpha `$ 233 . The axis of this outflow with respect to LkH$`\alpha `$ 233 is $``$62° (marked on Fig. 10), which agrees well with the estimate of the P.A. of the jet closer to the source. It is also possible that HH 813 is a bow shock facing back towards IRAS 22317+4024 (Fig. 11). Proper motion studies could conclusively determine whether this is the case.

### 3.5 The NGC 7129 region

A large number of YSOs are known in the NGC 7129 cluster which lies at a distance of 1.25 kpc (Shevchenko & Yakubov, 1989). Aside from optically visible young stars such as LkH$`\alpha `$ 234 , there are many embedded ones. For example, a 160 $`\mu `$m survey by Bechis et al. (1978) revealed two far infrared sources and although one of them is spatially coincident with LkH$`\alpha `$ 234 (FIRS 1), the other is 3 further south (FIRS 2) and optically invisible. Additional infrared (Weintraub, Kastner, & Mahesh, 1994; Cabrit, et al., 1997) and sub - millimetre sources (Fuente et al., 2001) are also known.

Associated with these sources one finds the usual host of phenomena typical of star formation including reflection nebulae (Bertout, 1987), molecular outflows (Edwards & Snell, 1983; Bertout, 1987; Mitchell & Matthews, 1994), HH objects (Ray, 1987; Miranda, Eiroa, & Gomez de Castro, 1993), HH jets (Ray, 1987; Ray et al., 1990; Cabrit, et al., 1997) and shocked H2 flows (Eislöffel, 2000).

Looking at the distribution of previously known HH objects in the region, and the new ones discovered here (see below), one gets the impression (see Figs. 13 and 14) that the primary outflow axis, or axes, is roughly in a northeast to southwest direction and centred on the cluster core. Caution however is necessary. Proper motions studies have shown that some HH objects like GGD 32 and HH 103 are not moving to the southwest, as one might suspect, but instead to the west (Ray et al., 1990). Moreover others, like GGD 34 (Gómez de Castro & Robles, 1999) and possibly GGD 33 (Cohen & Schwartz, 1983; Goodrich, 1986), have their own sources outside the cluster core.

The HH 815 complex (Fig. 15) is over 1 in size and is at a distance of 11 (4 pc) from LkH$`\alpha `$ 234 . The three emission regions in HH 815 (A, B and C) appear to form the edges of a large bow shock that is concave towards the cluster core. HH 816 (Fig. 16) may be another bow trailed by a series of faint HH knots i.e. HH 817 - HH 820. HH 821, $``$15′′ north of GGD 35, is aligned with HH 819 and HH 816 and so could form part of the same flow. We note also that HH 105, is on the same axis from the cluster core and that no source for HH 105, at least in its immediate vicinity, is known.

The most distant HH object to the southwest of the cluster discovered by us is HH 822 (Fig. 17). It is at a P.A. of 238° with respect to LkH$`\alpha `$ 234 and the morphology of its brightest component (HH 822 A), i.e. a bow concave towards the cluster, suggests it is part of an outflow that originated there. An additional knot close to HH 103 was also found and we shall refer to it as HH 103 A (Fig. 13). Given its location, it seems likely it is part of the same flow that drives HH 822. Moreover faint emission can be seen linking HH 103 A to HH 822 A reinforcing this conclusion. Finally we add that HH 822 B is at a distance of 10$`\stackrel{}{.}`$7 (3.9 pc) from LkH$`\alpha `$ 234 .

The large number of sources in this region makes it extremely difficult, without detailed kinematic studies, to determine the origin of individual HH complexes. As previously mentioned there are a number of low mass YSOs, like the one that drives GGD 34, present and this complicates our analysis. That said, it seems likely, purely on morphological grounds as well as location, that many of the newly discovered HH objects are driven by a source(s) in the cluster core. Because of their large distance from the core, however, it may prove impossible, even with good proper motion data, to determine their precise origin. Several possible candidates exist including LkH$`\alpha `$ 234 , IRS 6 and FIRS 1 - MM1. If we draw an imaginary axis through the core at a P.A. of 60° / 240° it roughly delineates the region where most of the new HH objects, HH 816 to HH 822, are located. Assuming we are dealing with one outflow here, that originates in the core, then its overall angular size is 21$`\stackrel{}{.}`$8 i.e. $``$8 pc in projected length.

There are three other HH objects which are not situated along the major axis marked in Fig. 13. HH 824 and HH 825 are located on either side of IRAS 21416+6556 at 143° / 323° with respect to this source, suggesting a possible bipolar outflow driven by IRAS 21416+6556. These objects are 2′ and 51″ respectively from IRAS 21416+6556. The driving source of HH 823 is unclear as is its association with any of the other outflows in this region.

## 4 Discussion

### 4.1 Overall lengths

It has been stated, albeit somewhat tongue-in-cheek, that the apparent length of optical outflows from YSOs used to be a function of CCD chip size! Early chips were small and sampled only a small angular patch of the sky. This, in combination with the episodic nature of the flows themselves, conspired to suggest outflow lengths measured in tens of thousands of AU rather than parsecs. In a number of cases however CCD mosaicing (Ray, 1987) did hint that some flows were at least in the parsec league. Large format CCDs in focal - plane mosaics can cover fields of view larger than 30′ which, at a distance of 1kpc, corresponds to more than 8.7pc. Flows can therefore be detected to beyond the peripheries of their parent cloud.

Table 2 lists parameters such as distance, source luminosity, outflow length, cloud size and degree of collimation for both the intermediate - mass sources discussed here and several well-studied outflows powered by both intermediate - and low - mass YSOs. In some regions, the projected lengths of the outflows are similar to the sizes of the clouds from which they emerge. This correlation of length scales is to be expected and there should also be a tendency for shocks to be seen near the cloud edges where extinction by dust is minimal.

If we assume an average tangential velocity for outflows of 100 kms-1, in 105 yrs material originally at the source will be transported $``$10 pc i.e. typically to the edge of the clouds we are studying. Using the evolutionary tracks of Palla & Stahler (1993), we see that such a period corresponds to about 10% of the time an average Herbig Ae/Be star spends in the pre-main sequence phase. It follows that only in the case of the more massive YSOs, i.e. those with the shortest evolutionary timescales, do HH outflows represent a “fossil record” of activity over the totality of the outflow phase.

### 4.2 Morphology

HH outflows are episodic: they are clearly not continuous phenomena. For the most part, the shocks we see are generated by supersonic jet material ramming into previously ejected slower gas. This process produces a series of “working surfaces”, radiative shock systems that fade with time and therefore with distance from their source. Only the strongest shocks survive to produce dramatic, often chaotic, structures on parsec scales. It has even been suggested that the FU Orionis phenomena may signal the dramatic change in output needed at the source to produce such a shock (Reipurth, 1989). Thus the spacing between HH objects increases with distance from the source at least amongst low mass YSOs (Reipurth & Bally, 2001). Such a trend is also visible here amongst the flows from intermediate mass YSOs such as LkH$`\alpha `$ 198 and 1548C27. Another phenomenon that occurs with HH outflows from low mass YSOs is that the shock structures become larger and apparently more chaotic with distance. Again this is something which is replicated in their intermediate mass YSO counterparts.

A phenomenon that is found in parsec-scale outflows from lower mass YSOs is that the flow often exhibits “S” or “C” shape symmetry (Reipurth, 1989) possibly as a result of jet precession or source motion through the parent cloud respectively. In our small sample, we do not find any clear-cut examples of either although, as previously mentioned, the parsec-scale outflow from 1548C27 IRS 1 may be “S-shaped”.

### 4.3 Collimation

As alluded to in the Introduction, one of the most striking differences between parsec-scale outflows from low and massive YSOs is the relative lack of collimation seen in the latter (Shepherd, Churchwell, & Wilner, 1997; Hunter, Phillips, & Menten, 1997; Shepherd, Watson, Sargent, & Churchwell, 1998). This point is well illustrated by the archetypal example of an outflow from a high luminosity YSO: the Orion IRC2 flow (Allen & Burton, 1993). Its opening angle is $``$ 90° (Burton & Allen, 1994) which is a sharp contrast to the angles (typically a few degress) seen in outflows from low - mass YSOs.

Table 2 list “final” opening angles for our small sample. Values are determined by dividing the width of the most distant HH object by its projected source separation. Note that as we are using projected separations the quoted values must be upper limits, however our observations only reveal the brightest portions of the shocks and fainter outer parts may be missed, so the opening angles may actually be greater than what is optically observed. Observed opening angles range from 0.9° to 12° suggesting a degree of collimation comparable to that seen in outflows from low - mass YSOs. This also suggests that the transition from well - collimated outflows to poorly - collimated outflows occurs at higher masses than the sources observed here.

### 4.4 The Frequency of Blow - outs

The true size of an outflow, in comparison to that of its parent cloud, is an important factor in determining whether the outflow’s energy and momentum is transported into the ISM or remains within the cloud itself. A cursory examination of our data shows a clear tendency for HH objects to lie close to the edges of the parent cloud or at least close to the edges of clumps. As already mentioned, such an effect is to be expected considering the low extinction near cloud peripheries and the very low ISM densities beyond the cloud boundaries.

More importantly, it is clear that actual outflow timescales are very long in comparison to apparent dynamical ones. Here dynamical timescales are derived by dividing the projected length of an outflow by its estimated tangential velocity. This, and the observation of HH objects near cloud boundaries, immediately suggests that most, if not all, of the flows studied here have blown-out of their parent cloud.

### 4.5 Are outflows a source of cloud turbulence?

There is plenty of evidence to suggest molecular clouds are turbulent (Ward - Thompson 2002 and references therein) and that the pressure generated by this turbulence is sufficient to prevent clouds from collapsing under their own gravity. It has been shown however that cloud turbulence, even in the presence of magnetic fields, decays too quickly compared to typical cloud lifetimes (Stone, Ostriker, & Gammie, 1998) and must therefore be somehow replenished. Could outflows be the primary source of turbulence in molecular clouds?

The parsec-scale HH outflows imaged here appear largely well-collimated and therefore one might think they could affect only a narrow cone of ambient material. One has to remember however that these flows are supersonic and that we are viewing only the most highly collimated outflow component. The same flows “imaged” in the CO J = 1$``$0 line would normally appear much less collimated, especially at low velocities. That these flows affect cloud structure on parsec scales is evident from features such as the large CO “cavity” in NGC 7129. No doubt they also affect cloud dynamics. Arce & Goodman (2001) for example has found that molecular outflows associated with parsec-scale HH flows can possess kinetic energies comparable to the turbulent and gravitational binding energies of their parent clouds. This proves however only that they are a potentially important source of turbulence. Unfortunately we do not understand the coupling between outflows and their ambient medium well enough to be sure. Numerical simulations (Downes & Cabrit, 2003) are helping to address this problem but further studies are required.

## 5 Conclusions

We have investigated the occurrence of parsec-scale outflows from intermediate - mass YSOs. As is the case with lower mass YSOs, such flows appear to be common and we report the discovery of four here with well defined sources. These include LkH$`\alpha `$ 198, 1548C27 IRS 1, LkH$`\alpha `$ 233 and IRAS 19395+2313. All of these YSOs, with the exception of the last, were previously known to have small-scale outflow activity.

The region surrounding the Herbig Ae star LkH$`\alpha `$ 234 is cluttered with outflows and candidate sources. Twelve new HH objects are reported on here, with many of them lying along a preferential axis centred near LkH$`\alpha `$ 234 and orientated in a northeast-southwest direction. Their morphology suggests that at least some are part of a large-scale flow centred on the core of the NGC 7129 cluster.

Parsec - scale outflows from intermediate - mass YSOs show a number of similarities to those from their low - mass counterparts. In particular :
\- their lengths and degree of collimation appear comparable,
\- they share the same morphological trends such as decreasing frequency, increasing dimension and increasing complexity of HH emission with distance from the source.

The lengths of these large-scale outflows are usually comparable to the clump size of their associated clouds. As their expected lifetimes are much larger than the apparent dynamical timescales, this suggests that many have “blown-out” of the cloud complex.

Finally, it is evident that the transition from highly collimated jet - like flows to poorly collimated wide - angle outflows such as OMC 1 must lie at higher masses and luminosities than the sources studied here.

###### Acknowledgements.

FMcG and TPR acknowledge support from Enterprise Ireland. JBs research was supported by NSF grant AST-9819820, NASA grant ANG5-8108 (LTSA), and NASA grant NCC2-1052 (CU Center for Astrobiology).