text:
# The impact of shocks on the chemistry of molecular clouds

## 1 Introduction

One of the manifestations of a newly formed low-mass protostar is the presence of a highly collimated and energetic outflow or jet. A natural consequence of the propagation of such high velocity outflows through the protostellar envelope and the ambient molecular medium are shockfronts (Reipurth & Raga, 1999). Shocks both heat and compress the gas and also trigger chemical reactions in the gas-phase, leading to a different chemistry than observed otherwise. Shock processing of dust grains may lead to the injection of atoms and molecules back into the gas, which further distinguishes the chemistry in the shocked region from that of a quiescent protostellar environment (e.g., Tielens, 1999). This paper presents a study of the physics and chemistry of the outflow associated with a well-known class 0 young stellar object NGC~1333-IRAS2 using high-resolution millimeter wavelength interferometer and multi-transition single-dish observations. The NGC~1333-IRAS2 outflow provides a unique opportunity to study the effects of outflows on ambient molecular clouds, as the main shock is well separated from the central protostar and shows a relatively simple morphology. The combination of single-dish and interferometry observations makes it possible to discuss the physical and chemical properties of the outflowing gas and to address the spatial differentiation of the chemistry in the outflow region resolved by the interferometer observations.

Studies of molecular abundances in regions of high outflow activity provide insight into the dependence of the chemical reaction networks on temperature and density. Furthermore it is important to recognize the effect of outflow-triggered chemistry in the inner protostellar envelope, to disentangle it from emission from a circumstellar disk or to address the effect of passive heating by the central protostar. In the central part of the protostellar envelope, thermal evaporation of dust grain mantles can lead to a distinct chemistry as is seen in the case of low-mass protostars (e.g IRAS 16293-2422, Ceccarelli et al., 2000a, b; Schöier et al., 2002; Cazaux et al., 2003).

NGC~1333-IRAS2 (also known as IRAS~03258+3104; hereafter simply IRAS2) is located in the NGC~1333 molecular cloud, harboring several class 0 and I objects, first identified through IRAS maps by Jennings et al. (1987). Continuum observations reveal that IRAS2 is a binary source with two components, IRAS2A and 2B, separated by 6500 AU (30″) (Sandell et al., 1994; Blake, 1996; Looney et al., 2000). IRAS2A is responsible for a highly collimated east-west outflow giving rise to a strong shock $``$15000 AU from the central continuum source (Fig. 1). A strong CO outflow in the north-south direction has also been observed (Liseau et al., 1988; Engargiola & Plambeck, 1999), which originates within a few arcseconds from IRAS2A (Jørgensen et al., 2003a).

Langer et al. (1996) mapped the entire NGC~1333 region in CS and identified two peaks in CS emission toward the IRAS2 outflow. They suggested that these are associated with red-shifted (eastern) and blue-shifted (western) bow shock components of the outflow. Sandell et al. (1994) reported bright CH3OH emission toward the red-shifted (eastern) outflow lobe and correspondingly a high CH3OH abundance, possibly enhanced by the shock. The more detailed structure of the CH3OH emission from the IRAS2A east-west outflow was discussed by Bachiller et al. (1998), who mapped the outflow positions at $``$ 3″ using the IRAM interferometer and 30 m single-dish telescope. Bachiller et al. derived the physical conditions in the shock interaction zone from LVG calculations and obtained a density of $`10^6`$ cm-3 and temperature of $`100`$ K.

Bachiller et al. also found that the observed methanol emission translates to a large enhancement of CH3OH by a factor $`300`$ in the IRAS2A outflow. CH3OH is thought to be released directly from the dust grain mantles and is often seen to be associated with protostellar outflows (e.g., Bachiller et al., 1995). Other often-used tracers of shocks associated with protostellar outflows are Si-bearing species, in particular SiO (e.g., Martin-Pintado et al., 1992; Codella et al., 1999; Garay et al., 2002). High abundances of these species may mark a clear distinction of shocked gas from unprocessed gas in the envelopes around low-mass protostars (Bachiller & Pérez Gutiérrez, 1997; Garay et al., 1998; Bachiller et al., 2001).

In this paper we present a study of the detailed chemistry of the shock associated with the IRAS2A outflow based on observations of a wide range of molecular lines at $`36`$″ resolution from the Owens Valley Radio Observatory (OVRO) and Berkeley Maryland Illinois Association (BIMA) millimeter interferometers, together with millimeter and submillimeter single dish observations from the Onsala 20 m telescope (OSO) and the Caltech Submillimeter Observatory 10.4 m telescope (CSO). Parts of the OVRO observations have previously been presented by Blake (1996). Sect. 2 describes the observations and reductions. The maps from the interferometry observations are presented and discussed in Sect. 3.1, while the single dish observations are treated in Sect. 3.2. The physical and chemical properties of the shock region are analyzed using statistical equilibrium calculations as described in Sect. 4 and molecular abundances are derived. Sect. 5 discusses the inferred chemistry and compares it to other well-studied outflow regions, to other types of star-forming environments and to available models for the chemistry in outflow regions. The main findings are summarized in Sect. 6. A companion paper (Jørgensen et al., 2003a) presents details of a millimeter-wavelength interferometer study of the environment surrounding the central protostellar system.

## 2 Overview of observations

The position of the shock in the eastern lobe of the outflow associated with IRAS2A ($`\alpha (2000)=03^\mathrm{h}29^\mathrm{m}00\stackrel{s}{.}0`$; $`\delta (2000)=31^{}14\mathrm{}19\stackrel{}{.}0`$) was observed with the Millimeter Array of the Owens Valley Radio Observatory (OVRO)111The Owens Valley Millimeter Array is operated by the California Institute of Technology under funding from the US National Science Foundation (grant no. AST-9981546). between October 5, 1994 and January 1, 1995 in the six-antenna L- and H-configurations. Tracks were obtained in two frequency settings at 86 and 97 GHz, and each track observed in alternatingly two fields: the bow shock at the end of the eastern outflow discussed in this paper and the position of the central protostellar source (Jørgensen et al., 2003a). The observed tracks cover projected baselines of 3.1–70 k$`\lambda `$ at 86 GHz and the lines observed are listed in Table 1. The lines were recorded in spectral bands with widths of 32 MHz ($`100`$ km s-1). H13CO+ $`10`$ and CS $`21`$ were observed in 128 spectral channels, the remaining line setups included 64 spectral channels. The complex gain variations were calibrated by observing the nearby quasars 0234+285 and 3C84 approximately every 20 minutes. Fluxes were calibrated by observations of Uranus and Neptune. The rms noise levels are 0.05 Jy beam-1 in the 250 kHz channels with a synthesized beam size of $`3.2\times 2.8`$″. Calibration and flagging of visibilities with clearly deviating amplitudes and/or phases was performed with the MMA reduction package (Scoville et al., 1993).

The millimeter interferometer of the Berkeley-Illinois-Maryland Association (BIMA)222The BIMA array is operated by the Universities of California (Berkeley), Illinois, and Maryland, with support from the National Science Foundation (grants AST-9981308, AST-9981363 and AST-9981289). observed the IRAS2A outflow position between March 4 and April 15, 2003. The array B- and C-configurations provided projected baselines of 2.7–71 k$`\lambda `$. The lines of HCO+ 1–0, HCN 1–0, N2H+ 1–0, and C34S 2–1 were recorded in 256-channel spectral bands with a total width of 6.25 MHz ($`20`$ km s-1). The complex gain of the interferometer was calibrated by observing the bright quasars 3C84 (4.2 Jy) and 0237+288 (2.3 Jy) approximately every 20 minutes. The absolute flux scale was bootstrapped from observations of Uranus. The rms noise levels are 0.2 Jy beam-1 in the 24 kHz channels, with a synthesized beam size of $`6.1^{\prime \prime }\times 5.0^{\prime \prime }`$ FWHM ($`7.6^{\prime \prime }\times 6.8^{\prime \prime }`$ for the C34S and N2H+ observations). The data were calibrated with routines from the MIRIAD software package (Sault et al., 1995).

In addition to the interferometry data, a number of molecular lines were observed toward the position of the red-shifted shock using the Caltech Submillimeter Observatory 10.4 m (CSO)333The Caltech Submillimeter Observatory 10.4 m is operated by Caltech under a contract from the National Science Foundation (grant no. AST-9980846). and Onsala Space Observatory 20 m (OSO)444The Onsala 20 m telescope is operated by the Swedish National Facility for Radio Astronomy, Onsala Space Observatory at Chalmers University of Technology. telescopes. The pointing was checked regularly and found to be accurate to a few arcseconds. The typical beam sizes are 45″–33″ for the OSO 20 m (86–115 GHz) and 26″–20″ for the CSO (217–356 GHz) observations. The data were calibrated using the standard chopper wheel method. The spectra were reduced in a standard way by subtracting baselines and by dividing by the main-beam efficiencies $`\eta _{\mathrm{mb}}`$ as given on the web pages for the two telescopes. $`\eta _{\mathrm{mb}}`$ ranges from 0.6 to 0.43 for frequencies of 86 to 115 GHz for the OSO 20 m and 0.67 to 0.62 for frequencies of 217 to 356 GHz for the CSO. An overview of all the observed lines (single-dish and interferometer) is given in Table 1.

## 3 Data

### 3.1 Interferometry

Fig. 2-3 show moment maps for the lines observed at OVRO (CS, SO, SiO and CH3OH) and BIMA (i.e., HCO+, HCN, N2H+ and C34S). In all maps the coordinates are given as offsets relative to the position of the central protostar, IRAS2A: $`\alpha (2000)=03^\mathrm{h}28^\mathrm{m}55\stackrel{s}{.}7`$, $`\delta (2000)=31^{}14\mathrm{}37\mathrm{}`$. Emission of SO2 and H13CO+ was not detected in the interferometer maps toward the outflow position.

Most of the observed lines show indications of material affected by the outflow with clear line wings spreading out to 10-15 km s-1 from the systemic velocity of 7 km s-1. One exception is N2H+ which shows narrow hyperfine components of approximately 1 km s-1 width (FWHM). Of the observed lines, SiO, CS and HCN show emission stretching furthest from the systemic velocity (out to $``$ 20 km s-1) while the remaining species show somewhat narrower profiles (wings stretching out to $``$ 10 km s-1 relative to the systemic velocity), with HCO+ showing most material closest to the cloud systemic velocity. SiO is not seen at low velocities as is also the case for CH3OH and SO.

In general the molecules observed at BIMA most clearly trace the extended low velocity material. This may in part be due to the different $`(u,v)`$ coverage of the two arrays, with the BIMA observations recovering more of the extended emission, but the species observed at OVRO may also be those that are predominantly enhanced by the outflow generated shock.

The N2H+ emission traces a ridge of material with a number of “cores” stretching from the north-east of the map toward the center and back again to the north-west. A dominant core is seen close to the center (offsets of (67″,-3″)) which is also picked up by the HCO+ maps. South of this core a “$`<`$”-shaped extension is seen in both HCO+ and N2H+. The low velocity CS emission picks up only this feature. A similar component was also seen in the CH3OH emission mapped by Bachiller et al. (1998) but is not evident in the CH3OH observations presented here, possibly due to lower sensitivity.

The high-velocity material is generally much less extended than the low-velocity material. The HCN, SiO and CS trace a narrow component stretching 30-40″ along the outflow propagation and 5-10″ in the direction perpendicular to this. The narrow component points directly to the “$`<`$”-shaped feature in the low velocity emission material. The HCN emission is slightly more extended than that of the two other species, again likely due to the different $`(u,v)`$ coverage of the observations from the two arrays. SO and CH3OH show a slightly weaker structure along the same narrow component.

The emission of CH3OH and SO is located downstream (west) of the outflow propagation direction compared to, e.g., the peak of SiO. Even farther downstream around offset (59″,-13″), HCN, CS and SiO show another strong feature where the N2H+ emission “pinches” the outflow. The HCO+ wing emission is seen only at this position and is found to be more extended, filling out the region void of N2H+. In fact the HCO+ emission can be traced all the way back to the central protostar as is also the case for CO (Fig. 1). It is striking how the HCO+ and, e.g., HCN wing emission trace significantly different components, implying a clear chemical differentiation.

Position-velocity diagrams for CS and SiO are presented in Fig. 5. Note the symmetry around the X-axis in these diagrams with low-velocity emission constituting a broad component of weak emission. For both species, the high velocity component is more pronounced toward the working surface of the outflow.

