text:
# Hydrodynamics of photoionized columns in the Eagle Nebula, M 16

## 1 Introduction

H$`\mathrm{II}`$ regions are common features of regions of massive star formation, such as Orion or the Eagle Nebula, M 16, in which ultraviolet photons from young massive stars photoionize and photodissociate the molecular gas clouds from which the stars originally formed. Emission from the ionized gas, in lines such as the Balmer series or \[O iii\]$`\lambda 5007`$Å, allows the regions to been seen at optical wavelengths. Images of these regions have many different forms, although they may often be characterised as a photoevaporating blister.

The basic stratification of gas in an H$`\mathrm{II}`$ region is well understood. Close to the central star or stars the gas is highly ionized, as a result of the strong photoionizing continuum. Atoms which recombine within this region are rapidly photoionized once again, and the gas is maintained at an equilibrium temperature close to $`10^4\mathrm{K}`$. At the edge of this region, where the photoionizing continuum spectrum has been essentially exhausted, is an ionization front (IF) which separates the ionized gas from atomic gas. The ionization front is generally narrow, compared to the overall scale of the region. Beyond the IF, the gas in the photodissociation region (PDR) is kept atomic by photons with wavelengths longer than 912 Å which do not have sufficient energy to ionize hydrogen, but are energetic enough to dissociate molecular hydrogen \[Bertoldi & Draine 1996\]. Eventually, the spectrum of photodissociating photons is also absorbed, and the gas can remain molecular.

However, this simple stratification belies the intricate structures often observed in H$`\mathrm{II}`$ regions. A feature of particular interest is the photoionized columns, or ‘Elephant Trunks’, seen at the edges of regions such as M 16, M 20 and NGC 3603. These intrusions of molecular gas into the high-pressure environment of the ionized nebula are potentially important sites for star formation, and indeed evidence has been mounting which associates young stellar sources with the tips of the columns in M 16 \[Hester et al. 1996, White et al. 1999\]. But the origin, structure and evolution of these columns has been the subject of controversy for many years. Do they result from pre-existing dense cores in the molecular gas or from instabilities in the IF/PDR structure? What is the timescale for their evolution?

### 1.1 The model of White et al.

White et al. \[White et al. 1999\] have studied the structure of the columns in M 16 using observations in a wide variety of wavebands. They suggest that the heads of the columns are pre-existing molecular clumps which are now in the process of dynamical collapse, driven by an overpressure in the ionized gas around the heads, which will soon result in the formation of stars. In this model, gas at the heads of the columns in M 16 started to collapse at most $`10^5\mathrm{yr}`$ ago, when the ionization front (IF) moving away from the exciting stars reached their surfaces. The velocity of the shock driven into the neutral material by the IF is approximately

$$v_\mathrm{s}\sqrt{\frac{P_{\mathrm{sn}}P_\mathrm{n}}{\rho _\mathrm{n}}}$$
(1)
where $`P_{\mathrm{sn}}/k_\mathrm{B}1.2\times 10^8\mathrm{cm}^3\mathrm{K}`$ and $`P_\mathrm{n}/k_\mathrm{B}3.5\times 10^7\mathrm{cm}^3\mathrm{K}`$ are the pressures inferred for the shocked neutral gas and for the gas at the head of the column, respectively, and $`\rho _\mathrm{n}2\times 10^5m_{\mathrm{H}_2}\mathrm{cm}^3`$ is the density of the gas in the head. White et al. calculate that $`v_\mathrm{s}1.5\mathrm{km}\mathrm{s}^1`$, and so that the shock will cross the head of the column in $`1.5\times 10^5\mathrm{yr}`$.

In consequence, White et al. argue that the columns have been subject to the incident radiation flux for a time no longer than this. For instance, the parameters given by White et al. imply that the flux of gas into the shock is $`5\times 10^{10}\mathrm{cm}^2\mathrm{s}^1`$ H nuclei, while the mass flux out through the IF is $`4\times 10^9\mathrm{cm}^2\mathrm{s}^1`$ H nuclei. Even allowing for the different areas of IF and shock, shocked gas should pile up between them to a total column of $`10^{23}\mathrm{cm}^2`$ in $`10^5\mathrm{yr}`$. If the clump had been subject to photoevaporation for more than $`10^5\mathrm{yr}`$, the shocked gas should now be observable in the CO or $`850\mu \mathrm{m}`$ data.

White et al. provide secondary support for the age they infer from detailed chemical and thermal modelling, which agrees with the observed properties of the molecular gas, and particular the low temperature of the densest molecular gas, and from the absence of diagnostics of protostars within the cold clumps. They argue that the distance between the ionization front and shock propagating into the clump will be far smaller than $`0.2\mathrm{pc}`$, so the dense material seen at the head of the column cannot itself be shocked gas.

Based on this model, they propose that the tips of the columns in M 16 may represent the best examples known to date for earliest stages of Class 0 protostellar development, which have not yet collapsed far enough to obey the full criteria for Class 0 \[André 1996\].

### 1.2 Other possibilities

Despite the arguments of White et al. \[White et al. 1999\], a lifetime of $`10^5\mathrm{yr}`$ still seems difficult to reconcile with the dynamical age of the other structures around the columns. Since the density of the ionized gas around the column is far smaller than that of the molecular gas, the main ionization front is likely to be D-type. As a result, the IF and its leading shock will be moving into unperturbed molecular material ahead at a velocity $`10\mathrm{km}\mathrm{s}^1`$. Indeed, since the age of the stellar cluster is $`12\times 10^6\mathrm{yr}`$, the average velocity of the ionization front over this lifetime can be no more than $`2\mathrm{km}\mathrm{s}^1`$. Even if the velocity of the IF were as high as $`10\mathrm{km}\mathrm{s}^1`$ as it propagated along the sides of the column, the IF must have passed the head of the column at least $`3\times 10^4\mathrm{yr}`$ ago, so an age substantially smaller than $`10^5\mathrm{yr}`$ seems unlikely. Similarly, Hester et al. \[Hester et al. 1996\] argue that the distribution of evaporating gaseous globules (EGGs) ahead of the column suggests that they are being photoevaporated on a timescale of a few $`\times 10^4\mathrm{yr}`$ as the column recedes, so the limit on their distribution relies on this photoevaporation timescale rather than on the lifetime of the finger itself. Finally, the presence of strikingly similar columns in other regions suggests that columns may survive for a significant fraction of the total age of the H$`\mathrm{II}`$ regions into which they intrude. None of these arguments excludes the picture developed by White et al.: however, the upper and lower limits on the lifetime are tight, and may come to exclude it in time.

In the present paper, we will explore longer-lived alternatives to the model of White et al. One is that the dense gas in the head is separated from the lower density gas in the column by a shock. The comparison between this case and the model of White et al. is illustrated schematically in Figure 1. White et al. rule out such a scenario on the basis of the narrowness of the ionization front/shock front layer. Indeed, the predicted flux of H nuclei through the shock, $`10^{10}\mathrm{cm}^2\mathrm{s}^1`$, means that, to have collected the observed $`30\mathrm{M}_{}`$ of gas in the head, it would have needed to propagate through approximately $`1\mathrm{pc}`$ of gas at the mean density of the neutral column over a timescale of $`4\times 10^5\mathrm{yr}`$: substantial fractions of both the size and lifetime of M 16. The collapse of the neutral material in the tangential direction may, however, decrease the time required to build up the mass \[Bertoldi 1989, cf\]. Structures of this form may result from density structures with initial contrast lower than those described by White et al., or from partial shadowing of the ionization front \[Williams 1999\]. The detailed structures found in numerical hydrodynamic simulations below are more complex than those shown in Fig. 1, but accord broadly with this scenario.

Another alternative is that the dense gas at the head of the column is instead confined by gravity, with the shock leading the ionization front having effectively bypassed the region and moved away downstream. White et al. find that self-gravity may be important at the head of the column. The size of the region dominated by the gravity of the observed molecular core is about $`0.13(M/30\mathrm{M}_{})(\sigma /1\mathrm{km}\mathrm{s}^1)^2\mathrm{pc}`$, where $`\sigma 1\mathrm{km}\mathrm{s}^1`$ is the effective sound speed of the gas (including contributions from turbulence and magnetic fields as well as the thermal sound speed) and $`M30\mathrm{M}_{}`$ is the mass of neutral gas in the head of the column.

In the following sections, we discuss numerical hydrodynamic simulations of these scenarios. These have been tailored to the central column in M 16, which has the most symmetrical appearance and for which the observational data have the highest spatial resolution. We follow Hester et al.’s nomenclature, and term it column II (it is referred to as $`\mathrm{\Pi }_2`$ by White et al.).

The structures which we find have some similarities to the equilibrium cometary cloud solutions discussed by Bertoldi & McKee \[Bertoldi & McKee 1990\], based on the collapse of isolated, near-spherical clumps. However, our time-dependent hydrodynamical models give useful information on flow stability and the presence of shocks. Similar structures have been found by García-Segura & Franco \[García-Segura & Franco 1996\], in their study of the development of H$`\mathrm{II}`$ regions in inhomogeneous environments, although as a result of the global scope of their simulations individual columns were poorly resolved. Mellema et al. \[Mellema et al. 1998\] model the formation of neutral tails behind isolated dense clumps in some detail. The models we present here have been developed for cases where the columns are close to the edge of an H$`\mathrm{II}`$ region, with parameters chosen to compare directly with those in M 16.

In Section 2, we discuss the derivation of suitable parameters for our numerical models from the observational data, and particularly the comparison of the pressures of ionized and neutral gas across the ionization front at the head of the column. In particular, we find that the observational evidence for an overpressure around the heads is weak. In Section 3, we discuss the numerical method we have used, and, in Section 4, the initial conditions for the simulations. In Section 5, we present the results of these simulations, and in Section 6 discuss their implications.

## 2 The Columns in the Eagle Nebula, M 16

### 2.1 Observational data

The columns in M 16 have been studied by many authors. Perhaps the best optical data is the HST narrow band imaging of Hester et al. \[Hester et al. 1996\]. We retrieved two H$`\alpha `$ images from this data set from the HST archive, and stacked and median-limited them to remove strong cosmic ray features and very fine-scale structures. The resulting image is shown in Figure 2. There are many striking irregularities in this image. Small clumps break the surface, and some appear fully separated from it. Close to the apex, the surface is puckered on a fine scale. Some of the clumps harbour K-band sources, and are identified as protostars by Hester et al. Away from the head, the surface irregularities seem to grow in width, with bright regions where their surfaces face towards the source of ionization. In addition, we note the low-magnitude striations which are seen along lines approximately perpendicular to the nearest surface.

White et al. \[White et al. 1999\] combined this optical data with observations of M 16 in molecular lines, millimetre and sub-millimetre continuum and the mid-infrared, to present a remarkably complete picture of properties of the gas and dust in and around the columns. Within the columns, molecular gas is probed by line emission and the associated dust is seen in millimetre/sub-mm emission, while H$`\alpha `$ and radio continuum observations are sensitive to structures in the surrounding ionized gas.

Levenson et al. \[Levenson et al.2000\] have observed excited H2 emission in M 16, to study the properties of the PDR which leads the ionization front. They find that their observations are consistent with a stationary photodissociation region, around $`2\times 10^{17}\mathrm{cm}`$ deep (although since the properties of the incident radiation field are not well determined, it is possible that shocks driven by the IF might also contribute to the collisional component of this emission).

### 2.2 Comparison with photoevaporated winds

In Figure 3, we show cuts through this image. The first is along the axis of the large-scale column which is vertical in Fig. 2. The values are averages for a strip between 345 and 355 on the $`x`$-axis. The second is perpendicular to this, averaged for a strip between 395 and 405 on the $`y`$-axis.