### 3.2 Single-dish

The observed single-dish spectra are presented in Fig. 6. In agreement with the interferometry maps, the N2H+ hyperfine lines show Gaussian profiles with no sign of outflow wings. The same applies for the C18O 3–2 observations. The 3 mm lines of HCO+, HCN, CS and SO show a “two component” line profile with a narrow peaked profile close to the systemic velocity of the cloud along with a clear red wing extension. SiO on the other hand does not show the narrow component but has a more or less abrupt increase slightly above the cloud rest velocity with a close to linear decline in strength toward the higher (red-shifted) velocities. For the CSO (0.8–1.4 mm) lines a similar trend is seen (except for C18O 3–2). For the CO 2–1 and 3–2 transitions two strong dips are seen at the cloud rest velocity and in the blue part of the spectrum likely due to self-absorption and possible switching onto non-source CO emission. This, however, does not affect the wing emission in the spectrum, which stretches out to $`30`$ km s-1 (20–25 km s-1 relative to the cloud rest velocity).

In Fig. 7 spectra for individual molecules are compared. The two component separation into a core and wing profile seems to be unique for each molecule: the widths of the core part of the lines (where observable in more than one transition) and the dependence of line strength with velocity in the wings, in particular the terminal velocities inferred for each molecule, are both independent of the observed transition. This gives a clear indication that two distinct components with different excitation conditions and chemical properties are observed and that the different molecules probe distinct parts of each of these components, within which the excitation properties do not vary significantly.

As for the interferometry data, the terminal velocities seem to indicate different dynamical components of the outflow region. SiO and CO probe material at the highest velocities relative to the cloud rest velocity of 20–25 km s-1, whereas HCO+ has the lowest degree of wing emission extending $``$ 10 km s-1 from the rest velocity. The remaining molecules fall somewhere in between. A similar trend was seen in the L1157 outflow by Bachiller et al. (2001), who also found CO and SiO to have significantly higher terminal velocities than HCO+, H2CO, SO and CS. They suggested that this could be related to differences in the formation mechanisms for the various species. Bachiller et al. also argued in favor of rather homogeneous excitation conditions in the L1157 outflow, since maps of various transitions for specific molecules were found to be very similar. In this case the main differences among different molecular species would be the result of the chemistry in the region.

### 3.3 Qualitative scenario

The morphology of the interferometry maps and the line-profiles of the single-dish observations can be explained in a simple picture as signatures of a highly collimated outflow ramming into a quiescent core or static cloud traced by the clumps of C18O and N2H+ emission. The action of the outflow leads to sputtering of silicon off dust grains to form SiO, which is totally absent in the quiescent core. At the same time the abundances of CH3OH, CS and SO increase, more likely due to evaporation of ice mantle material. A sequence in chemistry can be seen with SiO probing material at the highest velocities followed by CS and SO, thought to be a result of enhanced sulfur gas phase chemistry, and CH3OH, H2CO and HCN, which are likely to be direct results of grain mantle release. The difference in extent of the line wings indicates either that the molecule formation time scales are varying, with SiO being produced most rapidly in the shock, or that the more volatile species such as CH3OH do not survive at the highest velocities in the flow.

That N2H+ is observed only in the quiescent cloud material is explained if the temperature in the material affected by the outflow increases to $`20`$ K. At this temperature CO is released from grain mantles and becomes the dominant destruction channel of N2H+, lowering the abundance of this molecule. Comparison between the morphology of the N2H+ emission and that of other species indicates a clear interaction between the outflowing material and the ambient cloud. Fig. 2-3 show a cavity of N2H+ emission where the shocked gas (e.g., SiO) appears. The two peaks seen in SiO, HCN and CS also seem to be related to an increase in N2H+ emission. A natural question is whether the outflow shapes or is being shaped by the ambient material. Judging from the morphology of the larger scale emission in Fig. 1, it is striking to note the presence of large amounts of dust northeast of the central IRAS2A protostar. The two CO outflows identified by different authors (e.g., Liseau et al., 1988; Engargiola & Plambeck, 1999; Knee & Sandell, 2000) trace the edges of this dust condensation. The two perpendicular outflows seen toward IRAS2 may therefore reflect the conditions in the ambient cloud material rather than the intrinsic properties of the central protostellar system: the CO outflow could simply be deflected around the dense material traced by the N2H+ and continuum emission leading to the quadrupolar morphology. In either case, however, this does not change the interpretation of the shocked material in this paper. High velocity gas ($``$ 20 km s-1 relative to the systemic velocity) is present toward the eastern lobe as it is seen from most of the species observed in this paper and this indicates the presence of the shock.

## 4 Analysis

### 4.1 Line intensities

Despite the clear separation between the “quiescent” and “shocked” parts of the line profiles for most molecules, the disentanglement of the emission into a core and wing component is not unique. Each line was decomposed in two parts, with the emission integrated over velocities higher and lower (wing and core components, respectively) than where the profiles in Fig. 7 separate. For C18O and N2H+ Gaussians were fitted to each line giving widths of 1.5–2 km s-1 (FWHM), which compares well to the widths of the core parts of the remaining lines. Table 2 lists the resulting line intensities.

The hyperfine splitting of the HCN 1–0 line gives rise to three components within the same setting. The two weaker transitions are offset $`5`$ and 7 km s-1 relative to the main hyperfine line. Overlap in the line wings therefore makes the interpretation of this line difficult. If the emitting material is in local thermodynamical equilibrium (LTE) and the emission is optically thin, one should expect the hyperfine components to be in a ratio of 1:3:5. In Fig. 8 a comparison between the HCO+ and HCN spectra toward the shock position is shown. The HCN spectrum has been overplotted with a composite of three versions of the HCO+ spectrum shifted according to the frequency shifts of the hyperfine lines and scaled in the relative 1:3:5 proportions - leaving an overall “normalization factor” between the HCN and HCO+ line intensities as the only free parameter. The good agreement is remarkable and indicates that HCO+ and HCN trace the same material, especially at the lower velocities. The intensity of the main hyperfine component of the HCN 1–0 line is thereby found to be 0.5 times the intensity of the HCO+ 1–0 line. It should still be re-emphasized, that the interferometer maps show clear differences for the less extended HCN and HCO+ emission at higher velocities.

### 4.2 Tying interferometry and single-dish observations together

One question to address is whether the single-dish and interferometry observations trace the same material. The similar trends seen in the data from the two types of observations seem to support this, but it is also known that the interferometry observations lack sensitivity to extended structures on larger scales.

Fig. 9 presents comparisons between the single-dish observations and spectra taken from interferometry data sets, restored with beam sizes appropriate for the Onsala 20 m telescope. The interferometry spectra have been scaled to match the wings of the single-dish spectra. Only a small correction of the order of 30% needs to be applied to match the CS and SiO spectra, which is close to the calibration uncertainty. This confirms that the higher velocity emission is relatively compact (as the maps also suggest for the emission in the transverse direction of the outflow), making it less subject to the incomplete $`(u,v)`$ sampling of the interferometer observations.

Closer to the systemic velocity of the cloud ($`7`$ km s-1) the discrepancy between the single-dish and interferometry spectra increases. The CS interferometry observations pick up only a small fraction of the emission in the “core” part of the single-dish line. The dip seen in the single-dish CS spectra at the rest velocity of the cloud is a result of self-absorption, while for the interferometry observations it is caused by the interferometer resolving out extended emission close to the cloud systemic velocity. The $`(u,v)`$ sampling is also responsible for the lack of emission in the SiO interferometry spectra at velocities close to the rest velocity, although it is less significant for this molecule.

For the BIMA observations of HCN and HCO+, emission close to the systemic velocity is still mostly resolved out as indicated by the dips in the interferometer HCN and HCO+ spectra and as seen in the channel maps in Fig. 4. The slightly better $`(u,v)`$ coverage from BIMA, however, makes these lines less subject to resolving out at velocities different from the systemic velocity.

Nyquist sampled single-dish maps of the different molecular species would make it possible to combine the interferometer and single-dish data to create maps including the short-spacings. This would settle the issue of the differing $`(u,v)`$ coverage of the two arrays. In this paper we only have single pointing observations. The agreement between the single-dish and interferometry spectra in the line wings, however, justifies the discussion of the outflow component based on the morphology in the interferometry maps presented in Sect. 3.1. The agreement also makes it possible to use the interferometer maps to determine the spatial extent of the wing component, and thereby to estimate the beam filling factor for the single-dish observations. The wing part of the SiO and CS interferometry maps give a rough estimate of the extent of the outflow emission in the transverse direction of 5-10″, leading to filling factors ranging from 0.07 to 0.32 for the single-dish beam sizes. For the core component of the lines on the other hand, the interferometer observations are less useful because of the significant fraction of the low-velocity, extended emission that is resolved out. For the following discussion, a filling factor of unity is therefore assumed for the core component of the single-dish data, which seems realistic as the interferometry maps do reveal emission extended over scales larger than 20–30″.

### 4.3 Statistical equilibrium calculations

In order to derive the physical properties of the emitting gas and the column densities of the various molecules, one-dimensional statistical equilibrium calculations were performed using the 1D radiative transfer code called Radex. Radex uses the escape probability formalism to solve the statistical equilibrium equations for a medium with constant density and temperature (see Jansen et al., 1994; Schöier et al., 2003). For subthermally excited emission this approach is an improvement compared with the rotation diagram method, as demonstrated for interpretations of the CH3OH emission toward the IRAS2A outflow by Bachiller et al. (1998). The statistical equilibrium calculations also treats opacity effects in the correct way, again contrasting the rotation diagram analysis which relies on optically thin emission.

For each molecule, line intensities were calculated for varying column density, density of the main collision partner (H2) and kinetic temperature. The filling factors estimated on the basis of the interferometry observations as discussed above were adopted and the calculated line intensities were compared to the observed ones. Since our main interest is in the relative behavior of the lines some of the uncertainties in the assumptions, e.g., the filling factor of unity for the core component, will cancel out, if the differences between the observed lines are not too large.

The comparison with the observations was performed by calculating the $`\chi ^2`$-statistics for each set of parameters. The uncertainty in the derived line intensities due to the calibration and the disentanglement of the core and wing components was assumed to be 30%. The best fit models for the different species agree quite well in the $`(T_{\mathrm{kin}},n_{\mathrm{H}_2})`$ plane, so all lines are combined into a single $`\chi ^2`$ estimate to constrain the parameters. This is illustrated in Fig. 10 where the constraints on densities and temperatures have been plotted for given (optimal) values of the column densities for the individual molecules. For the core component an H2 density of $`10^6`$ cm-3 and temperature of 20 K is found to be consistent with the observations with a reduced $`\chi _{\mathrm{red}}^2`$ of 1.8 for 8 fitted lines. Note that the density and temperature are closely coupled, making individual determinations somewhat ambiguous, as illustrated in the $`\chi ^2`$ plot, where it is seen that a lower temperature and correspondingly higher density are equally probable. A lower density/higher temperature can also not be completely ruled out. For the wing component a best fit density of 2.0$`\times 10^6`$ cm-3 and temperature of 70 K is found with $`\chi _{\mathrm{red}}^2`$ of 3.8 for 13 fitted lines. The temperature is slightly lower than the value quoted by Bachiller et al. (1998), but still within the mutual uncertainties. The derived column densities assuming these temperatures and densities are given in Table 3.

The kinetic temperatures and densities may vary between regions traced by different molecular species. Both the velocity profiles of the various molecules and the structure of their emission in the interferometry maps show dissimilarities, indicating chemical differentiation possibly due to a combination of the shock evolution and variations in the physical conditions. On the other hand as illustrated in Fig. 7, the conditions for lines of a particular molecule are remarkably homogeneous over the entire shock velocity range. It is also found that the derived parameters do not vary much when specific molecules are included or not. The constraints put on our derived temperature and density are in rough agreement with measurements from other shocked regions from molecular outflows (e.g., Bachiller & Pérez Gutiérrez, 1997; Garay et al., 1998). As a first order approximation the two component structure therefore seems to describe the excitation conditions well. Moreover, the derived column densities are well constrained with the assumed temperature and density in this regime, as illustrated in Fig. 11.

## 5 Discussion

### 5.1 Comparison to other protostellar outflows and envelopes