In Fig. 4, we show the radial emission measure profile derived from isothermal wind models, as described by Henney & Arthur \[Henney & Arthur 1998\]. These authors found that such models gave a very good fit to H$`\alpha `$ emission profiles of gas photoevaporated from proplyds in the Orion nebula. We have here been rather less strenuous in the correction of the observational data for contamination by continuum and \[N ii\] emission, but our aim in this section is simply to get an indication of the properties which pertain at the head of column II.

Comparing the model profiles with Fig. 3, we see that their overall forms are very similar. The central brightness of the columns appears rather lower, relative to that at large distances, than the model would predict, but that may be explained by the decreasing intensity of illumination as the surface of the column points away from the head and by the occulting of background emission by the column. The sharpness of the inner edge of the profile suggests that the column is slightly closer to us than is the photoionizing star.

More fundamentally, if (following Hester et al.) we take the radius of curvature of the head of the column to be 90 pixels, then the breadth of the observed H$`\alpha `$ profile is difficult to account for. Not even the most extreme of the models in Figure 4 explains the width of the transverse light distribution. For comparison, Sankrit & Hester \[Sankrit & Hester 2000\] find that power-law models with a radius of curvature of $`2\times 10^{17}\mathrm{cm}150`$ pixels fit the emission profile in several emission lines, which confirms that the effective radius of curvature is larger than that apparent in the images. Megeath & Wilson \[Megeath & Wilson 1997\] also find that the radio emission measure is more extended in their observations of photoevaporating clumps in NGC 281 West than for a $`n_e1/r^2`$ density model.

One possible explanation for the M 16 results is that the head of the observed column is more elongated along the line of sight than perpendicular to it. Alternatively, if the column is directed into the plane of the sky \[Pound 1998\], the inner peak of the emission measure profile may be hidden from view, although the inclination does not appear to be that extreme.

### 2.3 Inferred density and radiation field

For the model ionized flow profiles, isothermal flow with initial Mach numbers of $`1`$, $`1.2`$, $`2`$ and constant velocity flow, the effective integration distance (i.e. the peak emission measure divided by the square of the peak density) is 0.8, 1.0, 1.3, and $`\pi /2`$ times the radius of curvature, respectively. This may be compared to the value of one-half assumed by Hester et al. (since the profiles are well resolved, the correction due to resolution effects is minimal). When the uncertainty in the local radius of curvature is included, it seems very likely that Hester et al.’s value of $`4000\mathrm{cm}^3`$ for the electron density at the head of column II may have been overestimated. Indeed, Hester et al.’s own photoionization modelling suggested that the density at the head of the column is $`2000\mathrm{cm}^3`$ \[Sankrit & Hester 2000, although this has subsequently been revised by\]. As a result, the ionized gas may not in fact have a greater pressure than the effective pressure in the neutral gas at the head of the column.

The density at the head of the column can also be used to determine the incident ionizing flux. The thickness parallel to the radiation field is roughly $`\mathrm{}0.04\mathrm{pc}`$ from Figure 3. Balancing the number of recombinations in the flow with the intensity of the incident ionizing radiation, $`J`$, suggests that $`J\alpha _\mathrm{B}n_\mathrm{e}^2l1.3\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$, for a peak electron density of $`n_\mathrm{e}2000\mathrm{cm}^3`$. It is difficult to follow White et al.’s analysis of their VLA data, but assuming that the integrated $`48\mathrm{mJy}`$ emission corresponds to that occulted by the head of the column, which has a radius of curvature of $`0.04\mathrm{pc}`$ \[Hester et al. 1996\], the incident ionizing intensity predicted by the formulae of Lefloch, Lazareff & Castets \[Lefloch, Lazareff & Castets 1997\] is $`J1.5\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$. These values are mutually consistent, but smaller than those found by White et al., and also smaller than the values, inferred from stellar censuses, of $`4.4\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$ \[Hester et al. 1996\] or $`2.6\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$ \[Pound 1998, from the total flux given by\]. However, this difference can easily be accounted for by uncertainties in the ionizing flux of the stars, as the result of dust absorption or photoionization between the ionizing stars and the columns, or by the geometry of the region.

In our models, we neglect the role which a PDR may have in the dynamics of the region. The depth of the PDR in M 16, as derived by Levenson et al. \[Levenson et al.2000\], is close to the radius of curvature of the finger tips, so this assumption is only marginally valid. However, since the column density inferred is close to the stationary models \[Bertoldi & Draine 1996, rather than being substantially affected by advection, cf.\], it seems likely that the jump conditions across the combined IF/PDR structure are determined principally by the flux incident on the IF.111If, in other cases, the depth of a PDR in pressure equilibrium with the photoionized flow were larger than the radius of a globule, the PDR would be important in development of the globule; if it were larger than the inter-globule separation, the development of photoionized columns might be suppressed, with the internal structure of the molecular cloud only observable in atomic and molecular species.

### 2.4 Derived parameters

Based on the observations, we will adopt the following parameters for our modelling of column II. The pressure in the ionized gas at the head of the clump is $`P_\mathrm{i}/k_\mathrm{B}3\times 10^7\mathrm{cm}^3\mathrm{K}`$; this ionized gas is assumed to have an isothermal sound speed of $`10\mathrm{km}\mathrm{s}^1`$. The density in the molecular gas in the core at the head is $`n(\mathrm{H}_2)=2\times 10^5\mathrm{cm}^3`$ while that in the column of gas behind is typically 10 times smaller than this. White et al. suggest that while the temperature of the gas is $`20\mathrm{K}`$ throughout the column, the effective pressure within the columns is dominated by the effects of small-scale magnetic fields and turbulence, as evidenced by the smooth profiles of molecular lines with widths $`2\mathrm{km}\mathrm{s}^1`$. We cater for this in our dynamical models by assuming that the neutral gas behaves as if it were isothermal with a temperature of $`200\mathrm{K}`$ (note that the turbulent contribution to the pressure in the molecular gas also reduces the importance of heating in the PDR for dynamical modelling). The total mass in the core of column II is $`31\mathrm{M}_{}`$ and its radius is $`0.085\mathrm{pc}`$. We assume a photoionizing flux of $`10^{11}\mathrm{cm}^2\mathrm{s}^1`$, impinging on the column parallel to its axis.

## 3 Numerical Method and Assumptions

We use the hydrodynamical code described by Williams \[Williams 1999\], which is based on the second-order Godunov algorithm by Falle \[Falle 1991\]. Ionization and recombination are treated by including the equations

$`n{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}`$ $`=`$ $`an(1x)J\alpha _\mathrm{B}n^2x^2`$ (2)
$`{\displaystyle \frac{\mathrm{d}J}{\mathrm{d}z}}`$ $`=`$ $`an(1x)J.`$ (3)
Here $`n`$ is the density of hydrogen nucleons, $`x`$ the ionization fraction and $`J`$ is the number flux of ionizing photons, assumed to propagate parallel to the $`z`$ axis. The photoionization cross-section of hydrogen, $`a`$, and the case B recombination coefficient, $`\alpha _\mathrm{B}`$, are taken as constants. The electron density, $`nx`$, is included as an advected variable in the hydrodynamic solver, and the ionization and recombination equations included by operator-splitting, implemented at the end of each step (and also for the intervening half-steps). The equations are integrated using an implicit scheme to advance $`x`$ in time and $`J`$ across each grid cell.

The code does not include a treatment of the diffuse ionizing radiation field (beyond using the case-B recombination coefficient). In the on-the-spot approximation, the diffuse field is found to be $`15`$ per cent of the local direct radiation field \[Cantó et al. 1998\]. The characteristic length for reabsorption of diffuse photons is $`\mathrm{}J/\alpha _\mathrm{B}n^20.05\mathrm{pc}`$ (dependent, or course, on spatial variations of $`J`$ and $`n`$), so the diffuse field might be expected to have some effect even on structures as large as the columns. More detailed treatments predict that the diffuse radiation is a smaller fraction of the direct field away from the edge of the global photoionized region \[Rubin 1968\]. Close to the edge, where relative size of the diffuse field is larger, the soft spectrum of the diffuse field will result in a smaller path length to absorption than for the direct field. The low contrast of the surfaces of the observed columns which do not point towards the ionizing source is empirical evidence that the diffuse radiation field would is adequately treated by the on-the-spot approximation in the present case.

Since the internal structures of IFs are not resolved in the simulations, most of the gas will be either almost completely ionized or fully neutral. Intervening regions with intermediate ionization are likely to result principally from numerical mixing, so physically detailed schemes which assume the gas within the computational cells is reasonably uniform are not appropriate. Instead, we set isothermal sound speed in the gas to

$$c_\mathrm{s}=(1+99x)^{1/2}\mathrm{km}\mathrm{s}^1.$$
(4)
This expression results from assuming that in cells with $`0<x<1`$, the gas is split into separate regions of fully ionized and fully neutral gas, in pressure equilibrium with each other, and so gives a better account of the sub-grid structure at the (dynamically significant) D-type IF.

## 4 Initial and boundary conditions

The main model presented here is calculated in cylindrical symmetry, with coordinates $`(r,z)`$. A uniform $`128\times 448`$ grid covers physical dimensions $`0<r<0.5\mathrm{pc}`$ and $`0<z<1.75\mathrm{pc}`$. The boundary conditions are reflective on the axis, $`r=0`$, as required by symmetry. On the boundary towards the source, $`z=1.75\mathrm{pc}`$, free-flow boundary conditions are used when the flow is outwards, mirror conditions when it is inwards, although as the outward flow is generally supersonic this has little influence on the flow inside the grid. Along the outer radial boundary, $`r=0.5\mathrm{pc}`$, the conditions beyond the grid are either free-outflow/zero-inflow if the ionization fraction is larger than 0.1, or mirror conditions if it is less than this value. This somewhat artificial criterion changes the behaviour at the boundary at the IF, in order to model an ionized flow which is divergent on large scales. The effect depends little on the specific limiting value, as away from the IF the gas is in general either almost entirely ionized or almost entirely neutral.

On the edge of the grid towards the molecular cloud, $`z=0`$, the boundary condition is set by assuming that the material retains its initial state, i.e. in most of our examples a flow of gas towards the IF. As the solution develops, a shock often becomes captured at the boundary, which reduces the inflow speed and increases the inflow density. The flux of neutral gas into the grid through the plane $`z=0`$ is, in all cases, small enough that it could not prevent the IF from entering if the entire grid were full of neutral material. In many cases, the density of the ionized gas becomes large enough to absorb all the ionizing photons within the size of the grid. Where the ‘molecular’ boundary covers the entire upstream edge of the grid, a plane-parallel IF is an equilibrium solution. If the ionized flow is allowed to escape from the side of the region, it is not obvious what the form of the equilibrium will be, if indeed there is one. One possibility is that the IF surface will remain smooth, and a steady solution will be found in which the variation in pressure at the IF is balanced by changes in flow velocity via the Bernoulli effect.

However, these IF are believed to be unstable to corrugation modes on both small and large scales, although for wavelengths larger than $`6\times 10^{15}n_3^1\mathrm{cm}`$ recombination decreases the growth rate of these instabilities decreases rapidly \[Kahn 1958, Axford 1964, Sysoev 1997\]. This wavelength is comparable to our grid scale of $`10^{16}\mathrm{cm}`$ for typical ionized gas densities at the IF of $`10^3\mathrm{cm}^3`$. We now present numerical models which illustrate what may happen for a range of initial conditions.