As mentioned in Sect. 1 only a few studies have addressed the chemistry in outflow regions in detail, the most thorough being that of the L1157 outflow (e.g Bachiller & Pérez Gutiérrez, 1997; Bachiller et al., 2001) and the molecular condensation ahead of the HH2 outflows (Girart et al., 2002). These regions differ significantly in their context: the CO maps L1157 show an outflow progressing through the protostellar envelope and extended cloud. This causes a number of bow shocks and condensations of enhanced density and varying chemistry along the outflow axis. The chemistry in HH2 is thought to be induced by irradiation of the molecular condensation through UV flux from the bright Herbig-Haro object. In such regions especially HCO+ and NH3 should be greatly enhanced according to the models of Viti & Williams (1999).

Table 4 lists the abundances found for the two components of the IRAS2A outflow, calculated as simple ratios between the column densities - taking the CO abundance relative to H2 at constant value of $`1\times 10^4`$. This means that the quoted abundances for the “wing component” are averaged over material with a large range of velocities. As indicated by the maps and the single-dish line profiles the emission from some of the molecules may be more concentrated toward material with lower velocities. In this sense, the quoted abundances are therefore lower limits to the abundances in the regions where these molecules are observed.

Table 4 also compares the derived abundances to those found for the L1157 and BHR71 outflows (Bachiller & Pérez Gutiérrez, 1997; Garay et al., 1998), the molecular condensation ahead of the HH2 object (Girart et al., 2002) and other types of protostellar environments - in particular the IRAS2A protostellar envelope (Jørgensen et al., 2003b), the “hot component” of the IRAS 16293-2422 envelope (Schöier et al., 2002) and the the “C” position of the dark cloud L134N (Dickens et al., 2000). The abundances in the IRAS 16293-2422 envelope were derived through detailed radiative transfer of the dust continuum and molecular line data. The CO, CS and HCN abundances quoted are averages over the entire envelope. Since the abundances in the outer part of the envelope may be lower due to freeze-out (see for example discussion for CO in Jørgensen et al. (2002)), the quoted numbers are likely to be lower limits to the abundances in the inner, warm region of the envelope.

In Fig. 12 the abundances in the two components of the IRAS2A outflow, in the IRAS2A envelope and L1157 outflow are compared. In particular SiO, CH3OH, SO and CS are significantly enhanced by factors 10-104 in the outflow regions compared to the envelope and quiescent dark cloud. The abundances in the IRAS2A wing component and L1157 agree very well for SO and SiO, but the abundances of HCN, H2CO and CS are lower in the IRAS2A outflow by factors of 10-100. The CH3OH abundances are very large for both outflows - between 5 and 10% of that of CO. This is especially significant in comparison with abundances found in the Orion hot core and the low-mass protostar hot core in IRAS 16293-2422. There CH3OH is thought to be enhanced through thermal evaporation off dust grain ice mantles, but its abundance is a factor of 10-100 lower than in the outflow regions. Also the SiO abundances are different between the hot core sources and the outflow regions, whereas the SO abundances are practically identical. For CS and HCN, the abundances have to be significantly higher in the inner regions of IRAS 16293-2422 in order to match the outflow abundances. Compared to the dark cloud, L134N, the CH3OH abundances again stand out as significantly increased, together with CS and SO. HCN and H2CO have similar abundances in the outflow region and the dark cloud whereas HCO+ shows slightly lower abundances in the outflow regions.

### 5.2 Dynamical time scales

The dynamical time-scale of the east-west outflow can be estimated by assuming that: (1) the terminal velocity of the outflow is equal to the maximal radial velocity of the SiO emission, i.e., 29 km s-1, (2) the rest velocity of the cloud corresponds to the velocity of the narrow non-shocked features ($`V_{\mathrm{LSR}}=7`$ km s-1), (3) the distance between the tip of the SiO outflow and the continuum position is the full extent of the outflow and (4) the inclination between the plane of the sky and the outflow is small, as its high degree of collimation and extent plus the relatively low observed velocities seem to indicate. With these assumptions one finds a dynamical age:

$$t_{\mathrm{dyn}}4\times 10^3\times \left(\frac{d}{220\mathrm{pc}}\right)\mathrm{years},$$

which is in good agreement with the dynamical timescale of

$$t_{\mathrm{dyn}}(37)\times 10^3\times \left(\frac{d}{220\mathrm{pc}}\right)\mathrm{years}$$

found by Bachiller et al. (1998) from their CH3OH maps.

The dynamical timescale calculated this way is subject to significant systematic errors. The true extent of the outflow will be $`L_{\mathrm{true}}=L_{\mathrm{obs}}\mathrm{cos}i`$, while the terminal velocity of the outflow relates to the observed maximum radial velocity as $`V_{\mathrm{true}}=V_{\mathrm{obs}}\mathrm{sin}i`$, modifying the dynamical timescale by a factor $`\mathrm{tan}i`$. Furthermore the dynamical timescale at best reflects the properties of the outflow at the present moment - changes in the flow velocities throughout the history will similarly change the outflow dynamical timescale. Thus, it is definitely not an unbiased indicator of the age of the driving protostellar source. Still, it agrees well with the timescales derived from comparison between the envelope structure and collapse models (Jørgensen et al., 2003a) and therefore does give an indication of the order of magnitude of the appropriate timescale to be used when discussing the chemical evolution of the shock in the following section.

### 5.3 Chemical evolution

The above analysis shows that a number of species are significantly enhanced in the outflow region. The various spectral signatures and morphologies in the interferometer maps together with the varying abundances found in the L1157 and IRAS2A outflows indicate different mechanisms regulating the abundances. For the sulfur-bearing species, i.e., CS and SO, the enhancement is thought to occur as a result of enhanced H2S formation as seen in hot cores (Pineau des Forêts et al., 1993; Charnley, 1997). CS and SO show different morphologies in the interferometry maps and different spectral signatures. In particular, CS is present in a narrow high velocity component, while SO only shows up at lower velocities and with more compact emission. In this context the differences and similarities between the abundances found in the L1157 and IRAS2A outflows are also interesting (see Fig. 12): SiO and SO are seen to be quite similar in the L1157 outflow and IRAS2A wing component, while CS, H2CO and HCN are found to be 10 to 100 times more abundant in the L1157 outflow. An explanation could be differences in the abundances of atomic carbon: the production of CS is increased with a higher abundance of atomic carbon, whereas the SO is more closely related to the oxygen abundance (e.g., van Dishoeck & Blake, 1998). In models for gas-phase chemistry, CS, HCN and H2CO are all enhanced with higher abundances of atomic carbon, thus potentially explaining the difference between the IRAS2A and L1157 outflows. The atomic carbon could either be produced in the shock itself or by photodissociation of CO in the pre-shocked gas (e.g., as in the case of IC 443, Keene et al., 1996).

The enhancement of SiO is thought to be caused by atomic silicon sputtering from the surfaces of dust grains and quickly forming SiO in the gas-phase through reactions with OH (e.g., Schilke et al., 1997; Caselli et al., 1997; Pineau des Forêts et al., 1997). CH3OH, H2CO and HCN, in contrast, are most likely enhanced through direct evaporation of ice mantles. Alternative explanations, e.g., gas-phase reactions between CH$`{}_{}{}^{+}{}_{3}{}^{}`$ and H2O forming CH3OH only produce CH3OH abundances of $`15\times 10^8`$ (e.g., Millar et al., 1991), significantly lower than those found in the outflow regions of $`10^610^5`$ (this paper, Bachiller et al., 1995; Bachiller & Pérez Gutiérrez, 1997; Bachiller et al., 1998). The spectral signatures of CH3OH compared to SiO indicate that it is formed at lower velocities. Furthermore, it is seen that CH3OH peaks slightly further downstream compared to, e.g., SiO in the interferometer maps. This is clearly illustrated in Fig. 13 where the intensities of SiO, CH3OH and HCO+ are compared along the outflow propagation axis. Garay et al. (2000) analyzed regions of different shock velocities in the outflow associated with NGC 2071 and found that CH3OH was most prominent in regions with low shock velocities ($`v_s10`$ km s-1). Garay et al. suggested that since molecules such as CH3OH, H2CO and HCN are more volatile, they would only be capable of surviving at much lower shock velocities than would, e.g., SiO. Applied to the IRAS2A outflow this explains both the differences in terminal velocities between SiO on the one hand and HCN, H2CO and CH3OH on the other, but also the different “onsets” along the propagation direction between CH3OH and SiO from the interferometer maps as seen in Fig. 13.

The characteristic core-wing structure of the observed lines is similar to that seen in a SiO survey of protostellar outflows by Codella et al. (1999), who argued that it could be an evolutionary effect with the SiO being produced at high velocities and subsequently slowed down toward lower velocities. They argued that it would take $`10^4`$ years to slow down an outflow-induced shock, which is similar to the time it would take SiO to be destroyed either through direct accretion onto dust grains (Bergin et al., 1998) or through reactions with OH, forming SiO2 (Pineau des Forêts et al., 1997). If this picture applies to the IRAS2A outflow, it is not surprising that SiO has low abundances in the “core component” - or ambient cloud. Since the SiO destruction timescale is similar to the dissipation timescale for the protostellar shock, SiO is almost completely destroyed in the slow-down phase and will therefore not be seen in the low-velocity/quiescent component. On the other hand, since SiO is created as a direct result of the shock impact (Pineau des Forêts et al., 1997; Schilke et al., 1997), it traces the highest velocities in the outflow, together with CO. The characteristic molecular depletion timescale at densities of $`10^6`$ cm-3 typical of the outflow region (Sect. 4.3) is on the order of 103 years, comparable to the outflow dynamical timescale, and could therefore explain why, e.g., SiO and CH3OH are not observed over the entire extent of the outflow back to the central protostar.

The differences between the hot core/warm envelope and outflow abundances of, e.g., CH3OH and SiO could be caused by differing time scales related to the densities in the differing regions: the density in the hot inner part of protostellar envelopes is higher by 2-3 orders of magnitude than what is found in the outflow regions. This will lead to more rapid destruction of molecules with “anomalous” abundances, e.g., SiO, either through accretion or reactions with other species and therefore also lower abundances in the envelope regions. Of course the mechanisms for producing the given molecules in the first place are also likely to be dependent on the environment, further complicating the picture.

The depletion timescale for CH3OH may also be taken as an important clock related to the HCO+ abundance. As noted previously, HCO+ stands out compared to the other molecules tracing material only in the aftermath of the shock. In the L1157 outflow, HCO+ was only found to be prominent in the part of the outflow close to the driving source. Through chemical models, Bergin et al. (1998) found that HCO+ should be destroyed after the passage of the shock through reactions with H2O ($`\mathrm{H}_2\mathrm{O}+\mathrm{HCO}^+\mathrm{H}_3\mathrm{O}^++\mathrm{CO}`$), but increases later as the water abundance reaches lower levels due to freeze-out. This is in fact seen in interferometer data as illustrated in Fig. 13: the emission of HCO+ and CH3OH is almost anticorrelated, with CH3OH being located closer to the “head” of the outflow and HCO+ showing up in the aftermath of the shock. As higher abundances of both CH3OH and H2O are expected to be results of grain mantle evaporation and the timescales for their freeze-out are similar, the HCO+ and CH3OH enhancements should indeed be anticorrelated as seen in Fig. 13.

N2H+, like HCO+, is expected to be destroyed by reactions with H2O. In contrast to HCO+, however, CO may also be important in destroying N2H+ (e.g., Bergin & Langer, 1997; Charnley, 1997). Observational studies of pre- and protostellar objects (e.g., Bergin et al., 2001; Tafalla et al., 2002; Jørgensen et al., 2003a) suggest that N2H+ is enhanced where CO is depleted. The narrowness of the N2H+ lines and the morphology of the emission in the IRAS2A region indicate that this molecule is indeed only tracing the ambient cloud material where CO may be depleted and not the outflowing material where CO is returned to the gas-phase.

## 6 Conclusion

A (sub)millimeter study of the shock associated with the NGC1333-IRAS2A outflow has been presented. Both single-dish and interferometry line observations are presented, which allows for a detailed discussion of both the physical and chemical properties in the shocked region and the spatial distribution of emitting species. The main findings are as follows:

1. Interferometer observations of the outflow region reveal a distinct morphology with a narrow high velocity feature in CS, SiO and HCN, while the low velocity part traces the more spatially extended material. HCO+ is seen only in the aftermath of the shock, whereas N2H+ shows very narrow lines and does not seem to be present in the outflowing gas.
2. Statistical equilibrium calculations show that the region can be divided into two distinct components: low velocity material in the ambient, quiescent cloud and high velocity material associated with the outflow. It is found that the conditions in the quiescent material are consistent with a temperature of $``$20 K and density of $`1\times 10^6\mathrm{cm}^3`$, while the shocked gas is warmer with a temperature of $``$70 K and density of $`2\times 10^6\mathrm{cm}^3`$. Within these components, however, the physical conditions are remarkably homogeneous as indicated by the similarities between the lineshapes for different transitions of various molecules and between the wing emission in spectra from interferometry and single-dish observations.
3. The chemistry in the outflow and quiescent regions are significantly different. CH3OH, SiO and the sulfur-bearing species are significantly enhanced in the outflowing gas, whereas HCN, H2CO, HCO+ and N2H+ show abundances more similar to those found in molecular clouds and protostellar envelopes. Compared to the well-studied L1157 outflow the CS, HCN and H2CO abundances are markedly lower. This could be due to differences in the amount of atomic carbon in the shocked or pre-shocked gas. Higher abundances of the shock-tracing molecules (in particular SiO and CH3OH) in outflows compared to the inner warm envelopes may be related to more rapid destruction of these molecules in the envelopes where the densities are higher.
4. A scenario is suggested where the highly collimated protostellar outflow is progressing into a region with a steep density gradient, focusing the shock and giving rise to the narrow morphology observed. This leads to a shock-induced chemistry, which can explain both the morphologies and qualitatively the abundances of the different molecules. In particular, CH3OH is seen to be greatly enhanced reaching an abundance of about 5% of the observed CO abundance in the shocked gas. CH3OH is, however, not observed to as high velocities as seen in SiO, possibly as a consequence of CH3OH being more volatile. Thus together with, e.g., HCN and H2CO, CH3OH is not able to survive at higher velocities.

This work illustrates the large impact of protostellar outflows in shaping the physical and chemical properties of their parental environment. The combination of high-resolution interferometer observations and single-dish spectra makes it possible to address the physical and chemical conditions in the shocked and ambient gas and to investigate the spatial variation and time-scales characteristic for the shock induced chemistry. So-far only a few shocks have been studied in great chemical detail. Similar systematic studies of a large number of different outflows will allow for a more detailed comparison between outflows and shocks of different velocities and energetics and in different environments. Future observations with facilities such as the SMA, CARMA, and ALMA will allow further studies of the variation of physical and chemical conditions in shocks through high resolution, high sensitivity multi-transition molecular line observations. Also high spatial resolution observations of H2O lines with Herschel-HIFI can confirm the anticorrelation between HCO+ and H2O. All such more detailed observational studies will serve as important starting points for more detailed physical and chemical models for shocks in protostellar environments.

###### Acknowledgements.

We thank the referee for a prompt and well-considered report. The research of JKJ is funded by the Netherlands Research School for Astronomy (NOVA) through a network 2 Ph.D. stipend and research in astrochemistry in Leiden is supported by a Spinoza grant. GAB acknowledges support from the NASA Origins of Solar Systems program. FLS further acknowledges financial support from the Swedish Research Council.


format: text

text:
# 1 Introduction

## 1 Introduction

One of the experimental techniques for fabricating submicron-sized well-defined pattern on solid surfaces is the low energy ion beam sputtering. Among the various ion beam induced morphologies, the formation of periodic ripple structures has attracted much research interest in recent years both for understanding the underlying physical mechanisms  as well as for potential applications in the field of nanotechnology, e.g. ripple pattern might be exploited as a template for quantum dot formation . Although a large number of works on the ripple morphology have been carried out in semiconductor materials, little work has been done on metal surfaces . The group of Valbusa  was pioneer in observing nanoscale-ripples on single crystals of $`Cu`$ and $`Ag`$, where they utilized the presence of Erlich-Schwoebel (ES) barrier at crystalline step edges to develope ripples along energetically favoured crystallographic directions. Since the diffusion-biased ripples are highly sensitive to the substrate temperature, the ripple morphology are found to be unstable at room temperature, because of the low ES barrier heights . Very recently, Sekiba et al.  have shown that in-situ oxidation of the rippled surface immediately after the formation may provide a long term stability at room temperature or higher.

Polycrystalline metal thin films have wide industrial applications in electronic, magnetic, and optical devices . It will, therefore, be interesting to explore the possibility of formation of correlated surface features such as wavelength selected ripple structures, at sub-micometre length scales, in these systems also. In this paper, we report the development of ripple topography at grazing ion beam sputtering on a number of thin metallic films at ambient temperature.

## 2 Experimental

Thin films of $`Au`$, $`Pt`$, $`Ag`$, $`Cu`$ and $`Co`$ were deposited by d. c. magnetron sputtering (Pfeiffer, PLS 500) onto commercially available polished $`Si(100)`$ wafers, previously degreased and cleaned. The base pressure in the deposition chamber was $`2\times 10^6`$ mbar. The film thicknesses were in the range of $`30`$ to $`200`$ nm. These samples were then sputtered with mass analyzed $`527`$ keV $`Ar^+`$ ions in a low-energy ion beam (LEIB) system developed in the laboratory . The angle of ion incidence with respect to the surface normal was varied between $`10^0`$ and $`80^0`$. The beam current through a $`4`$ mm diameter aperture was in the range of $`23`$ $`\mu A`$. The samples were exposed to ion doses in the range $`5\times 10^{15}`$ $`5\times 10^{17}`$ $`ions/cm^2`$, which were measured by a current integrator (Danfysik, model 554) after suppression of the secondary electron emission. The base pressure in the target chamber was less than $`5\times 10^8`$ mbar. The surface morphology of the ion-irradiated samples was examined by a Park Scientific AFM (Auto Probe CP) in the contact mode. All the measurements were carried out in air at room temperature.

## 3 Results and discussion

Atomic force microscopy of the as-deposited films shows that the initial surface topography contains characteristic bump-like structures, which are typical for thin films grown by the process of sputtering . The evolution of surface morphology as a function of the angle of ion incidence for a given ion dose and energy shows first the development of mound structure which tends to grow upto the angle of incidence $`\theta `$ $``$ $`50^0`$. Further increase of the incidence angle shows the beginning of distinct changes of the surface morphology which ultimately ends up with regular ripple structures at grazing ion incidence. Figs. 1(a) - (t) show some selected AFM images of the sputtered $`Co`$, $`Cu`$, $`Ag`$, $`Pt`$ and $`Au`$ films, respectively, at the angles of $`60^0`$, $`70^0`$, and $`80^0`$, because at these angles the gradual morphological transitions are clearly visible and this is a general feature for all the metallic films studied in the present sputtering conditions. At $`60^0`$ weakly pronounced ripples with wave vector parallel to the ion beam, especially in $`Co`$ and $`Cu`$ substrates, appear. At $`70^0`$ the morphology shows the development of arrays of tiny cones aligned along the projection of the ion beam direction. Finally, at $`80^0`$ regular ripple-like surface instability with the wave vector perpendicular to the ion beam direction is developed. In passing we mention that the AFM images of the ripples including that of the other morphological structures are found to be quite reproducible even after several weeks of bombardment.

For the quantitative analysis of the ripple morphology, we have calculated numerically the height-height autocorrelation function $`C(𝐫)=`$ $`[h(𝐫)h(0)]`$, where $`h(𝐫)`$ is the relative surface height at the position $`𝐫`$ and $``$ $``$ denotes an average over all positions and directions. As an illustration, Fig. 2 shows a typical AFM image of the rippled structure on a $`Au`$ film together with the corresponding two-dimensional autocorrelation function. The ripple wavelength $`\mathrm{\Lambda }`$ is defined as the separation between the central peak and the first secondary correlation maximum while taking linear scans of $`C`$. Fig. 3 shows a typical set of data for the ripple wavelength $`\mathrm{\Lambda }`$ as a function of the substrate material when sputtered with the same total ion dose and energy.

Valbusa et al.  showed that sputtering of metal surfaces involves two types of surface instabilities depending on the angle of ion incidence, $`\theta `$, the first one arising from erosion process and the other deriving from anisotropic surface diffusion. The erosion-induced surface instability dominating at grazing incidence $`\theta >\theta _c`$ ($`\theta _c50^070^0`$ ) leads to ripple structures aligned parallel to the ion beam projection, independent on the surface crystallinity or orientation. The erosive regime is supposed to govern by the Bradley and Harper (BH) theory , where the surface height evolution $`h(x,y,t)`$ can be described as

$$\frac{h}{t}=v_0+\gamma \frac{h}{x}+\nu _x\frac{^2h}{x^2}+\nu _y\frac{^2h}{y^2}K^4h+\eta ,$$
(1)
where $`v_0`$ is the surface erosion rate of the flat surface at normal incidence, $`\gamma `$ is related to the derivative of the sputtering yield with respect to the angle of ion incidence, $`\nu _x`$ and $`\nu _y`$ are the effective surface tensions generated by the surface erosion process, the constant $`K`$ is related to the surface diffusion which is activated by different physical processes, namely, thermal or ion beam induced or both  and, finally $`\eta `$ is the noise term associated with the randomness of the bombarding ions.

Eq. (1) can further be extended to crystalline materials by including the effects of anisotropic diffusion in different crystallographic directions as well as the existence of the ES barrier at the step edges . One of the consequences of the presence of the ES barrier is that the ripple structure in metal surfaces could be formed even at normal ion beam incidence ($`\theta `$ $``$ $`0^0`$) by tuning the surface temperature . However, for polycrystalline thin films, the grains are mainly randomly oriented and the grain sizes are usually much smaller compared to the film thickness . For such a system, the existence of ES barrier is improbable because of the lack of well-defined atomic steps at the surface . Therefore, the approximation of isotropic diffusivity as in amorphous materials seems to hold good also in polycrystalline metallic films. The stability of the ripple structure at room temperature indicates that the thermally activated diffusion energy barriers in thin polycrystalline films is comparatively higher than that in monocrystalline metal surfaces. Rossnagel and Robinson  measured the activation energy for adatom surface diffusion on various polycrystalline materials from the Arrhenius plot of the sputter cone spacings for several temperatures. The barrier heights lie typically in the range of $`0.31`$ eV. Similar plots for ripple wavelengths on $`Ag`$ single crystals yields activation energy around $`0.15`$ eV .

At room temperature the surface diffusion is driven by the collisional effects rather than pure thermal effect . The mobility of adatoms is believed to originate from the overlapping collision cascades due to multiple ion impact . Carter and Vishnyakov  proposed to add a ballistic smoothening term of the form $`A(E,\theta )^2h`$ in Eq. (1) in order to account the effect of recoiling-adatom diffusion induced by ion irradiation at a given energy $`E`$. More recently, Makeev et al.  showed that fourth-order derivatives of the surface height function $`h(x,y)`$ may also cause the smoothing effect which, however, does not involve real mass transport. For such a case the wavelength of the ripples, in present experimental conditions, can be derived as $`\mathrm{\Lambda }=2\pi \sqrt{\frac{2D_{yy}}{\nu _y}}`$, where $`D_{yy}`$ is the ion induced smoothing coefficient in the $`y`$ direction as defined in Eq. 52 of Ref. 1. The ratio $`D_{yy}/\nu _y`$ can be estimated from the values of the ion penetration depth $`a`$ and the longitudinal and lateral stragglings $`\sigma `$ and $`\mu `$ , respectively, using the computer code SRIM . Although the experimental data at fixed bombarding ion energy and dose follow the same trend as the theoretically calculated $`\mathrm{\Lambda }`$ in different substrate elements, the theory underestimes substantially the experimental wavelength values (cf., Fig. 3). The discrepancy is due to the fact that the ripple wavelength $`\mathrm{\Lambda }`$ is found to increases with the ion dose $`\varphi `$ as a power law $`\mathrm{\Lambda }\varphi ^n`$ with the exponent $`n=0.53`$, e.g. measured for $`Pt`$ films , whereas the BH model predicts no dependence of the ripple wavelength on ion dose or sputtering time.

In passing we should mention that the ripple formation does not depend on the initial surface topography, i.e. whether the initial surface is atomically flat or rough, the surface is always characterized by ripples at grazing ion beam sputtering. Such a result is due to the effect of the noise term in Eq. (1), where random arrival of the bombarding ions destroys the old surface morphology and generates new morphology which grows with the erosion time .