The models we have calculated have initial conditions split into two regions: ionized gas with a density $`n(\mathrm{H}\mathrm{total})=200\mathrm{cm}^3`$ towards the source of ionizing photons and moving towards it at $`10\mathrm{km}\mathrm{s}^1`$ relative to the neutral gas behind. This neutral gas has density, $`n(\mathrm{H}\mathrm{total})=2\times 10^4\mathrm{cm}^3`$ (n.b., not $`n(\mathrm{H}_2)`$, see below). Except where noted, in the initial conditions the neutral gas is assumed to move at $`2\mathrm{km}\mathrm{s}^1`$ towards the ionizing source in order to keep the ionization front broadly at rest within the grid, with the ionized gas assumed to stream towards the radiation source $`10\mathrm{km}\mathrm{s}^1`$ more rapidly than this. Pressure equilibrium holds where ionized and neutral gas are alongside each other. We take the ionizing flux incident at the top of the grid in the figures to be $`1.5\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$, since recombinations in the upstream flow are found to reduce the ionizing radiation incident on the head of the column somewhat. The Strömgren distance in the ionized gas, $`3\mathrm{pc}`$, is sufficiently great that a column the length of the computational grid can be maintained with almost complete ionization for the initial conditions.

The boundary between these regions is assumed either to be plane (to study the development for small perturbations) or to include a column of dense gas (to study the likely long-timescale development of the columns observed in many regions).

## 5 Results

We will now present the results of four different scenarios for the development of photoevaporated columns.

is based on the model of White et al. \[White et al. 1999\], including initial density stratification and gravity,

has an initial short cylindrical column,

has a 10 per cent deficit in photoionizing flux close to the axis,

has dense neutral gas initially filling the grid, with gravitational forces from a $`30\mathrm{M}_{}`$ object included.

#### Case I:

In Figure 5, we present the results of a first simulation in which the flow conditions modelled on those inferred by White et al. The initial conditions have a dense core of material at the head, and a lower density column behind. A gravitational potential from a $`30\mathrm{M}_{}`$ mass, smoothed at a radius of $`0.08\mathrm{pc}`$, is centred on the core of the head. The uniform gas at the head of the column initially falls in, as it is not in equilibrium with the gravitational field, but does not collapse far because of the smoothed core of the field.

By the time of the second frame, $`10^5\mathrm{yr}`$ after the start of the simulation, the shock at the head of the the column has passed through much of the dense gas. The shock proceeds down the column in later frames. Soon the ram pressure of the radiation field is sufficient to push the dense gas away from the gravitating core. The amount of gas with density above $`10^5\mathrm{cm}^3`$ in these models remains roughly constant, and so the gravitational field would in reality respond to this movement. Nevertheless, it is striking that the density field in the column remains close to that inferred by White et al. even when the head is a considerable distance away from the gravitating centre, rather than requiring the head of the column to be caught within $`10^5\mathrm{yr}`$ after the IF reaches it.

In the remaining examples, we investigate possibilities by which such structures may have formed, and what conditions are necessary for them to be long-lived.

#### Case II:

In Figure 6, we illustrate a case in which a short column of gas is initially present. The radius of the column is initially $`0.125\mathrm{pc}`$ and it stands proud by $`0.5\mathrm{pc}`$ ahead of the main body of molecular material. The column is again seen to be long lived, but highly variable. Certain broad features characterise the structure. The highest densities in the neutral gas are seen close to the head of the column, and reach values similar to those inferred for the column tip by White et al.

The surface of the column is convoluted. Neutral clumps split off from from the column at different stages (e.g. close to the base of the column in Figure 6f) and linear structures are seen where different streams of ionized gas interact. These are reminiscent of the observed clumps and striations in the flow (see Figure 2 and Section 2). The striations are sometimes seen to be triggered by surface inhomogeneities resulting from the discrete nature of the grid. However, the growth of these instabilities appears to be a generic feature of the development of obliquely-illuminated ionization fronts \[Vandervoort 1962\], even when account is taken of the effects of recombination in the upstream flow (Williams, in preparation).

In the frames shown in Figure 6, the dense column becomes more marked with time. Initially, the ionizing field generates a strong flow away from the head of the column, a low density region ahead of it and a net acceleration away from the ionizing source \[Oort & Spitzer 1954, the well-known ‘rocket effect’,\]. However, we see here that short lengthscale which characterises the photoevaporative flow from the head means that its pressure decreases rapidly, and eventually a termination shock forms. The convergence of this termination shock at the axis leads to the formation of an over-dense region at larger distances from the column. Such over-dense regions lead to the formation of some thin filaments connected to the main ionization front \[e.g., those seen close to the base of the column in Figure 6(b)\]. Filaments at the head of the column can allow the column to move forward once again, although the rigidity provided for the filaments by the reflecting boundary condition on-axis may have a significant influence by preventing non-axisymmetric instability modes.

#### Case III:

In Figure 7, we show the development of columns for a case where there is a 10 per cent weakening of the impinging ionizing radiation field for a small radius around the axis. A narrow feature in the radiation field such as this might formed by an upstream overdensity: a photoevaporating globule, for instance, or dense structures close to the stars caused by interacting stellar winds.

A strong column soon forms in this case. Once it has formed, its evolution is not greatly influenced by the small perturbation in the radiation field. The structure is broadly similar to those seen in the previous two cases, although the column has a somewhat higher aspect ratio and the shock remains within the computational grid for a longer period.

#### Case IV:

In Figure 8, we take uniform dense gas to nearly fill the grid, and include a gravitational field equivalent to a $`30`$-$`\mathrm{M}_{}`$ core smoothed at a radius of $`0.02\mathrm{pc}`$. In this case only, the neutral gas is initially at rest. As the simulation develops, mass first falls into the gravitational potential, to form a dense core close to hydrostatic equilibrium. A shock and ionization front driven by the impinging ionizing radiation propagate over this core, generating a column structure. This column is very long-lived, tending to form a cometary structure at late times. The length of this cometary structure may in reality be limited by the non-plane parallel ionizing radiation field resulting from diffuse emission and the multiple ionizing stars \[Pavlakis et al. 2001\].

This model seems similar to that discussed by White et al. However, it will be noted that by the time the ionization front reaches the dense clump, it interacts with the reflection of the main shock front from the dense core, rather than driving a separate shock itself. Shocks are always present within the solution, but are a characteristic of the dynamically disturbed nature of the flow.

Here we have imposed a separate potential, rather than calculating the gravity self-consistently, primarily as a computational convenience. However, it should be noted that in the Eagle columns, a reasonable fraction of the gravitational potential may be generated by protostellar clumps which have effectively decoupled from the flow. The small smoothing radius means that the gravitational field could not be generated by the self-gravity of a stable isothermal molecular gas core. It is possible that protostellar cores may be formed from gravitational instabilities in the swept up gas shell (Francis & Whitworth, in preparation), or be present in the initial molecular cloud.

Note, however, that if the heads of the columns were formed from individual isothermal cores confined by self-gravity, then these cores must have been surrounded by significantly hotter, less dense gas if the overall structure was not to be unstable to gravitational collapse. In this case, the source of the $`10^4\mathrm{cm}^3`$ gas in the column behind the column is unclear.

### 5.1 Discussion

In each of the cases presented in the present section, structures strikingly reminiscent of those observed by White et al. \[White et al. 1999\] are maintained within the computational grid for periods of several $`10^5\mathrm{yr}`$. The overall structures observed for models run for a range of similar parameters were broadly the same, suggestion that our results are not dependent on the details of our numerical or physical assumptions. This suggests that structures seen in M 16 are in a near-steady state, with the head of the column confined by the rocket effect at the front and by the ram pressure of the gas being swept up, by self-gravity, or a combination of the two. The structures can survive for several times $`10^5\mathrm{yr}`$, longer than the shock crossing time of the head of the column, and similar to the dynamical age of M 16 as a whole. However, the influence of the symmetry axis on the stability of these columns means that this interpretation is tentative, until fully three-dimensional models are available.

In Figure 9, we show the H$`\alpha `$ emissivity expected from the gas within Case I after $`5\times 10^5\mathrm{yr}`$ from the start of the simulation, calculated as in Steffen et al. \[Steffen et al. 1997\]. This image shows the presence of bright rims and strong linear features in the flow away from the surface.

Since the structures appear to be close to equilibrium, the observations contain limited information about the mechanisms by which they formed: as we have shown, a variety of initial conditions can lead to very similar final structures. However, the observations seem not to require highly stratified density structures to have been present in the original molecular cloud. Indeed, it seems possible that the columns result from instabilities of the rim of a champagne-type expanding H$`\mathrm{II}`$ region \[Tenorio-Tagle 1979, Bodenheimer et al. 1979\]. It is difficult to make substantial analytic progress in the study of such surface perturbations, as a result of the nonlocal nature of the system. As we have seen in the simulations, density enhancements resulting from fluctuations in the surface a significant time ago can result in long term variations in the ionizing flux incident on the surface a considerable distance from the original perturbation. While it is certainly true that the assumption of cylindrical symmetry enhances the variation of ionizing flux on the axis, it seems likely that the long lags between perturbation and response inherent in the system will in general lead to overstabilities of the flows, when studied on large scales.

The interpretation we have developed, however, has its own coincidental features. Most particularly, the mass of shocked molecular gas at the head of the columns is curiously close to the Jeans mass. There is no obvious dynamical reason for the mass of dense gas to settle on any particular value: indeed, we have seen the mass of dense gas to be highly variable in our simulations. A possible explanation for why the mass is close to the limit of collapse is that, in fact, the mass at the head increases as the column erodes, but once it reaches the limit of stability it collapses. The dense, self-gravitating globules produced will then no longer be subject to significant radiation or hydrodynamic forces, and can be left behind by the further photoevaporation of the column. In this manner, an ionization front can lead to stars forming in an episodic manner over a substantial period, compared to the one-off radiation-induced collapse of clumps which has been studied previously \[Oort 1954, Sandford et al. 1982, Bedijn & Tenorio-Tagle 1984, Bertoldi 1989, Bok, as referred to by\]. Note that the ionizing radiation incident on a clump will increase gradually when the main ionization front is of D-type, rather than suddenly as often assumed. The shock which precedes the IF may however result in a thin gas layer unstable to gravitational collapse \[Elmegreen & Lada 1977, Whitworth et al. 1994\] and may trigger the collapse of any pre-existing self-gravitating clumps \[Megeath & Wilson 1997\]. To develop these ideas further requires models in which the self-gravity of the molecular gas is included.

## 6 Conclusions

In this paper, we have presented hydrodynamical models of the development of photoionized columns, which agree well with the observations of the columns in M 16. Our models contrast with the results of White et al., who interpret observations of the columns in M 16 as requiring that they have been exposed to the radiation field for less than $`10^5`$ years, the dense material in the head of the column being a gravitationally-bound condensation from the initial molecular cloud.

It remains to determine unambiguously how columns such as those in M 16 form, and how they evolve, if indeed there is any single mechanism. It seems likely, however, that these structures are often long lived components of the environment of massive stars, and of great importance in the formation of lower mass stars in these environments, and so they warrant further detailed study.

## Acknowledgments

The work of RJRW is supported by the PPARC through the award of an Advanced Fellowship. We would like to thank the referee, Matt Redman, for a helpful and constructive report.


format: text

text:
# Bond-Stretching-Phonon Anomalies in Stripe-Ordered La1.69Sr0.31NiO4

## Abstract

We report a neutron scattering study of bond-stretching phonons in La1.69Sr0.31NiO4, a doped antiferromagnet in which the added holes order in diagonal stripes at 45 to the Ni-O bonds. For the highest-energy longitudinal optical mode along the bonds, a softening of 20% is observed between the Brillouin zone center and zone boundary. At 45 to the bonds, a splitting of the same magnitude is found across much of the zone. Surprisingly, the charge-ordering wave vector plays no apparent role in the anomalous dispersions. The implications for related anomalies in the cuprates are discussed.

There is resurgent interest in the role of phonons with respect to the high-temperature superconductivity found in layered copper-oxides Zhao et al. (1997). Particularly striking are the anomalies in high-energy optical modes observed by neutron scattering in La2-xSrxCuO4 Pintschovius and Reichardt (1998); Pintschovius and Braden (1999); McQueeney et al. (1999a) and YBa2Cu3O6+x Pintschovius and Reichardt (1998); Reichardt (1996); Petrov et al. ; McQueeney et al. . There have been various speculations as to whether the observed phonon anomalies might be related to instantaneous charge inhomogeneities, particularly those in the form of stripes Emery et al. (1999). Recently, detailed analyses of electron-phonon interactions in a dimerized stripe phase have been reported Castro Neto (2001); Park and Sachdev (2001).

One way to learn about the effect of charge stripes on lattice dynamics is to study a model system with well-defined stripe order. Here we present the first single-crystal study, to our knowledge, of the bond-stretching phonon modes in such a system, specifically La2-xSrxNiO4 with $`x=0.31`$. The stripe order in Sr-doped nickelates has been characterized in detail by neutron diffraction, and the most recent summary of results is given in Yoshizawa et al. (2000). For $`x0.22`$, the crystal structure is tetragonal, consisting of a body-centered stacking of NiO2 planes. Within the NiO2 planes, the charge stripes run diagonally along either  or $`[1\overline{1}0]`$ directions, at 45 to the Ni-O nearest-neighbor bonds, which extend along  and  directions. (In contrast, the charge stripes observed in superconducting La1.6-xNd0.4SrxCuO4 run parallel to the Cu-O bonds Ichikawa et al. (2000); however, for $`x0.06`$, the stripes inferred to exist in La2-xSrxCuO4 have the diagonal orientation of the nickelates Wakimoto et al. (1999); Fujita et al. .) The maximum transition temperatures for charge-stripe and magnetic order occur at $`x=\frac{1}{3}`$ Yoshizawa et al. (2000); Cheong et al. (1994). Well below the charge-ordering temperature, $`T_{\mathrm{co}}`$, the nickelates have very large resistivities Katsufuji et al. (1996, 1999), consistent with all of the added holes being localized in charge stripes.

We focus on the highest-energy longitudinal optical modes propagating along the  and  directions. These modes show only weak dispersion in stoichiometric La2NiO4 Pintschovius et al. . In contrast, we observe a softening of 20% along  on moving from the Brillouin zone center to the zone boundary, quite similar to that found in La2-xSrxCuO4 Pintschovius and Braden (1999); McQueeney et al. (1999a); along , a splitting of the same magnitude is observed over much of the zone and, in particular, at the zone boundary. These results are important for two reasons. 1) The observed anomalies must be associated with the local charge inhomogeneity. They are induced by the hole doping, and the holes are localized in the stripes. 2) There is no evidence that the charge-ordering wave vector, $`𝐪_{\mathrm{co}}`$, plays a special role. If the phonon anomalies were related to collective phase fluctuations of the charge stripes, as in a conventional charge-density-wave system Schulz (1978), then one might expect them to appear at $`𝐪_{\mathrm{co}}`$. With the absence of a collective signature, it seems likely that the dominant effects involve local interactions between charge and lattice fluctuations.

Our La2-xSrxNiO4 crystal, grown by the floating-zone method, is cylindrical, with a diameter of 6 mm and length of 30 mm; the Sr concentration of $`x=0.31`$ was confirmed by inductively-coupled plasma analysis (with an uncertainty of $`\pm 0.01`$). At room temperature, the lattice parameters of the tetragonal unit cell are $`a=3.84`$ Å and $`c=12.70`$ Å. The charge and spin ordering transitions were confirmed to be consistent with previous workYoshizawa et al. (2000) (magnetic and charge-order transitions of approximately 160 K and 235 K, respectively) by neutron diffraction measurements performed at the JRR-3M reactor in Tokai, Japan. The inelastic-neutron-scattering measurements of phonons were performed on the 1T triple-axis spectrometer at the Orpheé reactor of the Laboratoire Léon Brillouin in Saclay, France. The spectrometer is equipped with a Cu (111) monochromator and a pyrolytic graphite (PG) (002) analyzer, each of which is both vertically and horizontally focused. The analyzer was set to detect neutrons with a final frequency of 3.55 THz (14.7 meV). A PG filter was placed after the sample to minimize unwanted neutrons at harmonic wavelengths.

Figure 1(a) shows scans of scattered intensity versus excitation frequency for phonons propagating in the  direction, parallel to the in-plane Ni-O bonds. The scans are taken at wave vectors along the line $`𝐪=(h,0,0)`$ in the first Brillouin zone, moving from the zone center ($`h=0`$, bottom) to the zone boundary ($`h=1`$, top), where the components of the wave vector are measured in reciprocal lattice units, $`(2\pi /a,2\pi /a,2\pi /c)`$. (The effective zone boundary for a single NiO2 plane is at $`h=0.5`$; the fact that the actual boundary is at $`h=1`$ results from the body-centered stacking of the layers.) Each scan is dominated by a single, well-defined peak, that disperses from approximately 21 THz at zone center to less than 18 THz half-way across the zone, and then back up again. The lines through the data points are fitted gaussian peaks on top of a background that is taken to be independent of wave vector and frequency.

Figure 1(b) shows scans measured along the  direction. Because of twinning of the stripe domains, we simultaneously probe phonons propagating parallel and perpendicular to the stripes. Again, we focus on the mode that starts at about 21 THz at zone center; the mode at 14 THz involves bond-bending motion, and does not exhibit any notable doping-dependent behavior. \[The intensity variation seen in the figure for this mode is associated with its sensitivity to the $`L`$-component of the wave vector Q, which varies from point to point as noted in Fig. 1(d).\] On moving from zone center to zone boundary, we note that the highest-frequency mode varies little in frequency, but it loses intensity. Concommitant with this, signal appears in the 17–18 THz regime, growing in strength as the zone boundary is approached.

Using the gaussian fits as a smoothed version of the data, we present our results as intensity plots in Fig. 2(a) and (b). Here, the intensity has been multiplied by $`\nu /Q_{}^2`$ ($`Q_{}`$ is the component of Q along the phonon propagation direction) to correct for sensitivities of the neutron-scattering cross section. Dramatic differences are observed in the dispersions of the stripe-ordered system compared to those in stoichiometric La2NiO4, which are indicated by the gray circles Pintschovius et al. . Contrary to naive expectations based on the Peierls-distortion model Schulz (1978), there is no particular anomaly at the charge-ordering wave vector $`𝐪_{\mathrm{CO}}=(0.31,0.31,0)`$ \[see Fig. 2(b)\]. Instead, we find a splitting of the mode along  over a substantial part of the zone. The size of the dispersion in the  direction is about the same as the  splitting, and is essentially identical to the magnitude of the dopant-induced softening reported for superconducting La1.85Sr0.15CuO4 Pintschovius and Braden (1999); McQueeney et al. (1999a). The softened modes are consistent with a previous study of the phonon density-of-states in La2-xSrxNiO4 McQueeney et al. (1999b), where a dopant-induced peak was observed at $`75`$ meV ($`18`$ THz).

Integrating the corrected intensity data of Fig. 2(a) and (b) over frequency yields the open circles shown in Fig. 2(c) and (d), respectively. The experimental results are compared with intensities calculated from the interatomic-potential model used to describe the phonons in La1.9NiO3.93 Pintschovius et al. . (An overall normalization of the calculated to the measured intensities has been applied in each panel.) The general consistency between the measured and calculated integrated intensities suggests that the distribution of weight for the anomalously-softened bond-stretching modes does not extend significantly below 16 THz.

To interpret our results, let us consider some simple phenomenological models. The high-energy bond-stretching phonons involve the motion of the in-plane oxygens between the much heavier nickel atoms. Suppose we assume the nickel atoms to be infinitely heavy and consider only nearest-neighbor forces; then each oxygen acts like an Einstein oscillator, with no dispersion of its vibrational frequency. Correcting for the finite mass of the nickel would result in a slight decrease in frequency towards the zone boundary; however, accounting for Coulomb repulsion between neighboring oxygen ions would counter that with an increase in frequency at the zone boundary, where neighboring oxygens move opposite to one another. Ignoring these corrections, a simple Einstein model gives a rough approximation of the measured dispersions in stoichiometric La2NiO4 Pintschovius et al.  \[see Fig. 2(a) and (b)\].

To go further, we will restrict ourselves to a linear Ni-O chain model. The softening of phonons at the zone boundary can be described by introducing a force between nearest-neighbor O ions with a negative force constant, $`k_{\mathrm{br}}`$ Pintschovius et al. (1989, ) \[see Fig. 3(a)\]. This force phenomenologically incorporates a particular type of electron-phonon coupling. An example of the dispersion from such a model is shown in Fig. 3(c). The eigenvector for the softened zone-boundary mode, sketched in Fig. 3(e), involves linear breathing motion of the O about the Ni.

Next, suppose that $`k_{\mathrm{br}}`$ acts only between one pair of oxygens out of three \[see Fig. 3(b)\]. The result is a splitting of the modes, as illustrated in Fig. 3(d). Note that the intensity of the softened mode is strong at the zone boundary, but goes to zero at zone center. This effect is achieved without a dimerization of the lattice. As shown in Fig. 3(f), the eigenvector of the softened mode at zone boundary involves an isolated, local breathing motion of a pair of O about a Ni site. The behavior of this mode does not require coherence between motions in different unit cells.

It seems likely that both of these models have some degree of relevance to our observations; however, it is difficult to associate them in a consistent or unique way with specific details of the measured dispersions. Of greater interest is the nature of the electron-phonon coupling modelled by $`k_{\mathrm{br}}`$. In model calculations for La2-xSrxCuO4, Falter and Hoffmann Falter and Hoffmann (2001) have emphasized the importance of ionic charge fluctuations in response to atomic displacements and the irrelevance of Fermi-surface nesting effects of the type once proposed by Weber Weber (1987). Given the experimental evidence for poor electronic screening of phonons in optimally-doped cuprates Homes et al. (2000), extended Coulomb interactions should also be important.

In the present case of the nickelate, little metallic screening is expected due to the gap of $`>0.1`$ eV in the optical conductivity at low temperature Katsufuji et al. (1996); Pashkevich et al. (2000). The bond-stretching modes involve polar fluctuations of the negative oxygen ions against the positive nickel ions, and so should couple to charge fluctuations. The lowest-energy channel for charge fluctuations must involve the dopant-induced holes, and since the holes are segregated into stripes, charge fluctuations must be associated with some form of stripe fluctuations.

In a conventional charge-density-wave (CDW) system, one might expect the mode along  to couple to fluctuations of the phase of the CDW with respect to the lattice Schulz (1978). In the present case, the absence of a characteristic wave vector suggests that local charge fluctuations may be more relevant for the bond-stretching modes. Local fluctuations were predicted by Yi et al. Yi et al. (1998) in a random-phase-approximation treatment of Hartree-Fock stripes. Nevertheless, our results do not rule out the possibility of weak, low-energy, collective phase fluctuations of the stripes.

Do the effects we observe depend on static ordering of the stripes? Figure 4 shows a comparison of phonons at zone-center and zone boundary along  measured at 10 K and 250 K, which is above the charge-ordering temperature. One can see that at 250 K, the zone-center mode has softened slightly but remains strong and well-defined; however, at the zone boundary, the intensity remains spread out in frequency. Thus, the anomalies do not disappear with the static order.

Rather similar phonon softenings were observed in the sample identified as La1.9NiO3.93 Pintschovius et al. ; Pintschovius et al. (1989). Despite the specified stoichiometry, that sample exhibited superlattice peaks of the type $`(H/2,K/2,L/2)`$ with $`H`$, $`K`$, and $`L`$ odd, consistent with the stage-2 ordering of interstitials found in La2NiO4.105 Tranquada et al. (1994). The latter compound exhibits evidence of incipient stripe correlations, but no stripe order Tranquada et al. (1997). This is further evidence that local interactions are sufficient to yield the phonon anomalies.

To summarize, we have observed doping-induced anomalies in bond-stretching modes of La1.69Sr0.31NiO4, a compound in which the doped holes order in stripes at low temperature. The anomalies persist over a substantial portion of the Brillouin zone, with no obvious signature of collective stripe fluctuations. The splitting of the mode along  is similar to recent observations in superconducting Ba0.6K0.4BiO3 Braden et al. , where charge inhomogeneity is also expected to be relevant. These results provide support for a connection between the phonon anomalies observed in the cuprates and instantaneous charge inhomogeneity.

We gratefully acknowledge stimulating discussions with R. Werner and helpful comments from P. D. Johnson and S. A. Kivelson. This work was supported by the Materials Sciences Division, Office of Science, U.S. Department of Energy under Contract No. DE-AC02-98CH10886, and by the U.S.-Japan Cooperative Research Program on Neutron Scattering. JMT and KN wish to thank the staff of the Laboratoire Léon Brillouin for their hospitality during the experiments.


text:
# Chandra detection of X-ray Absorption Associated with a Damped Ly𝛼 System

## 1 Introduction

The optical-ultraviolet spectra of quasars show narrow absorption lines from a range of heavy elements which originate in the gaseous disks or halos of intervening galaxies (reviewed by Bechtold 2001). The highest column absorbers – the damped Ly-$`\alpha `$ systems and 21-cm absorbers – are thought to originate in the progenitors of large galaxies like the Milky Way at high redshift (Wolfe, Lanzetta, Foltz, & Chaffee 1995; Rao, Turnshek, & Briggs 1995; Storrie-Lombardi & Wolfe 2000 and references therein). At moderate redshift (z$`<`$1) they generally appear associated with spirals or low surface brightness galaxies (Steidel, Dickinson $`\&`$ Persson 1994; Steidel, Pettini, Dickinson $`\&`$ Persson 1994; Steidel et al. 1997; Le Brun, Bergeron, Boisse, & Deharveng 1997; Lane et al. 1998; Boisse et al. 1998). The damped Ly$`\alpha `$ absorbers have been the subject of many detailed studies of the metallicity of the absorbing gas, since the damping profile of Ly-$`\alpha `$ yields an accurate measure of the H I column density, a first step in measuring abundances, and the relatively high columns enable weak, unsaturated transitions of important ions to be detected. For some damped Ly$`\alpha `$ systems, 21-cm absorption can be detected as well; individual 21-cm velocity components can often be deblended more easily than the UV absorption components, allowing kinematics of the parent galaxy to be probed (Wolfe & Briggs 1981; Briggs & Wolfe 1983; Carilli et al. 1996; Lane, Briggs, & Smette 2000)

The detailed abundance pattern of the heavy elements in quasar absorbers reflect the nucleosynthetic history of the parent galaxy, and hence the origin of the gas. The interpretation of the UV absorption spectrum, however, is uncertain, despite very high quality spectra in many cases (Lu et al. 1996; Prochaska & Wolfe 1999; Pettini et al. 2000 and references therein). The problem is disentangling the effect of dust depletion from intrinsic metallicity variations, in the face of large uncertainties due to saturation of the lines of the most abundant elements, and the uncertain ionization corrections for the abundant elements.

In this paper, we report the first results of a program to measure abundances in quasar absorbers in a new way, namely by measuring the photoelectric absorption in the soft X-rays. The advantage of the X-ray absorption is that it will give a measure of the heavy element column density that is essentially independent of ionization, inclusion of an atom in a molecule, or depletion onto grains (Morrison & McCammon 1983; Wilms, Allen & McCray 2000). Here, we describe Chandra ACIS-S observations of two 21-cm and damped Ly$`\alpha `$ absorbers, at $`z=0.312`$ in PKS 1127-145 ($`z_{em}`$=1.17), $`z=1.77`$ in Q 1331+171 ($`z_{em}`$=2.08), and a third candidate damped Ly$`\alpha `$ absorber, Q 0054+1331 ($`z_{em}`$=0.171). The imaging of PKS 1127-145 is presented in a companion paper (Siemiginowska et al. 2001).

## 2 Observations and Analysis

We observed the quasars with a back-side illuminated CCD, ACIS-S (CCD number 7) on board the Chandra X-ray Observatory (Weisskopf et al. 1996; 2000). Details of the observations are listed in Table 1. We read out custom subarrays (1/8 subarray) and moved the target off-axis in order to mitigate $`\mathrm{`}\mathrm{`}`$pile-up”, which distorts the X-ray spectrum. The measured count rates indicate that pile-up is negligible: less than 1$`\%`$ of the source photons are affected (McNamara, 2001). For all objects, the position of the X-ray source agreed with the optical positions to better than 1 $`\mathrm{}`$. The right ascension and declination of the X-ray sources are listed in Table 1.

The data were reduced by the standard pipeline (version R4CU5UPD7.1 with CALDB v.1.7 for Q 0054+1331 and version R4CU5UPD13.3 with CALDB 2.1 for PKS 1127-145 and Q 1331+171). CALDB 2.1 was used to create the instrument responses and PSF files. We used CIAO 2.0 and Sherpa to analyze the data, and cross-checked the fits with XSPEC where feasible.

Source and background spectra were extracted from the event files assuming circular or elliptical regions for the quasars with dmextract and then fit in Sherpa. (Data for Q1331+171 were obtained far off-axis so an elliptical source region was required).

Results of spectral fitting are given in Table 2. We fit each spectrum with a power law, with fixed unredshifted absorption with the Galactic column, and redshifted absorbers, as indicated in Table 2. We did not use data with E$`<`$ 0.4 keV since large uncertainties exist in the calibration for soft energies at this time. We discuss each object in turn.

## 3 Q 0054+144, PHL 909, $`z_{em}=0.171`$

Q 0054+144 is a low redshift QSO, which, prior to the launch of Chandra, was known from $`ROSAT`$ observations to be bright in the X-rays (HEASARC ROSAT archive, Snowdon et al. 1995). It is radio quiet, and has a well-studied host galaxy, which has colors and morphology suggesting that it is an early-type galaxy (Bahcall, Kirhakos & Schneider 1996; McLure et al. 1999). A candidate damped Ly$`\alpha `$ absorber at $`z_{abs}=0.102`$ was suggested by Lanzetta et al. (1995) using $`IUE`$ data. Because Q 0054+144 appeared to be bright enough to allow a high quality grating observation with $`Chandra`$, we obtained a short exposure to measure the spectral energy distribution precisely, particularly in the critical soft X-ray region where the oxygen edges are expected. The results of the $`Chandra`$ observation spectral fits are listed in Table 2, and the best fit power law is shown in Figure 1.

The $`Chandra`$ X-ray spectrum is adequately fit with a power law and no intervening absorption. Subsequently, we became aware of a Hubble Space Telescope GHRS spectrum of Q 0054+144 taken with the G140L grating, which showed that no damped Ly$`\alpha `$ absorption is present at $`z_{abs}=0.102`$. The GHRS spectrum is shown in Figure 2. The most prominent absorption feature is a complex of absorption at the emission line redshift, seen in Ly$`\alpha `$ and N V $`\lambda \lambda `$1238,1242. Si II and Si III are weak or absent, suggesting that the system originates with material intrinsic to the quasar, although O VI is also absent. A low ionization metal absorption line system is present at $`z=0.102`$, but Ly$`\alpha `$ is not damped. Lines at 1354 and 1365 $`\AA `$ are probably Ly$`\alpha `$ lines at $`z=0.114`$ and 0.123, with metal absorption weak or absent.

Despite the lack of an intervening damped Ly$`\alpha `$ absorber or a strong X-ray absorber in the Chandra spectrum, Q0054+144 would be a good candidate for future X-ray grating observations, because of its high count rate, steep X-ray spectrum and associated ($`z_{abs}z_{em}`$) ultraviolet absorption system.

## 4 Q1331+170, $`z_{em}=2.08`$

Prior to our observations, Q 1331+170 had never been observed before in the X-rays, even though optically it is one of the brightest $`z=2`$ quasars known. There is a strong damped Ly$`\alpha `$ and 21-cm absorber at $`z=1.77`$ which contains unusually low ionization material, including C I (Carswell et al. 1979; Wolfe & Davis 1979; Songaila et al. 1994).

The short (3000 second) $`Chandra`$ observation yielded only 198 counts, and a 3-$`\sigma `$ upper limit to NH at $`z_{abs}`$=1.77 of 1.9$`\times `$ 1022 cm-2, assuming solar abundance. This value is larger than the observed N(HI) from Ly$`\alpha `$ of 1.6$`\times `$1021 cm-2. With a longer exposure, X-ray absorption should be detectable.

## 5 PKS 1127-145, $`z_{em}=1.18`$

PKS 1127-145 is a bright, gigahertz-peaked radio source with a VLBI and X-ray jet (Kellerman et al. 1998; Siemiginowska et al. 2001). The properties of the X-ray jet are discussed in the companion paper (Siemiginowska et al. 2001). A strong intervening Mg II and Mg I absorption system at $`z_{abs}=0.312`$ was discovered by Bergeron & Boisse (1991) who also identified two late type galaxies near the quasar line of sight at the redshift of the absorber. The nearer galaxy, separated by 9.6 $`\mathrm{}`$ from the quasar (27 h$`{}_{100}{}^{}{}_{}{}^{1}`$ kpc), has strong \[O II\], \[O III\] and H$`\beta `$ emission lines. An observation with FOS of the galaxy by Deharveng et al. (1995) did not detect Ly$`\alpha `$ emission, and although this galaxy is one of the reddest in the survey of Steidel (1993) it appears to have little extinction. Deharveng et al. (1995) then estimated that the oxygen abundance is approximately 30$`\%`$ solar, based on the nebular emission line ratios.

Lane et al. (1998) discovered H I 21-cm absorption at $`z=0.3127`$ and described the Ly$`\alpha `$ profile of the absorber, based on FOS data, which appears damped. The neutral hydrogen column density derived from the damped Ly$`\alpha `$ profile is N(HI)=(5.1$`\pm `$0.9)$`\times `$ 1021 cm-2. Further, they obtained an optical spectrum of a third galaxy with \[O III\] emission lines, and a separation, 3.9 $`\mathrm{}`$, or 11 h$`{}_{100}{}^{}{}_{}{}^{1}`$ kpc from the quasar. This may be the absorber, or the absorption could arise in tidal debris associated with this group of galaxies.

The X-ray spectrum of PKS 1127-145 shows significant absorption in excess of the Galactic value. Figure 3a shows the data with E$`>`$2 keV fit by a power law and Galactic absorption; the fit extrapolated for E$`<`$2 keV is poor, with systematic residuals. Figure 3b shows the data fit by a power law, Galactic absorption, and intervening absorption at the redshift of the 21-cm absorber, $`z=0.312`$.

The statistics are good enough to show significant features in the residuals. All of the features can be accounted for by calibration uncertainties, that is, they occur at energies where there are sharp features in the effective area curve due to instrumental properties. Improved calibration may therefore change our fits, but we have been conservative in our interpretation of the X-ray data, and we estimate that these changes will only make the conclusions we have drawn stronger. We fit the spectrum with both Morrison & McCammon (1983) and Wilms, Allen & McCray (2000) cross-sections, and obtain similar residuals for both.

If we fix the redshift of the absorber to be $`z=0.312`$ for PKS 1127-145, then the absorber has N(HI)=1.2 $`\pm `$ 0.08 $`\times `$ 1021 cm-2 if the abundance is solar, compared to N(HI) = 5.1 $`\pm `$ 0.9 cm-2 derived from the damped Ly$`\alpha `$ profile (Lane et al. 1998). The X-ray absorption implies that the abundance is about 23$`\%`$ solar, consistent with that found from the nebular emission lines.

With the low spectral resolution of the ACIS-S, the redshift of the X-ray absorption is not certain to be that of the intervening 21-cm absorption. The two parameters of the absorber, redshift and NH are correlated: higher redshift requires higher NH to fit the observed spectrum. Figure 4 shows the confidence levels for the two parameters, where the 1$`\sigma `$ contour has an elongated shape. When we allow the redshift of the absorber to vary, the best fit redshift is $`z0.3`$, but the chi-squared minimum is very shallow, and adequate fits are allowed at other redshifts. Better calibration of ACIS below 0.4 keV is needed to constrain the redshift better. A grating observation would give a definite answer about the location of the absorber, because of the possible detection of the oxygen edge in the high resolution data.

However, the most plausible redshift for the X-ray absorber is $`z=0.312`$. The UV spectrum shows that the $`z=0.312`$ absorber has by far the highest column of H I of any redshift system along the line of sight (Figure 5). There is an absorber at $`z_{abs}z_{em}`$, but it is weak, and has no associated O VI or N V absorption; it is likely a low column density Ly$`\alpha `$ forest cloud.

We can rule out with high confidence an origin for the absorption in diffuse Milky Way interstellar gas. We refit the X-ray data using a power law and $`z=0`$ absorption only, with the column density allowed to vary as a free parameter. The best fit reduced $`\chi ^2`$ = 1.31 (517 degrees of freedom), compared to $`\chi ^2`$ = 1.125 for the case with a $`z_{abs}=`$ 0.312 absorber. Further, the best fit column if all the absorption is at $`z=0`$ is 1.35 $`\pm `$ 0.07 $`\times 10^{21}`$ cm-2, compared to 0.383 $`\pm `$ 0.01 $`\times 10^{21}`$ cm-2 measured from high resolution 21-cm emission data by Murphy et al. (1996). The absorption column required to fit the X-ray data is more than 10$`\sigma `$ greater than the Galactic column, and is unlikely to be at $`z=0`$.

## 6 Limits on Abundances of Individual Elements for the PKS1127-145 Absorber

The spectral resolution of the ACIS-S is insufficient to resolve the individual edges of the elements. Moreover, we expect that most of the opacity is from oxygen, and if the absorption is at $`z=0.321`$, the K-edge of oxygen at rest energy 0.532 keV is redshifted to 0.405 keV, outside the part of the ACIS-S pass-band with reliable calibration. We tried a fit with $`z_{abs}`$=0.312 frozen, and two absorbers, one with He/H only, frozen at the primordial value, and the other with $`\mathrm{`}\mathrm{`}`$processed” metals, including oxygen (c.f. Madesjki et al. 1996). We allowed the total column of each to vary. The results are shown figure 5.

We found that we could fit the data with a model where all the absorption is the result of helium (and NH=$`3.31\pm 0.15\times 10^{21}`$ cm-2), about a factor of 3 greater than the column derived assuming heavy metals are present. It is implausible however that the X-ray absorber would be pure helium and hydrogen, since the $`z_{abs}=`$ 0.312 system shows strong metal absorption in the ultraviolet.

Even if the X-ray absorption opacity is due to helium and hydrogen only, we can still place a limit on the oxygen abundance relative to the H I measured by the damped Ly$`\alpha `$ profile. The 3$`\sigma `$ upper limit to the oxygen abundance is 17$`\%`$ solar. Improved calibration of ACIS-S may change these results significantly. Grating observations with Chandra will be able to measure the oxygen edge directly, and remove the ambiguity in the assumed ratio of the elements heavier than hydrogen.

## 7 Discussion

High quality spectra from ground and from HST have been used to investigate the detailed abundance patterns in damped Ly$`\alpha `$ absorbers (Lauroesch, Truran, Welty & York 1996; Pettini et al. 1999 and references therein). The observed pattern appears to be consistent with the Galactic halo pattern (Population II) with little dust depletion, although dust depleted solar ratios cannot be ruled out (Lu et al. 1996; Prochaska & Wolfe 1999), although uncertainties remain. Saturation of the absorption lines, particularly for the more abundant elements like C, N and O mean that column densities of these elements are difficult to determine (Lu, Sargent & Barlow 1998). Elements condense out of the gas phase onto dust grains, so that the absorption line columns are estimated to be low by factors of 3-5 for lightly refractory elements such as C, N and O, and depleted by factors of $``$1000 or so for Fe, Mg or Si (Jenkins, Savage & Spitzer 1986). The exact ratio depends on details of the grain chemistry, and the particular history of shocks which the material has suffered. Wide variations in the depletion factors are seen empirically in Milky Way ISM clouds (see Tielens 1998, and references therein).

One element which is not depleted onto dust, however, is zinc (Morton 1974; de Boer, Fitzpatrick, & Savage 1985; Sembach, Steidel, Macke, & Meyer 1995). and extensive surveys of the abundance of Zn II relative to hydrogen in damped Ly$`\alpha `$ absorbers (Meyer, Lanzetta, & Wolfe 1995, Pettini et al. 1997). Zinc shares the same nucleosynthetic origin as iron and so Zn/H gives a reliable indication of Fe/H. Iron unfortunately suffers heavy, variable depletion onto dust grains and therefore is impossible to use to measure abundances (deBoer & Lamers 1978).

The signature of Population II abundances is an enhancement of oxygen and other even-proton nuclei (Ne, Mg, Si, S, Ar and Ca) relative to the iron group elements (Fe, Zn, Cr) (Lauroesch et al. 1996; Timmes, Lauroesch & Truran 1995 and references therein). Such an enhancement is observed very clearly in the analysis of the oldest, metal-poor halo stars in the Milky Way where the oxygen group (by which we mean O, Ne, Mg, Si, S, Ar and Ca) is enhanced by a factor of 3 relative to solar ratios of iron and other even-Z iron peak nuclei (Zn, Fe, Cr, Ni). The reason is that oxygen, along with other $`\alpha `$-elements, is produced primarily by massive stars which evolve on short time scales compared with the lower-mass progenitors of Type Ia SNe which produce the iron group elements such as zinc. It is therefore expected that O/H will increase towards the solar value earlier than Zn/H, particularly given the high star-formation rates expected in the early stages of galaxy evolution.

The X-ray absorption detected with Chandra for PKS 1127-145 measures the abundance of the oxygen group elements directly. The X-ray absorbing column measures the metallicity of the gas independent of dust depletion and saturation effects, and is caused by exactly those elements (O, Ne, Mg, Si, S, Ar. Ca) which are expected to be enhanced in Pop II objects (see Morrison and McCammon (1983). By comparing metallicity derived from the X-ray absorbing column with that measured by the Zn II lines (iron group) we can see directly if the damped Ly$`\alpha `$ has a halo-type abundance pattern, or a disk-type abundance pattern. These observations also serve as an important check on the UV absorption line spectroscopy.

The existing HST spectra of PKS 1127-145 unfortunately are inadequate to detect Zn II at the strength predicted by our X-ray measurements. Ultraviolet transitions of S II are also interesting to search for, but fall within the Ly$`\alpha `$ forest; the strongest S II transition at $`\lambda `$1259 is also blended with the much stronger Si II $`\lambda `$1260 line. Future high resolution observtions in the UV may be able to detect S II lines for direct comparison to the oxygen abundance derived from the X-ray data. Spectroscopy to search for Zn II absorption with the echelle on STIS is planned.

ACKNOWLEDGMENTS

We thank the staff of the Chandra X-ray Center for their aid with the planning and analysis of these observations. We thank Greg Madejski for suggesting we consider the relative importance of absorption by helium. JB appreciates financial support from NSF grant AST-9617060 and Chandra Guest Observer grant GO0-1164X from NASA. AS, ME, TA and AD acknowledge support from NASA contract No. NAS8-39073 (ASC). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA, as well as data obtained with the the NASA/ESA Hubble Space Telescope, retrieved from the archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract No. NAS5-26555.


text:
# Proton Stability in Six Dimensions

\[

## Abstract

We show that Lorentz and gauge invariance explain the long proton lifetime within the standard model in six dimensions. The baryon-number violating operators have mass dimension 15 or higher. Upon TeV-scale compactification of the two universal extra dimensions on a square $`T^2/Z_2`$ orbifold, a discrete subgroup of the 6-dimensional Lorentz group continues to forbid dangerous operators.

hep-ph/0107056 YCTP-P5-01
\]

The stability of the proton is one of the most intriguing problems in physics. The current limit on the proton lifetime is $`1.6\times 10^{33}`$ years if the dominant decay mode is $`pe^+\pi ^0`$, and $`1.6\times 10^{25}`$ years independent of the decay mode . The conserved quantity associated with proton stability, the baryon number $`B`$, is the charge corresponding to a global $`U(1)_B`$ symmetry, or one of its discrete subgroups. It is important to explain why $`B`$ is conserved to such a high degree of accuracy. In this letter we demonstrate that the spatial symmetry of two universal extra dimensions, accessible to all the standard model fields, naturally constrains the proton lifetime to be greater than the experimental limits.

The minimal standard model, including operators of dimension four and lower, automatically conserves $`B`$. This theory, however, is valid only up to some scale $`\mathrm{\Lambda }_{\mathrm{SM}}`$ where new physics becomes manifest. The effect of physics at scales above $`\mathrm{\Lambda }_{\mathrm{SM}}`$ is described most generally by higher dimension operators, suppressed by the appropriate power of $`\mathrm{\Lambda }_{\mathrm{SM}}`$. Among these are $`B`$ violating operators of dimension six: $`(\overline{q}_L^cq_L)(\overline{q}_L^cl_L)`$, $`(\overline{d}_R^cu_R)(\overline{q}_L^cl_L)`$, $`(\overline{q}_L^cq_L)(\overline{u}_R^ce_R)`$, $`(\overline{u}_R^cd_R)(\overline{e}_R^cu_R)`$, where only one Lorentz invariant form is displayed for each four-fermion term; the weak-eigenstate quarks, $`q_L`$, $`u_R`$, $`d_R`$, and leptons, $`l_L`$, $`e_R`$, may belong to any generation. The experimental limits on the proton lifetime impose a bound on the coefficients of these operators: $`C_{\mathrm{SM}}/\mathrm{\Lambda }_{\mathrm{SM}}^2\begin{array}{c}\\ <\end{array}10^{24}\mathrm{TeV}^2`$. If $`\mathrm{\Lambda }_{\mathrm{SM}}\begin{array}{c}\\ >\end{array}10^{16}`$ GeV, then the coefficients $`C_{\mathrm{SM}}`$ can safely be as large as order one at that scale. However, the naturalness of the Higgs sector in the standard model restricts $`\mathrm{\Lambda }_{\mathrm{SM}}`$ to lie in the TeV range. The $`C_{\mathrm{SM}}`$’s must therefore be extremely small — below $`10^{23}`$ or so.

Various explanations for the smallness of $`B`$ violation have been proposed. In the minimal supersymmetric standard model, a discrete symmetry is needed to banish $`B`$-violating terms of dimension three and four from the superpotential. The discrete symmetry may be a remnant of a $`U(1)`$ gauge group , or else one could hope that it is preserved by quantum gravity. In the case of large extra dimensions, the smallness of $`C_{\mathrm{SM}}`$ has been suggested to arise from the localization of the quarks and leptons at different points inside a thick brane .

Here we concentrate on the chiral 6-dimensional standard model, in which all standard model fields propagate in two extra dimensions. The current bound on the compactification scale of two universal extra dimensions is $`1/R\begin{array}{c}\\ >\end{array}500`$ GeV , suggesting a rich phenomenology . The standard model in two universal extra dimensions is especially appealing because of the properties of the Lorentz group in six dimensions, as well as the constraints imposed by anomaly cancellation. The 6-dimensional standard model is an effective theory valid up to a scale $`M_s5/R`$, above which the 6-dimensional $`SU(3)_C\times SU(2)_W\times U(1)_Y`$ gauge interactions become non-perturbative.

In six dimensions, the $`SO(1,5)`$ Lorentz symmetry has two irreducible spin-1/2 representations generated by $`\mathrm{\Sigma }^{\alpha \beta }/2`$, where $`\mathrm{\Sigma }^{\alpha \beta }i[\mathrm{\Gamma }^\alpha ,\mathrm{\Gamma }^\beta ]/2`$ in terms of the anticommuting $`8\times 8`$ matrices $`\mathrm{\Gamma }^\alpha `$, $`\alpha =0,1,\mathrm{},5`$. The product $`\mathrm{\Gamma }^0\mathrm{\Gamma }^1\mathrm{}\mathrm{\Gamma }^5`$ has eigenvalues $`\pm 1`$ defining the two 6-dimensional chiralities. The irreducible gauge \[$`SU(3)_C\times SU(2)_W\times U(1)_Y`$\] and gravitational anomalies cancel for only two chirality assignments, up to an overall sign :

$$𝒬_+,𝒰_{},𝒟_{},\{\begin{array}{ccc}_+,& _{},& 𝒩_{},\\ & \mathrm{or}& \\ _{},& _+,& 𝒩_+,\end{array}$$
(1)
where $`𝒬_+,𝒰_{},𝒟_{}`$ and $`_\pm ,_{},𝒩_{}`$ are the 6-dimensional quarks and leptons, respectively, and a generational index is implicit. We assume that only one of these two assignments applies to all fermion generations, since then the $`SU(2)_W`$ global anomaly demands 3 mod 3 generations . Each fermion generation includes a gauge singlet field, $`𝒩_{}`$, such that the gravitational anomaly cancels. The 6-dimensional chiral quark and lepton fields have four components, and if two dimensions are appropriately compactified on an orbifold, then their zero-modes may be identified with the left- and right-handed standard model fermions $`q_L,u_R,d_R,l_L,e_R`$, as well as three right-handed neutrinos $`\nu _R`$. We first describe the constraints on $`B`$ violation from 6-dimensional Lorentz invariance, and then show that the constraints are equally tight even with the lesser symmetry remaining after orbifold compactification.

We begin by observing that $`SO(1,5)`$ has an $`SO(1,3)\times U(1)_{45}`$ subgroup: the combination of the 4-dimensional Lorentz symmetry associated with the $`x^\mu `$ coordinates, $`\mu =0,1,2,3`$, and the symmetry under rotations of the $`x^4`$, $`x^5`$ coordinates. Any chiral 6-dimensional fermion, $`\mathrm{\Psi }_\pm `$, may be decomposed under this subgroup as

$$\mathrm{\Psi }_\pm =\mathrm{\Psi }_{\pm }^{}{}_{L}{}^{}+\mathrm{\Psi }_{\pm }^{}{}_{R}{}^{},$$
(2)
where $`L`$ and $`R`$ are $`SO(1,3)`$ chiralities, defined by the projection operators

$$P_{\pm }^{}{}_{L}{}^{}=P_{}^{}{}_{R}{}^{}=\frac{1}{2}\left(1\mathrm{\Sigma }^{45}\right).$$
(3)
Since $`\mathrm{\Sigma }^{\alpha \beta }/2`$ are the generators of the spin-1/2 representations of $`SO(1,5)`$, the $`SO(1,3)`$-chiral fermions have $`U(1)_{45}`$ charges given by the eigenvalues of $`\mathrm{\Sigma }^{45}/2`$: $`1/2`$ for $`\mathrm{\Psi }_{\pm }^{}{}_{L}{}^{}`$, and $`\pm 1/2`$ for $`\mathrm{\Psi }_{\pm }^{}{}_{R}{}^{}`$. In Table I we list the charges of the quarks and leptons. The covariant derivative, $`D_\alpha `$, transforms as a vector under $`SO(1,5)`$, so that $`D_4+iD_5`$ has $`U(1)_{45}`$ charge $`1`$.

To construct the $`B`$-violating operators in six dimensions, we first note that since all quark fields carry $`B=+1/3`$, $`SU(3)_C`$ gauge invariance demands that for any such operator the number of quark fields (minus the number of charge-conjugated fields) is a multiple of three. Only operators constructed of fields that have zero modes can induce at tree-level the $`B`$-violating processes searched for in experiments so far. From Table I we see that all quark fields containing zero modes have $`U(1)_{45}`$ charge $`1/2`$, so that an operator with $`|\mathrm{\Delta }B|1`$ has a $`U(1)_{45}`$ charge given by $`(3/2)\mathrm{\Delta }B`$ plus the sum of the charges of all lepton fields included in that operator. Hence, all $`B`$-violating nucleon-decay operators allowed by Lorentz invariance in six dimensions must involve at least three quarks and three leptons.

Using this key fact, it is straightforward to find the lowest dimension operators, invariant under $`SO(1,5)`$ and standard-model gauge transformations, that (after compactification) can induce $`B`$-violating nucleon decays. We consider first the $`_+`$ chirality assignment. The operators then first appear at mass-dimension 16 in the 6-dimensional theory, each involving six fermions and one covariant derivative $`\overline{)}D=D_\beta \mathrm{\Gamma }^\beta `$:

$`\left(\overline{}_+𝒟_{}\right)^2\left(\overline{𝒩}_{}\overline{)}D𝒟_{}\right),\left(\overline{}_{}\overline{)}D𝒟_{}\right)\left(\overline{𝒩}_{}\mathrm{\Gamma }^\alpha 𝒟_{}\right)^2,`$ (4)
$`\left(\overline{}_+𝒟_{}\right)\left(\overline{𝒩}_{}𝒬_+\right)\left(\overline{𝒩}_{}\overline{)}D𝒟_{}\right),`$ (5)
$`\left(\overline{𝒩}_{}\overline{)}D𝒰_{}\right)\left(\overline{𝒩}_{}\mathrm{\Gamma }^\alpha 𝒟_{}\right)^2,\left(\overline{𝒩}_{}𝒬_+\right)^2\left(\overline{𝒩}_{}\overline{)}D𝒟_{}\right),`$ (6)
where we have exhibited only the Lorentz covariant bilinears with the smallest number of $`\mathrm{\Gamma }`$ matrices; permutations of $`\overline{)}D`$ and $`\mathrm{\Gamma }^\alpha `$ are also allowed. The $`SO(1,3)`$ chiralities of the zero modes imply that only the terms in (6) which do not contain $`\mathrm{\Gamma }^{4,5}`$ can induce nucleon decay. Each of the above dimension-16 operators enters with a coefficient proportional to $`M_s^{10}`$ leading to a strong (and adequate) suppression of $`B`$ violation after compactification to four dimensions. We estimate their effects after first considering dimension-17 operators, whose contributions to proton decay turn out to dominate.

There is only one proton-decay operator of dimension 17 (modulo Fierz transformations, and the insertion of two or three $`\mathrm{\Sigma }^{\alpha \beta }`$ matrices), invariant under $`SO(1,5)`$ and standard model gauge transformations (in particular $`U(1)_Y`$), that does not involve the gauge singlets $`𝒩_{}`$:

$$𝒪_{17}=\frac{C_{17}}{M_s^{11}}\left(\overline{}_+𝒟_{}\right)^3\stackrel{~}{},$$
(7)
where $`\stackrel{~}{}`$ is the charge-conjugated Higgs doublet in six dimensions. Dimension-17 operators involving $`𝒩_{}`$ have a similar form, and will be discussed below. All other $`B`$-violating nucleon-decay operators require more derivatives, $`\stackrel{~}{}`$ fields, or fermion bilinears.

Upon integrating over the compact dimensions, $`x^4`$, $`x^5`$, restricting the fields to their zero modes, and replacing the zero-mode Higgs doublet by its VEV ($`v_h174`$ GeV), $`𝒪_{17}`$ gives rise in the 4-dimensional theory to the operator

$$\frac{v_hC_{17}}{A_{45}^{5/2}M_s^{11}}\left(\overline{\nu }_Ld_R\right)\left(\overline{l}_Ld_R\right)^2,$$
(8)
where $`A_{45}`$ is the area spanned by $`x^4`$ and $`x^5`$. For the $`T^2/Z_2`$ square orbifold of radius $`R`$ constructed in , $`A_{45}=2\pi ^2R^2`$. This operator is non-vanishing only if the generational indices of the three $`\overline{}_+𝒟_{}`$ bilinears in $`𝒪_{17}`$ are not all identical. Hence, $`𝒪_{17}`$ induces proton decays into $`e^{}\pi ^+\pi ^+\nu \nu `$ or $`\mu ^{}\pi ^+\pi ^+\nu \nu `$ with the two neutrinos belonging to different generations. Final states with more pions or a $`K^+`$ are kinematically suppressed. The proton width is given by:

$$\mathrm{\Gamma }(pe^{}\pi ^+\pi ^+\nu \nu )\frac{v_h^2C_{17}^2}{A_{45}^5M_s^{22}}\left(\frac{m_p}{3}\right)^{11}\mathrm{\Phi }_5F(\pi \pi ),$$
(9)
where $`m_p`$ is the proton mass. If the process occurs predominantly via the decay of the $`d`$ valence quark , whose constituent mass $`m_p/3`$ is used to set the scale in Eq. (9), then the form factor $`F(\pi \pi )`$, which describes the probability for the two pions to be formed, is expected to be of order unity. The kinematical phase-space factor, $`\mathrm{\Phi }_5`$, defined here as the dimensionless quantity that relates the squared amplitude to the decay width, is tiny for a five-body decay. We estimate $`\mathrm{\Phi }_5\begin{array}{c}\\ <\end{array}(4\pi )^7\times O(10^4)`$, where the factors of $`4\pi `$ account for angular integration and the additional suppression is due to integration over the magnitudes of final state momenta. The process $`p\mu ^{}\pi ^+\pi ^+\nu \nu `$ has an even smaller $`\mathrm{\Phi }_5`$. Thus, $`𝒪_{17}`$ yields a finite proton lifetime:

$$\tau _p\frac{10^{35}\mathrm{yr}}{C_{17}^2}\left[\frac{(4\pi )^710^4}{\mathrm{\Phi }_5F(\pi \pi )}\right]\left[\frac{1/R}{0.5\mathrm{TeV}}\right]^{12}\left[\frac{RM_s}{5}\right]^{22},$$
(10)
where the quantities in square brackets are of order one or larger. There are no published experimental limits for five-body proton decays. For comparison, the searches for $`pe^{}\pi ^+\pi ^+`$  set a limit of $`\tau _p>3\times 10^{31}`$ yr. Thus, for $`C_{17}O(1)`$, the proton lifetime is orders of magnitude longer than the current experimental limits.

The $`𝒪_{17}`$ operator also induces $`B`$-violating neutron decays, including $`ne^{}\pi ^+\nu \nu `$, $`\mu ^{}\pi ^+\nu \nu `$. The four-body decays proceed via the fusion of two $`d`$ valence quarks, which is suppressed by their wave function overlap inside the neutron . This is compensated though by the larger constituent mass, $`2m_p/3`$, of the $`dd`$ pair. Compared with $`\mathrm{\Gamma }(pe^{}\pi ^+\pi ^+\nu \nu )`$ given in Eq. (9), the neutron partial width in these modes is enhanced by $`\mathrm{\Phi }_4/\mathrm{\Phi }_5`$, where $`\mathrm{\Phi }_4\begin{array}{c}\\ <\end{array}(4\pi )^5\times O(10^3)`$ is the four-body phase-space factor. Thus, the inverse width of the $`B`$-violating neutron decays induced by $`𝒪_{17}`$ is smaller than $`\tau _p`$ by three orders of magnitude. The experimental limit on the partial mean life of $`ne^{}\pi ^+`$ is $`6.5\times 10^{31}`$ yr . If a dedicated search for $`ne^{}\pi ^+\nu \nu `$ were to yield a comparable limit, then only a weak bound of $`C_{17}\begin{array}{c}\\ <\end{array}O(1)`$ would be imposed, with uncertainties especially due to the sensitive dependence on $`M_s`$ and $`R`$.

These are striking results. $`B`$-violating nucleon decays are adequately suppressed even with the scale of $`B`$ violation in the TeV range, providing only that the coefficient $`C_{17}`$ is not larger than order unity. In fact, $`C_{17}`$ is likely to be significantly smaller than unity because the operator $`𝒪_{17}`$ is composed of chirality flipping bilinears and therefore can be expected to arise with strength proportional to small Yukawa couplings.

We now return to the other operators. For the dimension-16 operators of Eq. (6), the presence of the covariant derivative leads to an additional factor of $`m_p/3`$ in the 4-dimensional amplitude in place of a $`A_{45}^{1/2}v_h/M_s`$ factor. Thus, there is a supression factor of $`10^3`$ in the widths relative to those arising from dimension-17 operators. As we mentioned above, all dimension-17 operators other than $`𝒪_{17}`$ involve at least one $`𝒩_{}`$ field. The tiny neutrino masses suggest that the couplings of its zero mode, $`\nu _R`$, are small. Thus, the coefficients of these operators may naturally be substantially smaller than $`C_{17}`$. Nonetheless, some could be relevant because they allow new decay modes: $`p\pi ^+\nu \nu \nu `$, $`n\pi ^0\nu \nu \nu `$, and other final states with mesons and three neutrinos. \[$`n\nu \nu \nu `$ requires a three-quark fusion  whose small probability is not compensated by the gain in phase space.\] The experimental limits on such modes are of order $`10^{32}`$ yr . There exists an operator of this type with a single $`𝒩_{}`$ field, $`(\overline{}_+𝒟_{})^2(\overline{𝒩}_{}𝒬_+)\stackrel{~}{}`$, and a few with more $`𝒩_{}`$ fields. The bounds on the coefficients of these operators are as loose as those on $`C_{17}`$.

Consider next the $`_{}`$ chirality assignment. It leads to a different set of proton-decay operators, invariant under $`SO(1,5)`$ and gauge transformations: the dominant ones are six-fermion operators of dimension 15,

$$\frac{1}{M_s^9}\left[C_{15}(\overline{}_+^c𝒰_{})+C_{15}^{}(\overline{}_{}^c𝒬_+)\right]\left(\overline{}_{}^c\mathrm{\Gamma }^\alpha 𝒰_{}\right)^2.$$
(11)
Other six-fermion operators involve $`𝒩_+`$ fields, but again, their coefficients may be small. The operators (11) lead to proton decay into three anti-leptons (one or two are positively charged) and a number of mesons. For the dominant decay mode, $`pe^+e^+\pi ^{}\overline{\nu }`$, we estimate

$$\tau _p\frac{10^{26}\mathrm{yr}}{C_{15}^2}\left[\frac{(4\pi )^510^3}{\mathrm{\Phi }_4F(\pi )}\right]\left[\frac{1/R}{0.5\mathrm{TeV}}\right]^{10}\left[\frac{RM_s}{5}\right]^{18},$$
(12)
where $`F(\pi )`$ is a form factor of order unity. Despite the reduced phase space, the process $`p\mu ^+e^+\pi ^{}\overline{\nu }`$ is more constraining due to better data on decays into muons : $`\tau _p/\mathrm{Br}(p\mu ^+X)>10^{31}`$ yr, where Br$`(p\mu ^+X)`$ is the branching fraction for inclusive decays. We derive a constraint, $`C_{15}\begin{array}{c}\\ <\end{array}10^2(R\times 0.5\mathrm{TeV})^5`$, which is rather loose given that $`C_{15}`$ may naturally be as small as some Yukawa couplings. Nonetheless, the prospects for observing $`B`$ violation look better compared with the $`_+`$ chirality assignment. Recall that no known theoretical argument determines which of the two chirality assignments of Eq. (1) is preferred.

We have so far considered proton decay arising from operators that respect the $`SO(1,5)`$ symmetry. But since the compactification of two dimensions breaks $`SO(1,5)`$, including its $`U(1)_{45}`$ subgroup, we must next study whether the reduced symmetry allows operators that could induce proton decay at an unacceptable level. A simple and symmetric choice for the compactification of the $`x^4`$ and $`x^5`$ dimensions is a $`T^2/Z_2`$ orbifold of equal radii . The “torus” $`T^2`$ is taken to be a square of size $`2\pi R`$ in the $`x^4,x^5`$ plane, with periodic boundary conditions, $`\phi (x^4,x^5)=\phi (x^4+2\pi R,x^5)=\phi (x^4,x^5+2\pi R)`$ for any field $`\phi `$. The square has a $`Z_4`$ symmetry, being invariant under $`\pi /2`$ rotations in the $`x^4,x^5`$ plane. Now, since the generator of the $`U(1)_{45}`$ rotations acting on fermions is $`\mathrm{\Sigma }_{45}/2`$ \[the spin-1/2 representation is “double-valued” under $`2\pi `$ rotations\], upon compactification on a square, the $`U(1)_{45}`$ symmetry is broken down to a $`Z_8`$ group whose elements are given by $`U^n`$ with $`U=\mathrm{exp}[i(\pi /2)\mathrm{\Sigma }_{45}/2]`$ and $`n=0,1,\mathrm{},7`$.

Compactification on $`T^2`$ would allow only vectorlike zero-mode fermions in the 4-dimensional theory, so it is necessary to introduce the observed $`SO(1,3)`$ chirality of the standard model fermions by orbifolding the square. Each field is taken to be either even or odd under the $`Z_2`$ orbifold transformation, defined by a $`\pi `$ rotation in the $`x^4,x^5`$ plane. The assigned $`Z_2`$ parity of the quarks and

leptons can be read from Table I: a field has a zero mode if and only if it is even. The fundamental region in the $`x^4,x^5`$ plane of the square orbifold is shown in Fig. I.

The important point is that the square $`T^2/Z_2`$ orbifold does not break further the $`Z_8`$ symmetry of the 4-dimensional effective theory. This is because the $`Z_2`$ orbifold projection commutes with the $`Z_8`$ transformation. A fermion $`\psi `$ of charge $`z`$ under $`Z_8`$ transforms as follows:

$$U\psi (x^4,x^5)=\sigma ^ke^{iz\pi /2}\psi (|x^5|,\sigma x^4),$$
(13)
where $`k=1`$ ($`k=2`$) if $`\psi `$ is odd (even) under $`Z_2`$. The transformation of $`\psi (x^4,x^5)`$ depends on the sign of $`x^5`$, through $`\sigma =+1`$ ($`\sigma =1`$) for $`x^5<0`$ ($`x^50`$). Note also that the presence of orbifold fixed points (see Fig. 1) does not break the $`Z_8`$ symmetry.

Having established that the residual symmetry of the effective theory is $`SO(1,3)\times Z_8`$, where $`Z_8`$ is the subgroup of $`U(1)_{45}`$ defined above, it follows from $`SU(3)_C`$ gauge invariance and from the charges given in Table I that all operators constructed only of zero-mode fields satisfy the selection rule $`\frac{3}{2}\mathrm{\Delta }B\pm \frac{1}{2}\mathrm{\Delta }L=0\mathrm{mod}\mathrm{\hspace{0.17em}4}`$. In particular, no proton-decay operators with less than 6 fermions are allowed and our previous operator analysis remains valid (with the dominant $`B`$-violating effects given by dimension-17 (-15) operators in the $`_+`$ ($`_{}`$) chirality assignment). Furthermore, it is also clear that in this model there are no $`n\overline{n}`$ oscillations (which require a $`\mathrm{\Delta }B=2`$ operator) and no Majorana masses ($`\mathrm{\Delta }L=2`$).

The remarkable suppression of $`B`$ violation established here is due to the $`Z_8`$ symmetry. But if, for example, $`x^4`$ and $`x^5`$ were compactified on a rectangular orbifold ($`T^2/Z_2`$ with different radii), then the $`Z_8`$ would be reduced to a $`Z_4`$. It would follow that dimension-10 operators such as $`(\overline{𝒬}_+^c\mathrm{\Gamma }^\alpha 𝒬_+)(\overline{𝒬}_+^c\mathrm{\Gamma }_\alpha _+)`$, which do not include terms with only zero-modes, would induce proton decay via loops. \[For the $`_{}`$ chirality assignment, the leading $`B`$-violating operators have dimension 12, e.g. $`(\overline{𝒬}_+^c\mathrm{\Gamma }^\alpha 𝒬_+)(\overline{}_{}^c\mathrm{\Gamma }_\alpha 𝒟_{})\stackrel{~}{}`$.\] It is therefore important to explore whether the $`Z_8`$ symmetry is indeed left intact by the compactification. Here we offer only a few comments regarding the choice of vacuum. It has been noted that the square configuration is an extremum of the Casimir energy of bulk fields, including gravity, while the rectangular one is not . Although the lowest order computation indicates that the square configuration is a saddle point of the effective potential (the minimum corresponds to a “rhombus” that has a $`Z_6`$ symmetry, still forbidding $`B`$-violating operators with less than three quarks and three leptons), it could be that a complete computation will reveal this configuration to be a minimum. Alternatively, it is possible to freeze the shape of the torus to be square by a $`Z_4`$ orbifold identification, in which case the $`Z_8`$ symmetry remains unbroken. Either of these possibilities would leave intact our conclusions based on the $`T_2/Z_2`$ square orbifold regarding $`B`$ violation.

We finally mention that some of the fields describing gravitational fluctuations are also charged under the $`Z_8`$ (the components of the metric along $`x^4`$ and $`x^5`$). However, we need not worry about the effects of these fields for nucleon decay if the radius stabilization mechanism makes them heavier than the proton or if they are projected out by a $`Z_4`$ orbifold identification.

In conclusion, we have shown that the combination of standard-model gauge invariance, 6-dimensional Lorentz invariance, and compactification of the two extra dimensions on a $`T^2/Z_2`$ orbifold of equal radii naturally suppresses proton decay to acceptable levels, even with the scale of baryon number violation in the TeV range.

We are grateful to Hsin-Chia Cheng and Erich Poppitz for illuminating conversations. This work was supported by DOE under contract DE-FG02-92ER-40704.