Finally, we have also bombarded clean $`Si(100)`$ wafers in the angular range $`10^0`$ and $`80^0`$ under identical conditions. Here ripples are formed only at the ion incidence angle of $`60^0`$ and the wave vector of the ripples is found to be parallel to the ion beam direction (Fig. 4). The most interesting observation, however, is that the $`Si`$ ripples are generated at doses $`>`$ $`10^{17}`$ $`ions/cm^2`$ in agreement with that of others, e.g. , in contrast to the metallic ripples which begin to develop at much lower doses of about $`10^{15}`$ $`ions/cm^2`$.

## 4 Conclusion

In conclusion, grazing ion beam sputtering at room temperature can induce ripple structures on thin metallic films similar to that reported on single crystal metal surfaces at low temperature . Unlike the latter, such ripples are found to be quite stable at ambient conditions. The present experiment also indicates that the ripples do not form, and the surface undergoes kinetic roughening so long as $`\theta <\theta _c`$, where $`\theta _c50^060^0`$. Finally, it is noted that, independent of the initial morphology of the surface, the ripples in metallic films start to generate at ion doses as low as $`10^{15}`$ $`ions/cm^2`$ which is roughly two orders of magnitude smaller than that for the ripples formed on $`Si`$ surfaces.

## 5 Acknowledgement

The authors thank Mr. A. Das for technical assistance during the AFM measurements.

## References

* M. A. Makeev, R. Cuerno, and A. -L. Barabási, Nucl. Instrum. and Meth. B 197 (2002) 185.
* G. Li, J. Zhang, L. Yang, Y. Zhang, and L. Zhang, Scripta mater. 44 (2001) 1945.
* U. Valbusa, C. Boragno, and F. Buatier de Mongeot, Materials Science and Engineering C 23 (2003) 201.
* U. Valbusa, C. Boragno, and F. Buatier de Mongeot, J. Phys.: Condens. Matter 14 (2002) 8153.
* D. Sekiba, S. Bertero, R. Buzio, F. Buatier de Mongeot, C. Boragno, and U. Valbusa, Appl. Phys. Lett. 81 (2002) 2632.
* International Technology Roadmap for Semiconductors (Semiconductor Industry Association, San Jose, CA, 2001). \[http://public.itrs.net/Files/2001ITRS/Home.htm\].
* P. Karmakar and D. Ghose, in: Proc. DAE-BRNS Indian Particle Accelerator Conference-2003, Eds. S. C. Bapna, S. C. Joshi and P. R. Hannurkar (Allied Publishers Pvt. Ltd., New Delhi, 2003), p. 169.
* G. S. Bales, R. Bruinsma, E. A. Eklund, R. P. U. Karunasiri, J. Rudnick, and A. Zangwill, Science 249, 264 (1990).
* R. M. Bradley and J. M. E. Harper, J. Vac. Sci. Technol. A 6 (1988) 2390.
* P. Sigmund, J. Mater. Sci. 8 (1973) 1545.
* G. Costantini, Thesis, University of Genova, December 1999.
* R. M. Bradley and J. M. E. Harper, J. Vac. Sci. Technol. A 6 (1988) 2390.
* G. Costantini, S. Rusponi, F. Buatier de Mongeot, C. Boragno, and U. Valbusa, J. Phys.: Condens. Matter 13 (2001) 5875.
* J. H. Jeffries, J. -K. Zuo, and M. M. Craig, Phys. Rev. Lett. 76 (1996) 4931.
* S. Wei, B. Li, T. Fujimoto, and I. Kojima, Phys. Rev. B 58 (1998) 3605.
* S. M. Rossnagel and R. S. Robinson, J. Vac. Sci. Technol. 20 (1982) 195.
* G. Carter, J. Phys. D: Appl. Phys.34 (2001) R1-R22.
* G. K. Wehner, Appl. Phys. Lett. 43 (1983) 366.
* G. Carter and V. Vishnyakov, Phys. Rev. B 54 (1996)17647.
* Z. F. Ziegler, IBM Research, SRIM-2000.40 (PC version), Yorktown Heights, NY, 1999.
* P. Karmakar and D. Ghose, to be published.
* S. Habenicht, K. P. Lieb, J. Koch and A. D. Wieck, Phys. Rev. B 65 (2002)115327.
* G. Carter, V. Vishnyakov and M. J. Nobes, Nucl. Instrum. and Meth. B 115 (1996) 440.

## Figure captions

Fig. 1. AFM images of the unbombarded and $`16.7`$ keV $`Ar^+`$ sputtered $`Co`$, $`Cu`$, $`Ag`$, $`Pt`$ and $`Au`$ surfaces at different angles of incidence $`\theta `$ as indicated. The bombarding dose $`\varphi `$ for $`Co`$, $`Cu`$, $`Ag`$, and $`Au`$ is $`1\times 10^{17}`$ $`ions/cm^2`$, while that for $`Pt`$ is $`5\times 10^{16}`$ $`ions/cm^2`$. The ion beam direction is from the bottom to the top.

Fig. 2. (a) AFM image of the rippled surface on $`Au`$ film after $`16.7`$ keV $`Ar^+`$ ion sputtering at $`\theta =80^0`$ and $`\varphi =1\times 10^{17}`$ $`ions/cm^2`$; the ion beam direction is from the bottom to the top. (b) showing the corresponding 2D-autocorrelation function. (c) showing the 1D-autocorrelation function along the marked line in (b) in order to determine $`\mathrm{\Lambda }`$.

Fig. 3. The ripple wavelength $`\mathrm{\Lambda }`$ versus the substrate element after $`16.7`$ keV $`Ar^+`$ ion sputtering at $`\theta =80^0`$ and $`\varphi =1\times 10^{17}`$ $`ions/cm^2`$.

Fig. 4. AFM image of the ripples formed on a $`Si(100)`$ surface after $`16.7`$ keV $`Ar^+`$ ion sputtering at $`\theta =60^0`$ and $`\varphi =5\times 10^{17}`$ $`ions/cm^2`$. The ion beam direction is from the bottom to the top.


text:
# Quasi-Periodic Changes in the Parsec-Scale Jet of 3C 345

## 1. Introduction

The 16th magnitude quasar 3C 345 ($`z=0.595`$) can be regarded as one of the archetypical sources for studies of superluminal motion with components moving along strongly curved trajectories on pc-scales. Several components with apparent velocities of $`220c`$ have been observed in the jet of 3C 345 (Zensus, Cohen, & Unwin 1995; Lobanov 1996). Periodic outbursts occurring every $`10.1\pm 0.8`$ yr were found by Zhang, Xie, & Bal (1997) from the B-band historical light curve.

The ejection directions of the jet components from the core (the optically thick region) vary, and the component trajectories differ significantly from each other (e.g. Zensus et al. 1995). The curvature of the component trajectories increases towards the core, while the trajectories straighten at larger core separations. This makes the innermost region the most interesting for studying the kinematics of the jet components.

## 2. Observations

Seven VLBA observations at 22 GHz and nine VLBA observations at 43 GHz were made between 1997.40 and 1999.64. Additionally, four space-VLBI (VSOP) observations were performed at 1.6 GHz and four at 5 GHz between 1998.22 and 1999.69. We also reanalysed 22 GHz data from Ros, Zensus, & Lobanov (2000) at epochs 1996.41 and 1996.82 and 43 GHz data from A. Marscher in epochs 1998.15, 1998,41, 1998.58, 1998.76, and 1998.94. Altogether the data base contains 14 epochs at 43 GHz, 9 epochs at 22 GHz, 4 epochs at 5 GHz and 4 epochs at 1.6 GHz. A detailed description of this work can be found in Klare (2003).

The observations allowed us to study in detail 5 jet components (C7–C11) in the inner 3 mas. Three of them were ejected during the time period covered by the new observations and could be followed from the beginning of their evolution, much more closely than it was done for older components. To investigate the long-term behavior of the jet, the new observations were complemented by all data available in the literature, covering the period from 1979 to 2000.

Model fitting by Gaussian components was performed in DIFMAP to parameterise the jet structure and study its evolution. The core of the jet is used as the stationary reference point, and the offset of the jet components is measured with respect to it. The observed proper motion of a component in right ascension and declination is fit by polynomials. To improve the fits, we determined the core shift with respect to compact optically thin jet features, allowing us to combine observations at different frequencies. This procedure improved significantly the accuracy of the trajectory determination for the younger components C7, C8 and C9. This led to a much better determination of their physical jet parameters.

## 3. Similarity and quasi-periodic changes in the jet

There is a clear trend of the trajectories of the jet components C4 to C7 in 3C 345 as shown in Figure 1. C6 was a peculiar faint component and disappeared at about 1 mas core distance. The new data for C8 show for the first time a similarity in shape and curvature of the trajectories between two different components, namely C5 and C8. The younger jet component C8 passed through similar stages such as the turning point in its trajectory 8–10 years later than C5.

Polynomial fits to jet component trajectories for C4, C5, C7, C8 and C9 in 3C 345. Remarkable is the similar shape of the trajectories C5 and C8. The numbers in brackets are epochs at core distances of about 0.25 mas.

The direction of the component ejection (P.A. near the core) changes regularly, with a period of 8–10 years (shown by the dashed line in Figure 2). Additionally, there is an overall trend (shown by the two dot-dashed lines) of $`2.6^{}\text{year}^1`$.

The maxima in flux density of the different jet components as a function of time are shown in Figure 3. A 9.0 yr period is observed here. Comparison of the jet trajectories with the component flux density maxima reveals a correlation. The turning-points in the trajectories coincide with the component flux density maxima. Furthermore, a detailed kinematical analysis of the component trajectories show that the most plausible scenario is that of an increasing Lorentz factor along the trajectories (Klare 2003; Lobanov & Zensus 1996, 1999). From these follows, that the viewing angles of the components decrease shortly after ejection to $`1^{}`$ and turn back to higher values. The component flux density peaks are thus caused by Doppler boosting. This means, that a periodic geometrical change of the jet causes the observed 9.0 yr period in the component flux density maxima.

## 4. Discussion

The results presented above emphasize the existence of an underlying quasi-periodic process in the pc-scale jet of 3C 345 with a time period of about 9 yr. A possible scenario, which can explain the observed quasi-periodicity, is the binary black hole (BBH) system (e.g., Begelman, Blandford, & Rees 1980). The BBH model was applied to describe the optical light curve and the kinematic and flux density evolution of the jet component C7 in 3C 345 (Lobanov & Roland 2002). The model yields an orbital period of 170 yr and a precession period of 2500 yr. The $`2.6^{}\text{year}^1`$ of the initial P.A. shown in Figure 2 may be related either to the orbital or to the precession period in this system. The observed 9 yr period is too short to be related to the orbital period in the BBH system. The 9 yr period is more likely to be caused by the rotation of the jet or the rotation of the accretion disc. Further work has to be done to determine the origin of the 9 yr period.

### Acknowledgments.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We gratefully acknowledge the VSOP Project, which is led by the Japanese Institute of Space and Astronautical Science in cooperation with many organizations and radio telescopes around the world. We thank A. Marscher for sharing with us the VLBA data at 43 GHz in 1998.

## References

Begelman, M.C., Blandford, R.D., Rees, M.J. 1980, Nature, 287, 307

Klare, J. 2003, Ph.D. thesis, Rheinische Friedrich-Wilhelms-Universität Bonn, Germany: http://thesis.jensklare.com

Lobanov, A.P. 1996, Ph.D. thesis, New Mexico Institute of Mining & Technology Socorro, USA

Lobanov, A.P., Zensus, J.A. 1996, in Energy Transport in Galaxies and Quasars, ed. Hardee, P.E., Bridle, A.H., Zensus., J.A., ASP Conf. Series, Vol 100, 109

Lobanov, A.P., Zensus, J.A. 1999, ApJ, 521, 509

Lobanov, A.P., Roland, J. 2002, in Proc. 6th EVN Symp., ed. Ros, E., Porcas, R.W., Lobanov, A.P., Zensus, J.A., (Bonn: MPIfR), 121

Ros, E., Zensus, J.A., Lobanov, A.P. 2000, A&A, 354, 55

Zensus, J.A., Cohen, M.H., Unwin, S.C. 1995, ApJ, 35, 433

Zhang, X, Xie, G.Z., Bal, J.M. 1997, A&A, 330, 469


text:
# Casimir effect for the scalar field under Robin boundary conditions: A functional integral approach

## I Introduction

The Casimir force between two uncharged macroscopic bodies in vacuum is widely regarded as arising from the zero-point fluctuations intrinsic to any quantum system. In the case of two flat parallel plates a distance $`a`$ apart the pressure is proportional to $`a^4`$ . In recent years several groups have performed high-precision measurements of the Casimir force between a flat plate and a spherical surface (lens) or a sphere , and also between two parallel flat plates . In the last case the original Casimir formula was confirmed to $`15\%`$ accuracy.

Due to its fundamental character the Casimir effect has applications in many areas of physics, ranging from the theory of elementary particles and interactions  to atomic and molecular physics . Besides, it has analogues in condensed matter physics, for instance in fluctuation induced forces  and boundary critical phenomena . More recently, its relevance to the design and operation of micro- and nano-scale electromechanical devices has been emphasized .

Usually the details of the interaction between the vacuum fluctuations of the quantum field and the macroscopic bodies are neglected, and replaced by classical boundary conditions (BC) at the boundary of the latter. While Dirichlet and Neumann BC have been extensively studied over the past, the more general case of Robin BC has attracted little attention.

A field $`\varphi `$ is said to obey Robin boundary condition at a surface $`\mathrm{\Sigma }`$ if its normal derivative at a point on $`\mathrm{\Sigma }`$ is proportional to its value there:

$$\frac{}{n}\varphi (x)=c\varphi (x),x\mathrm{\Sigma }.$$
(1)
Neumann and Dirichlet boundary conditions are particular cases of Robin boundary condition: the first one corresponds to $`c=0`$; the other is obtained in the limit $`c\mathrm{}`$ (assuming that $`_n\varphi `$ is bounded).

The mixed case of Dirichlet-Robin (DR) BC were considered in  for a 2D massless scalar field as a phenomenological model for a penetrable surface, with $`c^1`$ playing the role of the finite penetration depth. Recently, the Casimir energy for a scalar field subject to Robin BC on one or two parallel planes was computed in .

Here we also compute the Casimir energy for a scalar field under Robin BC on two parallel planes. However, we introduce a rather different approach which seems more amenable to an eventual computation of radiative corrections in models containing interactions. Its starting point is the introduction of suitable boundary terms in the action, which allows us to compute the partition function of the system without the explicit imposition of Robin BC on the fields. In spite of that we show that the two-point Green function does satisfy those boundary conditions. We also find agreement with the main results of Ref.  for the Casimir energy, which was computed there using the more conventional approach of summing the zero-point energy of the normal modes of the field. In particular, we show that for the mixed case of DR boundary conditions the Casimir energy as a function of $`a`$ develops a minimum, i.e., there is a configuration of stable equilibrium.

Finally, we study the renormalization at one-loop order of the two-point Green function in the $`\lambda \varphi ^4`$ theory submitted to Robin BC on a plate. Our analysis differs from previous ones  in two aspects: (i) we keep $`c`$ arbitrary, instead of considering only the the particular cases $`c=0`$ or $`c=\mathrm{}`$, and (ii) we perform the regularization entirely in momentum space. This procedure avoids dealing with distributions and test functions, which are unavoidable in the mixed coordinate-momentum space regularization used in .

## II The modified action

Let us consider, for simplicity, a real scalar field living in the half-space $`z0`$,Conventions: $`\mathrm{}=c=1`$, $`x=(𝐱,z)`$, where $`𝐱:=(x^0,\mathrm{},x^{d1})`$ and $`z:=x^d`$. satisfying the (Euclidean) equation of motion

$$^2\varphi +U^{}(\varphi )=0$$
(2)
and subject to Robin boundary condition at $`z=0`$,

$$_z\varphi c\varphi |_{z=0}=0.$$
(3)
(We shall assume that $`c0`$, in order to avoid the possible appearance of tachyons in the theory.) One can easily verify that Eqs. (2) and (3) are consequences of the stationary action principle applied to the Euclidean action

$$S[\varphi ]=d^dx\left\{_0^{\mathrm{}}𝑑z\left[\frac{1}{2}(_\mu \varphi )^2+U(\varphi )\right]+\frac{1}{2}c\varphi ^2(𝐱,0)\right\},$$
(4)
where

$$d^dx:=\underset{\beta \mathrm{}}{lim}\underset{L\mathrm{}}{lim}_0^\beta 𝑑x^0_{L/2}^{L/2}𝑑x^1\mathrm{}_{L/2}^{L/2}𝑑x^{d1}.$$
(5)
Indeed, computing $`\delta S:=S[\varphi +\eta ]S[\varphi ]`$ up to second order in $`\eta `$ we obtain

$$\delta S=d^dx\left\{(_z\varphi +c\varphi )\eta |_{z=0}+_0^{\mathrm{}}𝑑z[^2\varphi +U^{}(\varphi )]\eta \right\}+O(\eta ^2),$$
(6)
which implies (2) and (3) if $`\varphi `$ is a stationary point of $`S`$.

Until now, we have been talking about a classical field. What happens when one quantizes the theory? In the usual functional integrals approach one has to integrate over all field configurations obeying certain boundary conditions. If such is the case, is it necessary to retain the surface term in the action? Bordag et al argue that it is, in order to ensure the Hermiticity (more precisely, the self-adjointness) of the fluctuation operator $`\widehat{}:=^2+U^{\prime \prime }(\varphi _\mathrm{c})`$, where $`\varphi _\mathrm{c}`$ is the solution to Eqs. (2) and (3). Saharian  has also argued in favor of such a surface term: without it, the vacuum energy evaluated as the sum of the zero-point energy of each normal mode of the field does not agree with the result obtained by integrating the vacuum energy density.

Here we propose a different approach. We retain the surface term in the action, but we shall not impose any boundary condition at $`z=0`$ on the field configurations to be integrated over. Somewhat surprisingly, if we treat $`U(\varphi )`$ as a perturbation, the two-point Green function of the unperturbed theory does satisfy Robin BC at $`z=0`$, i.e.,

$$(_zc)\varphi (x)\varphi (x^{})_0|_{z=0}=0.$$
(7)
Let us prove it. First of all, we write the partition function of the unperturbed theory as§§§We follow here a procedure very similar to the one employed in Ref.  in the context of quantum field theory at finite temperature.

$$Z_0=[D\varphi _1]_{\varphi (𝐱,0)=\varphi _1(𝐱)}[D\varphi ]\mathrm{exp}(S_0),$$
(8)
where $`S_0`$ is given by Eq. (4) without $`U(\varphi )`$. Note that we are integrating over all field configurations satisfying the boundary condition $`\varphi (𝐱,0)=\varphi _1(𝐱)`$, and then we integrate over all configurations of the surface field $`\varphi _1(𝐱)`$. In other words, we integrate over all possible boundary conditions at $`z=0`$. $`\varphi (𝐱,z)`$ also satisfies periodic BC in the x-coordinates, with period $`\beta `$ in the $`x^0`$-direction and $`L`$ in the others (note, however, that we have the limits $`\beta \mathrm{}`$ and $`L\mathrm{}`$ in mind).

We now decompose $`\varphi `$ as a sum of two fields: $`\varphi =\varphi _0+\eta `$, where $`\varphi _0`$ satisfies

$$^2\varphi _0(x)=0,\varphi _0(𝐱,0)=\varphi _1(𝐱),\varphi _0(𝐱,\mathrm{})=0.$$
(9)
Note that, because of the boundary conditions imposed on $`\varphi _0`$, $`\eta `$ vanishes at $`z=0`$ and $`z=\mathrm{}`$. Eq. (9) can be solved using Fourier transform; the result is

$$\varphi _0(𝐱,z)=\frac{d^dk}{(2\pi )^d}e^{i𝐤𝐱}\phi _1(𝐤)e^{kz},$$
(10)
where $`k=|𝐤|`$ and $`\phi _1(𝐤)`$ is the Fourier transform of $`\varphi _1(𝐱)`$.

In terms of the fields $`\phi _1`$ and $`\eta `$, the partition function $`Z_0`$ becomes the product of two independent functional integrals: $`Z_0=Z_AZ_B`$, where

$$Z_A=[D\phi _1]\mathrm{exp}\left\{\frac{d^dk}{(2\pi )^d}\frac{1}{2}(c+k)\phi _1(𝐤)\phi _1(𝐤)\right\},$$
(11)
$$Z_B=_{\eta (𝐱,0)=0}[D\eta ]\mathrm{exp}\left\{d^dx_0^{\mathrm{}}𝑑z\frac{1}{2}(_\mu \eta )^2\right\}.$$
(12)
Let us now compute the two-point Green function. As a consequence of the factorization of $`Z_0`$, one has $`\varphi (x)\varphi (x^{})_0=\varphi _0(x)\varphi _0(x^{})_A+\eta (x)\eta (x^{})_B`$. According to Eq. (11) we have

$$\phi _1(𝐤)\phi _1(𝐤^{})_A=\frac{1}{c+k}\delta ^{(d)}(𝐤+𝐤^{}).$$
(13)
Combining this result with Eq. (10) we obtain

$$\varphi _0(x)\varphi _0(x^{})_A=\frac{d^dk}{(2\pi )^d}e^{i𝐤(𝐱𝐱^{})}\frac{e^{k(z+z^{})}}{c+k}.$$
(14)
On the other hand,

$$\eta (x)\eta (x^{})_B=\frac{d^dk}{(2\pi )^d}e^{i𝐤(𝐱𝐱^{})}D_\eta (𝐤;z,z^{}),$$
(15)
where $`D_\eta `$ satisfies

$$(_z^2+k^2)D_\eta (𝐤;z,z^{})=\delta (zz^{}),D_\eta (𝐤;0,z^{})=D_\eta (𝐤;\mathrm{},z^{})=0.$$
(16)
One can easily verify that the solution to (16) is given by

$$D_\eta (𝐤;z,z^{})=\frac{1}{k}\mathrm{sinh}(kz_<)\mathrm{exp}(kz_>),$$
(17)
where $`z_<(z_>)=\mathrm{min}(\mathrm{max})\{z,z^{}\}`$. Collecting terms, we finally obtain

$$\varphi (x)\varphi (x^{})_0=\frac{d^dk}{(2\pi )^d}e^{i𝐤(𝐱𝐱^{})}\left[\frac{e^{k(z+z^{})}}{c+k}+\frac{1}{k}\mathrm{sinh}(kz_<)\mathrm{exp}(kz_>)\right].$$
(18)
One can easily verify that $`\varphi (x)\varphi (x^{})_0`$ indeed satisfies Robin BC at $`z=0`$, Eq. (7).

## III Casimir energy

Let us now apply our procedure to the computation of the Casimir energy of a free massless scalar field $`\varphi `$ subject to Robin BC on two parallel plates located at the planes $`z=0`$ and $`z=a`$,

$$_z\varphi c_1\varphi |_{z=0}=0,_z\varphi +c_2\varphi |_{z=a}=0(c_1,c_20).$$
(19)
The Euclidean action for such a system is given by

$$S=d^dx\left\{_0^a𝑑z\left[\frac{1}{2}(_\mu \varphi )^2\right]+\frac{1}{2}c_1\varphi ^2(𝐱,0)+\frac{1}{2}c_2\varphi ^2(𝐱,a)\right\},$$
(20)
and its partition function is given by

$$Z=[D\varphi _1][D\varphi _2]_{\varphi (𝐱,0)=\varphi _1(𝐱)}^{\varphi (𝐱,a)=\varphi _2(𝐱)}[D\varphi ]\mathrm{exp}(S).$$
(21)
The Casimir energy can be extracted from $`Z`$ using the identity

$$E_0=\underset{\beta \mathrm{}}{lim}\frac{1}{\beta }\mathrm{ln}Z.$$
(22)
As we did in the previous Section, we shall write $`\varphi `$ as the sum of two terms: $`\varphi =\varphi _0+\eta `$, where $`\varphi _0`$ is the solution to the classical equation of motion, $`^2\varphi _0=0`$, that obeys the boundary conditions $`\varphi _0(𝐱,0)=\varphi _1(𝐱)`$ and $`\varphi _0(𝐱,a)=\varphi _2(𝐱)`$ (consequently $`\eta (𝐱,0)=\eta (𝐱,a)=0`$). By Fourier transforming in the $`𝐱`$-coordinates one can explicitly solve for $`\varphi _0`$, obtaining

$$\varphi _0(𝐱,z)=\frac{d^dk}{(2\pi )^d}\frac{e^{i𝐤𝐱}}{\mathrm{sinh}ka}\left[\phi _1(𝐤)\mathrm{sinh}k(az)+\phi _2(𝐤)\mathrm{sinh}kz\right]$$
(23)
where $`k=|𝐤|`$ and $`\phi _j(𝐤)`$ is the Fourier transform of $`\varphi _j(𝐱)`$, $`j=1,2`$.

Expressing $`S[\varphi ]`$ in terms of $`\varphi _0`$ and $`\eta `$, we obtain $`S=S_A+S_B`$, where

$$S_A=d^dx\left\{_0^a𝑑z\left[\frac{1}{2}(_\mu \varphi _0)^2\right]+\frac{1}{2}c_1\varphi _0^2(𝐱,0)+\frac{1}{2}c_2\varphi _0^2(𝐱,a)\right\}$$
(24)
$$S_B=d^dx_0^a𝑑z\left[\frac{1}{2}(_\mu \eta )^2\right]$$
(25)
Since $`\varphi _0`$ is a functional solely of $`\varphi _1`$ and $`\varphi _2`$, the partition function $`Z`$ can be written as the product of two terms: $`Z=Z_AZ_B`$, where

$$Z_A=[D\varphi _1][D\varphi _2]\mathrm{exp}(S_A),Z_B=_{\eta (𝐱,0)=0}^{\eta (𝐱,a)=0}[D\eta ]\mathrm{exp}(S_B).$$
(26)
It follows from Eq. (22) that the Casimir energy is given by the sum of two terms, $`E_0=E_A+E_B`$, of which the second one is the Casimir energy of a field subject to Dirichlet boundary conditions on the plates. Since it is a well known result , we shall just quote the result:

$$E_B(a)=\frac{L^{d1}}{a^d}\mathrm{\Gamma }\left(\frac{d+1}{2}\right)(4\pi )^{(d+1)/2}\zeta (d+1).$$
(27)
Let us now compute $`E_A`$. Inserting (23) into (24) we can rewrite $`S_A`$ in terms of $`\phi _1`$ and $`\phi _2`$ as

$$S_A=\frac{1}{2}\frac{d^dk}{(2\pi )^d}\phi _i(𝐤)M_{ij}(𝐤)\phi _j(𝐤),$$
(28)
$$M_{ij}(𝐤)=(c_i+k\mathrm{coth}ka)\delta _{ij}k\mathrm{csch}ka(1\delta _{ij}).$$
(29)
Changing the variables of integration in $`Z_A`$ to $`\phi _j`$, we obtain

$$Z_A=[D\phi _1][D\phi _2]\mathrm{exp}(S_A[\phi _1,\phi _2])=\underset{𝐤}{}\stackrel{1/2}{det}M(𝐤),$$
(30)
hence

$`E_A^{\mathrm{bare}}(a)`$ $`=`$ $`\underset{\beta \mathrm{}}{lim}{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{𝐤}{}}\mathrm{ln}detM(𝐤)`$ (31)
$`=`$ $`{\displaystyle \frac{L^{d1}}{2}}{\displaystyle \frac{d^dk}{(2\pi )^d}\mathrm{ln}\left[c_1c_2+k^2+(c_1+c_2)k\mathrm{coth}ka\right]}.`$ (32)
The expression above diverges, hence requires renormalization. This is achieved by subtracting from it the quantity

$$\delta E_A=\underset{a\mathrm{}}{lim}E_A^{\mathrm{bare}}(a)=\frac{L^{d1}}{2}\frac{d^dk}{(2\pi )^d}\mathrm{ln}[(c_1+k)(c_2+k)],$$
(33)
which can be interpreted as part of the self-energy of the plates. (An analogous subtraction is necessary in the calculation of $`E_B`$.) Since $`\delta E_A`$ does not depend on the distance between the plates, it does not contribute to the force between them. Its subtraction from $`E_A^{\mathrm{bare}}`$ is thus permissible as long as one is interested — as we are — only in the Casimir force. The result of the subtraction is given by

$$E_A(a)=\frac{L^{d1}}{2}\frac{d^dk}{(2\pi )^d}\mathrm{ln}\left[1+\frac{2(c_1+c_2)k}{(c_1+k)(c_2+k)}\frac{1}{e^{2ka}1}\right].$$
(34)
Performing the angular integration and adding the result to $`E_B`$, Eq. (27), we finally obtain the Casimir energy for the massless scalar field under Robin boundary conditions:

$$E_0(c_1,c_2;a)=E_B(a)+\frac{L^{d1}}{(4\pi )^{d/2}\mathrm{\Gamma }\left(\frac{d}{2}\right)}_0^{\mathrm{}}𝑑kk^{d1}\mathrm{ln}\left[1+\frac{2(c_1+c_2)k}{(c_1+k)(c_2+k)}\frac{1}{e^{2ka}1}\right].$$
(35)
As a check of this result, we note that the integral vanishes if $`c_1=c_2\mathrm{}`$ or $`c_1=c_2=0`$, thus reproducing the correct result for Dirichlet-Dirichlet and Neumann-Neumann boundary conditions. Dirichlet-Neumann boundary conditions ($`c_1\mathrm{}`$, $`c_2=0`$) can also be treated exactly: in this case, the integral in Eq. (35) becomes

$$I(a):=_0^{\mathrm{}}𝑑kk^{d1}\mathrm{ln}\mathrm{coth}ka.$$
(36)
Integrating by parts yields

$$I(a)=\frac{2a}{d}_0^{\mathrm{}}𝑑k\frac{k^d}{\mathrm{sinh}2ka}=\frac{4a}{d}\underset{n=0}{\overset{\mathrm{}}{}}_0^{\mathrm{}}𝑑kk^de^{2(2n+1)ka}=\left(2\frac{1}{2^d}\right)\frac{\mathrm{\Gamma }(d)}{(2a)^d}\zeta (d+1).$$
(37)
Inserting this result into Eq. (35) and using the identity $`\mathrm{\Gamma }(2z)=(4\pi )^{1/2}\mathrm{\hspace{0.17em}2}^{2z}\mathrm{\Gamma }(z)\mathrm{\Gamma }(z+1/2)`$  we finally obtain

$$E_0(\mathrm{},0;a)=\frac{L^{d1}}{a^d}\left(1\frac{1}{2^d}\right)\mathrm{\Gamma }\left(\frac{d+1}{2}\right)(4\pi )^{(d+1)/2}\zeta (d+1),$$
(38)
which agrees with the correct result , thus giving us another check of Eq. (35).

Next in simplicity are the following three cases: (i) $`c_1=c_2=c`$, (ii) $`c_1=\mathrm{}`$, $`c_2=c`$, and (iii) $`c_1=0`$, $`c_2=c`$. We shall denote them RR, DR and NR, respectively (R = Robin, D = Dirichlet, and N = Neumann). In all these cases, changing the variable of integration in Eq. (35) to $`q=k/c`$ allows us to rewrite it as

$$E_0^\alpha (c,a)=L^{d1}c^d_\alpha (ca)(\alpha =\mathrm{RR},\mathrm{DR},\mathrm{NR}),$$
(39)
where

$$_\alpha (x)=\frac{\mathrm{\Gamma }\left(\frac{d+1}{2}\right)\zeta (d+1)}{(4\pi )^{(d+1)/2}x^d}+\frac{1}{(4\pi )^{d/2}\mathrm{\Gamma }\left(\frac{d}{2}\right)}_0^{\mathrm{}}𝑑qq^{d1}\mathrm{ln}\left[1+\frac{f_\alpha (q)}{e^{2qx}1}\right],$$
(40)
with

$$f_{\mathrm{RR}}(q)=\frac{4q}{(1+q)^2},f_{\mathrm{DR}}(q)=\frac{2q}{1+q},f_{\mathrm{NR}}(q)=\frac{2}{1+q}.$$
(41)
The graphs of $`_\alpha (ca)`$ in three spatial dimensions are depicted in Fig. 1. We can conclude from it that the Casimir force between the plates: (i) is purely attractive in the RR case (i.e., $`c_1=c_2=c`$); (ii) is repulsive at short distances and attractive at long distances in the DR case, and (iii) is attractive at short distances and repulsive at long distances in the NR case.

To understand those behaviors, let us consider a free field subject to Robin BC at $`z=0`$, i.e., $`_z\varphi (𝐱,0)=c\varphi (𝐱,0)`$. If we write the $`z`$-dependent part of $`\varphi `$ as $`\phi (z)=\mathrm{sin}(kz+\delta )`$, the previous equation becomes $`\mathrm{tan}\delta =k/c`$. It follows that $`\delta 0`$ as $`k0`$, and $`\delta \pi /2`$ as $`k\mathrm{}`$. In terms of $`\phi (z)`$, this is equivalent to say that Robin BC tends to Dirichlet BC at low momentum, and to Neumann BC at high momentum. In the jargon of renormalization group theory, $`c=\mathrm{}`$ is an infrared and $`c=0`$ is an ultraviolet attractive fixed point.

Let us now return to Fig. 1. According to the analysis above, the RR curve should behave as the DD curve in the infrared (i.e., $`a\mathrm{}`$) and as the NN curve in the ultraviolet (i.e., $`a0`$). In both cases, the Casimir force is purely attractive, and so it is in the RR case. The DR curve should behave as the DD curve as $`a\mathrm{}`$, and as the DN curve as $`a0`$; indeed, this is what we observe: attraction at long distances and repulsion at short distances. The analysis of the NR curve is similar.

Such considerations suggest an interesting possibility. Let us suppose that $`0<c_1c_2<\mathrm{}`$. Then $`E_0(c_1,c_2;a)E_0(\mathrm{},\mathrm{};a)`$ as $`a\mathrm{}`$, and $`E_0(c_1,c_2;a)E_0(0,0;a)`$ as $`a0`$. In both these limits, therefore, the Casimir force is expected to be attractive. However, since the crossover from a Dirichlet-like to a Neumann-like BC takes place at different scales for each plate, there could be a range of distances for which one has a Dirichlet-like BC at one plate and a Neumann-like BC at the other, thus leading to a repulsive force between them. That such a possibility can indeed occur is shown in Fig. 2.

## IV Interacting field

There is still another reason for the introduction of boundary terms in the action: the perturbative treatment of a renormalizable interacting theory requires boundary counterterms in addition to the usual ones . Thus, from a conceptual point of view, it is more natural to treat boundary conditions as resulting from the interaction of the field with a background than to impose them a priori. We shall illustrate this point with the calculation of the first order correction to the two-point Green function of the $`\lambda \varphi ^4`$ theory in the presence of a flat boundary, where the field is submitted to Robin BC.

The renormalized Euclidean Lagrangian density of the theory is given byThe factor of $`\frac{1}{2}`$ in the boundary term in the action (4) is due to the fact that $`_0^{\mathrm{}}\delta (z)f(z)𝑑z=\frac{1}{2}f(0)`$.

$$=\frac{1}{2}(_\mu \varphi )^2+\frac{1}{2}m^2\varphi ^2+c\delta (z)\varphi ^2+\frac{\lambda }{4!}\varphi ^4+_{\mathrm{ct}},$$
(42)
with the field $`\varphi `$ living in the half-space $`z0`$. For simplicity, we shall assume $`m=0`$ and $`c0`$. $`_{\mathrm{ct}}`$ contains the renormalization counterterms — the usual ones,

$$_{\mathrm{ct},\mathrm{bulk}}=\frac{\delta Z}{2}(_\mu \varphi )^2+\frac{\delta m^2}{2}\varphi ^2+\frac{\delta \lambda }{4!}\varphi ^4,$$
(43)
plus boundary conterterms, which we shall exhibit later.

As we have seen in Sec. II, the unperturbed two-point Green function (Feynman propagator) is given by \[see Eq. (18)\]

$$G_0(x,x^{}):=\varphi (x)\varphi (x^{})_0=\frac{d^dk}{(2\pi )^d}e^{i𝐤(𝐱𝐱^{})}𝒢_0(k;z,z^{}),$$
(44)
where, in the massless theory,

$$𝒢_0(k;z,z^{})=\frac{e^{k(z+z^{})}}{c+k}+\frac{1}{k}\mathrm{sinh}(kz_<)\mathrm{exp}(kz_>).$$
(45)
If we neglect boundary counterterms, the first order correction to $`𝒢_0`$ is given by

$$𝒢_1(k;z,z^{})=_0^{\mathrm{}}𝑑w𝒢_0(k;z,w)\mathrm{\Sigma }_1(w)𝒢_0(k;w,z^{}),$$
(46)
where the one-loop self-energy, including the mass counterterm, is given by

$`\mathrm{\Sigma }_1(w)`$ $`=`$ $`{\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{d^dq}{(2\pi )^d}𝒢_0(q;w,w)}+\delta m^2`$ (47)
$`=`$ $`{\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{d^dq}{(2\pi )^d}\left(\frac{e^{2qw}}{c+q}+\frac{1e^{2qw}}{2q}\right)}+\delta m^2.`$ (48)
(An ultraviolet cutoff $`\mathrm{\Lambda }`$ is implicit in the integral above and wherever necessary below.) To fix $`\delta m^2`$ we impose that $`\mathrm{\Sigma }_1(w)`$ be finite for $`w>0`$ and vanishes for $`w\mathrm{}`$ (the latter condition is equivalent to requiring that the physical mass is zero when one is infinitely far from the plate). These conditions imply

$$\delta m^2=\frac{\lambda }{2}\frac{d^dq}{(2\pi )^d}\frac{1}{2q},$$
(49)
hence

$$\mathrm{\Sigma }_1(w)=\frac{\lambda }{2}\frac{d^dq}{(2\pi )^d}\left(\frac{1}{c+q}\frac{1}{2q}\right)e^{2qw}.$$
(50)
Although $`\mathrm{\Sigma }_1(w)`$ is now a well defined function of $`w`$ for $`w>0`$, this is not enough to ensure the finiteness of $`𝒢_1`$. Indeed, $`\mathrm{\Sigma }_1(w)w^{(d1)}`$ for $`w0`$, causing the divergence of the integral in Eq. (46) for $`d2`$. This problem can be (partially) solved by the boundary counterterm

$$_{\mathrm{ct},\mathrm{boundary}}^{(1)}=\delta c\delta (z)\varphi ^2.$$
(51)
It adds a new term to the self-energy $`\mathrm{\Sigma }_1(w)`$, turning it into $`\stackrel{~}{\mathrm{\Sigma }}_1(w)=\mathrm{\Sigma }_1(w)+2\delta c\delta (w)`$. Equation (46) then gets replaced by

$$𝒢_1(k;z,z^{})=_0^{\mathrm{}}𝑑w𝒢_0(k;z,w)\mathrm{\Sigma }_1(w)𝒢_0(k;w,z^{})\delta c𝒢_0(k;z,0)𝒢_0(k;0,z^{}).$$
(52)
Let us assume that $`0<z<z^{}`$; then Eq. (52) yields

$$𝒢_1(k;z,z^{})=\delta c\frac{e^{k(z+z^{})}}{(c+k)^2}\frac{\lambda }{2}\frac{d^dq}{(2\pi )^d}\left(\frac{1}{c+q}\frac{1}{2q}\right)[I(0,z)+I(z,z^{})+I(z^{},\mathrm{})],$$
(53)
where

$$I(a,b):=_a^b𝑑we^{2qw}𝒢_0(k;z,w)𝒢_0(k;w,z^{}).$$
(54)
Since we are interested in the large-$`q`$ behaviour of the integrand, we may restrict our attention to $`I(0,z)`$$`I(z,z^{})`$ and $`I(z^{},\mathrm{})`$ are exponentially supressed in that limit. Performing the integration, we obtain

$`I(0,z)`$ $`=`$ $`e^{k(z+z^{})}{\displaystyle _0^z}𝑑we^{2qw}\left[{\displaystyle \frac{e^{kw}}{c+k}}+{\displaystyle \frac{1}{k}}\mathrm{sinh}(kw)\right]^2`$ (55)
$`\stackrel{q\mathrm{}}{}`$ $`{\displaystyle \frac{e^{k(z+z^{})}}{(c+k)^2}}\left[{\displaystyle \frac{1}{2q}}+{\displaystyle \frac{c}{2q^2}}+O\left({\displaystyle \frac{1}{q^3}}\right)\right].`$ (56)
If we further multiply this result by the term in parentheses in Eq. (53), we conclude that the integrand in that equation behaves for large $`q`$ as

$$\frac{e^{k(z+z^{})}}{(c+k)^2}\left[\frac{1}{4q^2}\frac{c}{4q^3}+O\left(\frac{1}{q^4}\right)\right]\stackrel{q\mathrm{}}{}\frac{e^{k(z+z^{})}}{(c+k)^2}\left[\frac{1}{4q(q+c)}+O\left(\frac{1}{q^4}\right)\right].$$
(57)
Therefore, if we choose

$$\delta c=\frac{\lambda }{8}\frac{d^dq}{(2\pi )^d}\frac{1}{q(q+c)}+\delta \overline{c},$$
(58)
with $`\delta \overline{c}`$ finite, $`𝒢_1(k;z,z^{})`$ becomes finite for $`d<4`$.

To see why the boundary conterterm (51) only solves the problem partially, let us take $`z=z^{}=0`$ in Eq. (52). A straightforward calculation then shows that

$`𝒢_1(k;0,0)`$ $`=`$ $`{\displaystyle \frac{1}{(c+k)^2}}\left[\delta c+{\displaystyle \frac{\lambda }{2}}{\displaystyle \frac{d^dq}{(2\pi )^d}\left(\frac{1}{c+q}\frac{1}{2q}\right)\frac{1}{2(q+k)}}\right]`$ (59)
$`=`$ $`{\displaystyle \frac{\delta \overline{c}}{(c+k)^2}}+{\displaystyle \frac{\lambda }{8(c+k)}}{\displaystyle \frac{d^dq}{(2\pi )^d}\frac{1}{q(q+c)(q+k)}}.`$ (60)
The result is finite for $`d=2`$, but diverges for $`d=3`$. In this case, another boundary counterterm is needed. As we show below, it is given by

$$_{\mathrm{ct},\mathrm{boundary}}^{(2)}=\delta b\delta (z)\varphi (_nc)\varphi ,$$
(61)
with $`_n`$ the interior normal derivative (in our case, $`_n=_z`$). Indeed, such a counterterm gives rise to an extra contribution to $`𝒢_1(k;0,0)`$, given by

$`\mathrm{\Delta }𝒢_1(k;0,0)`$ $`=`$ $`\delta b{\displaystyle _0^{\mathrm{}}}𝑑w\delta (w)(_w2c)𝒢_0^2(k;0,w)`$ (62)
$`=`$ $`\delta b{\displaystyle _0^{\mathrm{}}}𝑑w\delta (w)(_w2c){\displaystyle \frac{e^{2kw}}{(c+k)^2}}`$ (63)
$`=`$ $`{\displaystyle \frac{\delta b}{c+k}}.`$ (64)
If we choose

$$\delta b=\frac{\lambda }{8}\frac{d^dq}{(2\pi )^d}\frac{1}{q^2(q+c)}+\delta \overline{b},$$
(65)
with $`\delta \overline{b}`$ finite, we cancel the divergence in $`𝒢_1(k;0,0)`$ for $`d=3`$. Indeed,

$`𝒢_1(k;0,0)+\mathrm{\Delta }𝒢_1(k;0,0)`$ $`=`$ $`{\displaystyle \frac{1}{c+k}}\left[\delta \overline{b}{\displaystyle \frac{\delta \overline{c}}{c+k}}{\displaystyle \frac{\lambda k}{8}}{\displaystyle \frac{d^dq}{(2\pi )^d}\frac{1}{q^2(q+c)(q+k)}}\right]`$ (66)
$`=`$ $`{\displaystyle \frac{1}{c+k}}\left[\delta \overline{b}{\displaystyle \frac{\delta \overline{c}}{c+k}}+{\displaystyle \frac{\lambda }{16\pi ^2}}{\displaystyle \frac{k}{ck}}\mathrm{ln}\left({\displaystyle \frac{k}{c}}\right)\right](d=3).`$ (67)
The boundary counterterm (61) is ineffective if $`z,z^{}>0`$; in this case, $`\mathrm{\Delta }𝒢_1(k;z,z^{})`$ is identically zero. (As a consequence, $`𝒢_1(k;z,z^{})`$ is insentive to the choice of $`\delta \overline{b}`$ if $`z,z^{}>0`$.) On the other hand, that counterterm is necessary if $`z=0`$ and $`z^{}>0`$ (or vice-versa); in this case, it is possible to show that the same choice (65) for $`\delta b`$ also ensures the finiteness of $`𝒢_1(k;z,z^{})`$.

In order to fix the finite part of the boundary counterterms, $`\delta \overline{b}`$ and $`\delta \overline{c}`$, one has to impose a pair of renormalization conditions. A natural choice, for it is satisfied at tree level, is given by

$$𝒢(\kappa ;0,0)=(c+\kappa )^1,$$
(68)
$$\frac{d}{dk}𝒢(k;0,0)|_{k=\kappa }=(c+\kappa )^2,$$
(69)
where $`\kappa `$ is some arbitrary, but nonzero, mass scale.

We end this section with some remarks:

(A) Because of the boundary counterterms, the renormalized two-point Green function does not satisfy Robin BC at $`z=0`$ — except for $`c\mathrm{}`$ (Dirichlet BC); in this case, the boundary condition is preserved at each order in perturbation theory .

(B) Inclusion of another plate at $`z=a`$ does not affect significantly the overall picture. New ultraviolet divergences arise, but the theory is made finite with the same type of counterterms used in the case of a single plate.

## V Conclusions

In this work we have computed the Casimir energy of a free massless scalar field subject to independent Robin boundary conditions on two parallel plates in $`d`$ spatial dimensions. It was shown that for mixed Dirichlet-Robin BC the Casimir energy as a function of the distance $`a`$ between the plates displays a minimum. We managed to understand the behavior of the Casimir energy as a function of $`c`$ relying on an analogy with the renormalization group flows expected from infrared/ultraviolet fixed points. It was found that the Dirichlet BC ($`c\mathrm{}`$) is analogous to an attractive infrared fixed point whereas the Neumann BC ($`c=0`$) resembles an attractive ultraviolet fixed point. This interpretation is consistent with the numerical results shown in Fig. 1 for the Casimir energy as a function of $`ca`$, and suggests that a crossover from Dirichlet-like to Neumann-like behaviors at the plates may lead to a repulsive force between them.

We also provided a detailed analysis of renormalization for the two-point Green function $`G(x,x^{})`$ at first order in $`\lambda `$ in the $`\lambda \varphi ^4`$ theory. For simplicity we worked out the case of Robin BC at a single flat boundary. We have shown that, in addition to the usual “bulk” counterterms, one boundary counterterm is necessary to render $`G`$ finite in the bulk, and a second one is necessary if at least one of its arguments lies on the boundary. This analysis is a necessary step in the computation of radiative corrections to the Casimir energy, which we intend to present elsewhere.

In view of the widespread interest on the Casimir effect and its possible technological applications it is worthwhile to seek alternative computational tools which may go beyond or complement the existing ones. For instance, in  a resummation scheme to compute the leading radiative corrections to the Casimir energy was suggested. On the other hand, the method outlined in the present work seems well suited to the computation of radiative corrections either in a perturbative setting or eventually via application of semi-classical methods. Its starting point is the indirect implementation of the BC by means of appropriate terms in the action functional, which is then reexpressed in terms of two kinds of fields: a field $`\eta (𝐱,z)`$ satisfying Dirichlet BC at the surfaces, and two surface fields $`\varphi _j(𝐱)`$ localized on the planes, depending only on the remaining $`d`$ transverse coordinates. The advantage of this procedure is that the functional integration over the surface fields is unconstrained, i.e., one does not have to enforce explicitly the Robin BC on the fields.

The implementation of BC via local terms in the action is usually employed in studies of boundary critical phenomena . In that context, it can be shown that Dirichlet and Neumann BC correspond to the so-called ordinary ($`c\mathrm{}`$) and special transitions ($`c=0`$), respectively. The Robin BC is relevant in the study of the crossover between those universality classes, for which, however, the computations become much more involved. It is relevant also for the analysis of the ordinary transition; in this case, however, one may resort to an expansion in powers of $`c^1`$ . We expect that the methods proposed here can be useful to the study of the crossover for the relevant case of two flat planes. This is presently under investigation.

###### Acknowledgements.

L.C.A. would like to thank the Mathematical Physics Department of USP at São Paulo for its kind hospitality. L.C.A. is partially supported by CNPq, grant 307843/2003-3. R.M.C. is supported by CNPq and FAPERJ.