text:
# Quantum Conductance Steps in Solutions of Multiwalled Carbon Nanotubes

## Abstract

We have prepared solutions of multiwalled carbon nanotubes in Aroclor 1254, a mixture of polychlorinated biphenyls. The solutions are stable at room temperature. Transport measurements were performed using a scanning–tunneling probe on a sample prepared by spin–coating of the solution on gold substrates. Conductance steps were clearly seen. An histogram of a high number of traces shows maximum peaks at integer values of the conductance quantum $`G_0=2e^2/h`$, demonstrating ballistic transport at room temperature along the carbon nanotube over distances longer than $`1.4\mu m`$.

carbon nanotubes, quantum conductance, ballistic transport

Since their discovery Iijima (1991) theoretical calculations and experimental measurements suggest the enormous potential of carbon nanotubes as the building blocks of future molecular nanoelectronic devices Saito et al. (1998); Venema (1999); Chico et al. (1996); Pérez-Garrido and Urbina (2002). The electronic properties of carbon nanotubes depend strongly on their chirality and diameter. All armchair $`(n,n)`$ nanotubes and those with $`nm`$ multiple of three are metallic, the other ones are semiconducting with an energy gap inversely proportional to their diameter. This was confirmed experimentally by Wildoer and collaborators Wildöer et al. (1998) using a scanning tunneling microscope (STM). They were able to obtain atomically resolved images and to perform local spectroscopy on their samples using the tip of the microscope as local probe.

Electronic transport measurements on individual single–wall carbon nanotubes (SWNT) deposited on metallic electrodes demonstrated their potential as molecular quantum wires Ebbesen et al. (1996); Tans et al. (1997, 1998). Other researchers have done experiments on individually manipulated multiwalled nanotubes (MWNT): Langer and collaborators found effects of localization and universal conductance fluctuations at $`T=20mK`$ Langer et al. (1996), and Frank et al. found quantization of the conductance at room temperature Frank et al. (1998). In the absence of scattering, the momentum relaxation length and the localization length are much larger than the wire length, and the transport is ballistic. This is the case of carbon nanotubes in the absence of lattice defects, where the wavefunction of the electron is extended over the nanotube, and there are only two channels which contributes to the electronic transport giving $`G=2G_0`$ (where $`G_0=2e^2/h`$ is the conductance quantum) Saito et al. (1998); Chico et al. (1996). If we consider elastic scattering, howewer the conductance of the system is described by the Landauer formula, which applies if the wavefunction can spread over the whole sample and is given by

$$G=\frac{2e^2}{h}M𝒯=\frac{2e^2}{h}\underset{\alpha \beta }{\overset{M}{}}|t_{\alpha \beta }|^2,$$

where $`𝒯`$ is the transmission probability for a channel to go from one electrode to another. This probability is given by the sum over the number of channels, $`M`$, of transmission probabilities from $`\alpha `$ to $`\beta `$ channel, $`|t_{\alpha \beta }|^2`$ (for a single–walled nanotube is $`M=2`$). The elastic scattering affects the transmission probabilities and thereby reduces the conductance, which is no longer exactly quantized.

Molecular electronics based on the ensemble of thousands of nanoscale devices cannot be a real alternative to current microelectronics if the devices have to be individually manipulated using local probes. Some groups have tried a new approach using self–assembly of the macromolecules into controlled nets of devices or regular structures that could be later processed using nanolithographic methods, but carbon nanotubes have resisted this kind of treatment because of their difficulty to be diluted in appropriate solvents. Chen and collaborators studied both ionic and covalent solution–phase chemistry with concomitant modulation of the SWNT band structure Chen et al. (1998), and Liu Liu et al. (1998) and Vigolo Vigolo et al. (2000) have prepared stable colloidal suspensions in water with the help of surfactants allowing for collective manipulation and orientation of the nanotubes. Howewer, the interaction between carbon nanotubes and their organic solvents or between carbon nanotubes and polymers in composites or other nanostructured materials is not yet fully understood.

We prepared our samples with purified arc–discharge soot material containing carbon nanotubes, which were first submitted to a sonication bath in 1,2–dichloroethane to unravel the material into nanotube bundles and single nanotubes. After evaporation of the 1,2–dichloroethane the remaining black powder was dissolved in a mixture of poly–chlorinated biphenyls (Aroclor 1254), using sonication and heat treatment to improve the dilution of the carbon material in the viscous solvent. Aroclor 1254 is a solvent of strong temperature–dependent viscosity. Whith our treatment, we arrived at an homogeneous dark solution of carbon nanotubes in Aroclor 1254 which is highly stable: no precipitation has been observed for months at room temperature. Scanning electron microscopy showed that our nanotubes have diameters ranging from 15 to 25nm, and lengths ranging form 0.8 to $`2\mu m`$, with a more or less Gaussian distribution centered at $`1.4\mu m`$. Our estimations of the critical density, $`d_c`$, for different solutions of carbon nanotubes with length $`l=1.5\mu m`$, give the following values: if they were single walled (10,10) nanotubes, $`d_c=1.43\times 10^6g/cm^3`$, if they were multiwalled with ten layers, then $`d_c=1.12\times 10^5g/cm^3`$.

In order to prepare our samples for the transport measurements, we used a solution of density $`c=0.005g/cm^3`$ which is at least $`10^3`$ times bigger than the overlaping concentration $`C^{}`$ of the nanotubes. The presence of overlaping in our samples is important to provide a percolative path between the nanotubes in solution and the gold substrate. We dropped the solution mentioned before on a gold substrate and then spin coated it to reach an homogeneous layer of about $`0.5\mu m`$. We used a scanning probe microscope (SPM) to make contact with the nanotubes and measure the conductance. The platinum tip was mechanically etched and mounted on a rigid STM piezo. We recorded the current as a function of time and tip distance to the sample. The applied bias voltage was kept fixed, so that conductance could be determined as a function of the position of the nanotube contact. We start the experiment with a coarse approach of the tip to the sample. Once electrical contact has been established, the tip is piezo controlled by the SPM and driven in and out cyclically. The vertical travel used is $`2\mu m`$. Data from sequences of 50 to 100 cycles were recorded automatically. When the tip makes contact with the sample, one nanotube can stick to it, remaining stuck while the tip cycles, and thus providing an electrical contact through the nanotube with the gold substrate, directly or through a percolative path with other nanotubes of the solution. While the tip is cycling, the contact of the nanotube with the tip remains constant, otherwise the nanotube contact with the substrate is changing while it is being lifted by the tip. When the tip retracts a distance larger than the length of the nanotube, the contact is broken and the conductance goes to zero. We can record sequences of more than fifty cycles, but usually it is difficult to reach more than seventy cycles without loosing contact with the nanotube.

In figure 1 we plot the conductance $`G`$ versus time for a portion of a sequence in which two steps of different length are observed. This is representative of a high number of steps, which can be used to construct an histogram for every sequence, shown in figure 1B. The histogram is taken over the first fifty cycles. In this figure, a peak centered at the conductance value of $`2G_0`$ can be seen, along with the broadening expected from the noise in our measurement. Note that if the nanotubes behaved as diffusive conductors, we would expect the conductivity to vary as a function of the distance between the nanotube contact points to the substrate and the tip. Howewer, what we observe in our experiment is a series of plateaus, meaning that the conductivity remains constant for periods of time, which in our experiment are equivalent to tip displacements ranging from 0.8 to 1.4$`\mu m`$. This can be explained by changes in the contact point between nanotube and substrate, which sometimes might even involve mechanical deformation (bending) of the nanotube. The value of the conductance remains fixed while the current path is changing, thus the nanotube behaves like a ballistic conductor. Furthermore, the conductance appears to be quantized, taking values only at even integer multiples of the conductance quantum $`G_0`$. This result suggests that the nanotube behaves like a quantum wire, with only two conductance channels. We find this to be true, regardless of the diameter and length of the nanotubes involved. We attribute the different plateau lengths observed to differences in the length of nanotube between tip and substrate, which could be due to the adherence of a different nanotube or to a occasional displacement of the previously stuck one along the platinum tip. The presence of the solvent may enhance nanotube adherence to the tip. Howewer, as the contact heats up, the viscosity of the Aroclor decreases modifying the nanotube–tip structure in an unpredictable way. Nevertheless, the value of the conductance remains fixed at $`2G_0`$. The persistence of this value indicates that only the outer layer contributes to the electronic transport. The inner layers are isolated from the outer because the resistivity in the direction perpendicular to the tube axis is very high, as one could expect from similarities with graphite.

A second kind of results obtained with the same experimental procedure can be seen in figure 2. A series of steps at even integer values of the conductance quantum $`G_0`$ appears when the current along the nanotube is plotted versus the position of the tip. Since in single–walled metallic nanotubes only two channels contribute to electronic transport, we think that the new steps come from transport through other nanotubes. It is possible that two or more nanotubes get stuck to the tip, thus providing more than one connection, also ballistic and quantized, for the current between tip and substrate to flow. When the tip makes contact, sometimes the staircase begins with a step at a value of 4 and occasionally 6 $`G_0`$. When the tip begins to retrieve, the step remains fixed for a while and suddenly changes reducing to the next even integer value of $`G_0`$. We have never observed a step which changes in values higher than $`2G_0`$. This provides support to the idea that only the outer shell of the MWNT contributes, because if two (or more) consecutive metallic shells were contributing, it would provide a step of $`4G_0`$ (or higher) when this nanotube losses contact with the gold substrate. An histogram of twenty cycles is shown in figure 2B, showing peaks at even integer values.

Some of the observed steps present some slight but clear deviations from the quantized values, even after considering the noise of our measurements. This effect is more pronounced at the end of the step, when the current flows through longer distances along the nanotube. We attribute this deviation to the elastic scattering produced by structural defects of the nanotube atomic lattice, which affects the transmission coefficients and deviates the conductance from the quantized value as described by the Landauer formula. Our nanotubes were subjected to ultrasonic bath and this procedure could have produced some structural damage in the nanotubes. Such specific defects may affect the conduction at special energies werheas the transport remains unaffected for different energies.

An alternative explanation of the elastic scattering could be found in the interaction between shells of the multi–walled nanotube. The shell which is contributing to the electronic transport can be affected by its interaction with other shells, introducing an effective potential equivalent to structural disorder. Since we can not resolve the chirality of the successive shells of the MWNT, we could not study this effect. We think that the interaction between shells can be in part responsible for the discrepancy between our results and those found by S. Frank and coworkers (Ref. Frank et al. (1998)). They found a quantized conductance at one unit of the conductance quantum $`G_0`$ instead of the expected $`2G_0`$. They attached a carbon fiber with protruding carbon nanotubes to a gold wire and used it as the tip of a scanning probe microscope. This fiber was found in the soft material inside the hardshelled deposit of an arc-discharge and no aditional treatment was performed to the fiber. We think that different shells of their MWNT can have electrical contact with the gold wire, increasing the intershell interaction even if only the outer shell have contact with the mercury melt. Recent theoretical works Roche et al. (2001a, b) have studied the effect of intershell coupling and incommensurability in MWNT and they found that the mixing of quantum channels results in an enhanced contribution of backscattering and quantum interference effects which gives a quantum correction that reduces the conductivity. This effect could explain the different kinds of electronic transport reported for MWNT: ballistic, diffusive and even insulating behaviours. This may be the case of Ref. Frank et al. (1998), but the exact reduction of one unit of the conductance quantum remains unexplained, unless the mechanism mentioned above eliminates only one of the two channels in the outer shell of the MWNT.

We also measured the current–voltage characteristics for some samples. When an electrical contact was obtained, we varied the bias voltage in a range between 0.1 and 1V. The I/V plot showed a linear behavior, which indicates that the nanotubes are indeed metallic. Howewer, at higher voltages ($`V0.6V`$) we observed some instabilities even if there was no sign of current saturation, those inestabilities may arise from the contribution of different subbands even if a single shell is contributing to transport. For a bias voltage of 1V, the current density is on the order of $`J10^6Acm^2`$, depending on the diameter of the nanotube. The dissipated power is $`1mW`$. The electrical breakdown of the nanotube occurs at higher values of bias voltage, when the non linear regime is reached. It has been reported that for MWNT this breakdown occurs at powers of $`300\mu W`$ Collins et al. (2001). The onset of saturation and the eventual breakdown are linked to a common dissipative process involving the excitation of phonons, temperature raising in the nanotube and a process of oxidation and loss of successive carbon shells. A simple heat transfer analysis, taking the nanotube as a solid cylinder with a bulk thermal conductivity of $`\kappa =3000Wm^1K^1`$Kim et al. (2001) gives for a nanotube of $`1.4\mu m`$ long and 20nm in diameter a temperature $`T=1850K`$, which is not possible, because nanotubes start to oxidate at $`T1000K`$. If we take other reported values such as $`\kappa =25Wm^1K^1`$ Yi et al. (1999), which possibly apply for nanotubes with more defects, or higher diameters, as is our case, the temperature would reach $`T22000K`$. Clearly heat is not dissipated in the nanotube, further supporting that the transport is ballistic.

In summary, we have prepared stable solutions of multiwalled carbon nanotubes in Aroclor 1254. Transport measurements using a scanning probe showed ballistic transport through distances longer than $`1.4\mu m`$ and conductance steps at even integer values of the conductance quantum $`G_0`$.

###### Acknowledgements.

We wish to acknowledge Professor John Schrag for supplying Aroclor 1254, and for helpful discussions and suggestions. This work has been partially supported by Fundación Séneca under Grant No.PI-60-00858-FS-01.


format: text

text:
# Simulations of gaseous disc-embedded planet interaction

## 1 Introduction

Planet formation is believed to occur in discs of gas and dust around young stars. Observations of circumstellar discs in multiple wavelengths reveal the initial conditions for planet formation (McCaughrean, Stapelfeldt & Close, 2000; Wilner & Lay, 2000). Detection of extrasolar planets has provided a diverse set of final states, including a statistically significant number of planets with masses comparable to that of Jupiter on orbits very close to their parent star, so-called ‘Hot Jupiters’ (for a review see Marcy & Butler, 2000, and exoplanets.org for an up-to-date catalogue). Formulating an acceptable theory of giant planet formation in such an environment has proved difficult (Wuchterl, Guillot & Lissauer, 2000, and references therein). Instead planet migration has become the standard theory to explain these interesting companions.

The interaction between a gaseous disc and an embedded protoplanet is a complex, highly non-linear process. Several phenomena have been classified, including the accretion of planetesimals, the accretion of gas, the exchange of angular momentum via tidal forces and planet formation via gravitational instability. How these behaviours combine is determined by the initial conditions, for example the mass and temperature of the disc and the location of an embedded planet. We can sort all possible interactions into a few categories: linear migration via tidal interaction, gap formation and viscous migration, and planet formation via gravitational instability. All these interactions have previously been investigated separately. Here we examine how the interactions are related and how the initial conditions determine which behaviour will be exhibited.

The dynamics of low-mass planets in gaseous discs have been explored both theoretically and numerically. The exchange of orbital angular momentum between disc and embedded satellite via spiral density waves excited at Lindblad resonances was discussed in Goldreich & Tremaine (1980). Ward (1997) gave arguments suggesting that this will almost always lead to inward migration of the embedded planet, and derived an expression for the rate of this migration in the case where a planet does not open a gap in the disc (Type I migration). More recently, Tanaka, Takeuchi & Ward (2002) gave an updated formula for the total torque on a planet, taking into account corotation resonances and mentioning three-dimensional (3D) effects. Previous numerical simulations have tested some of the theoretical assumptions, usually finding general agreement (e.g. Nelson et al., 2000; D’Angelo et al., 2002).

In Type II migration the planet exerts a tidal torque greater than the disc viscous force, opening a gap. The migration time-scale is slowed, and accretion on to the planet is expected to slow or stop. Criteria for gap opening have been derived (Takeuchi, Miyama & Lin, 1996; Rafikov, 2002) and examined numerically using two-dimensional (2D) grids (Bryden et al., 1999; D’Angelo et al., 2002). The accretion process has been simulated (Kley, 1999; Lubow, Seibert & Artymowicz, 1999) finding mass-doubling times for Jupiter-mass planets in light discs to be of the order of $`10^5`$ yr. For a review of Type I and II migration, see Thommes & Lissauer (2002).

The standard picture of migration is situated in a calm, stable disc. What if a single planet forms in a marginally stable disc? First hinted at by Armitage & Hansen (1999), the spiral density waves excited by a compact object (e.g. a rocky core formed via accretion of planetesimals) can lead to clump formation that would not otherwise occur. The fact that an embedded planet affects the disc density and flow locally has been known for some time; here we find that it can also have far-reaching (spatially and temporally) global effects.

Regardless of the embedded planets, some discs are themselves gravitationally unstable. This scenario of planet formation has been explored by Boss (1997, 2001) and recently by Mayer et al. (2002). In unstable discs many planets of Jupiter-mass or greater can be formed very quickly (hundreds of years). The newly formed proto-giant planets will violently tear through the disc, and their interaction as point masses will dominate the subsequent evolution of the system. This scenario is attractive because it easily produces planets with properties similar to observed extrasolar planets. Furthermore, because the time-scale of formation is so small, any system even marginally unstable is likely to form planets during the much longer disc lifetime.

Much of the previous numerical work in this area has been limited to two-dimensional approximations, using grids with a relatively small spatial range, and ignoring the self-gravity of the disc. Some three-dimensional simulations have been performed recently (Bate et al., 2003; D’Angelo et al., 2003), although these have still ignored self-gravity. The smoothed particle hydrodynamics (SPH) formalism was designed to be spatially adaptive. Our implementation uses variable time-steps, giving us temporal adaptivity. This allows us to explore the interaction with a much greater dynamic range. Our simulations are fully three-dimensional, have free boundary conditions, include the self-gravity of the disc and allow the planet to migrate and accrete freely.

Armed with these useful tools, we have simulated planets embedded in gaseous discs, for initial planet masses ranging from 10 Earth-masses to 2 Jupiter-masses, and disc masses ranging from 0.01 to 0.25 $`\mathrm{M}_{}`$. Previous investigations have focused on just one of the interactions described above, usually looking at lower-mass discs. Here we attempt to relate the interactions, and see how the choice of initial conditions determines the classification of the outcome. We draw an outline of interesting regions of parameter space, and present a few early results from each distinguishable region. We also discuss where our simulations differ from previous efforts and suggest that this is due to the more realistic system we are able to model. In future papers we will quantitatively address the specific phenomena we can simulate, such as gap opening criteria, gas accretion through the gap, Type I and II migration time-scales and the likelihood of triggered planet formation.

The plan of the paper is as follows. In Section 2 we describe the initial conditions of our models. The numerical technique is briefly discussed in Section 3 and our main findings are detailed in Section 4. Details of the results are discussed in Section 5 and we conclude with Section 6.

## 2 Initial Conditions

We use power-law distributions for density and temperature profiles in our simulations, and the disc stability criterion and mass ratio are dimensionless, so our results should scale well. None the less, it is convenient to use conventional units when discussing our results. The central star was of 1 $`\mathrm{M}_{}`$, and its potential was softened with a cubic spline on a length-scale of 0.4 AU. All simulations were performed in heliocentric coordinates, taking into account the back-reaction of the disc and planet on the star. Initially, the disc extends from 1 to 25 AU and is free to expand.

#### Temperature profile

Observations of protoplanetary discs do not put tight constraints on the temperature profile. Therefore, we adopt a power law for the radial temperature profile and choose an index which is equivalent to a constant aspect ratio $`H/r`$ (where $`H`$ is the scaleheight of the disc at a particular radius). Most of our simulations used $`H/r=0.05`$ (which gives $`T(5.2\mathrm{AU})=102`$ K). This profile is similar to that suggested by the solar nebula model of Hayashi, Nakazawa & Nakagawa (1985) and is common in other simulations of these phenomena. The temperature is capped at 5000 K within 0.1 AU of the star and limited below by 17 K outside of 30 AU.

#### Density profile

The minimum-mass solar nebula model suggests a power law for the vertically averaged surface density of $`\mathrm{\Sigma }(r)r^{3/2}`$. The vertical structure of the disc is determined by the balance between the vertical component of the gravity of the central star and the vertical pressure support of the disc. If we assume no vertical variation in the temperature, the vertical density structure can be solved, yielding a scaleheight expressed as a function of the radius and the temperature at that radius. An expression for the density distribution is then

$$\rho (r,\theta ,z)=\mathrm{\Sigma }(r)\rho _{\mathrm{vert}}(r,z)=\left(\mathrm{\Sigma }_0r^{3/2}\right)\left(\rho _0e^{z^2/H(r)^2}\right)$$

In a gravitationally stable disc we found this profile to be stable over many dynamical times, so we chose this initial distribution for all of our simulations. Our standard disc, 0.1 $`\mathrm{M}_{}`$ extending from 1 to 25 AU, thus has $`\mathrm{\Sigma }(5.2\mathrm{AU})=150`$ g cm<sup>-2</sup>. For the gas particles of the disc, the initial positions are randomly drawn from the desired distribution. Initial velocities are chosen so that particles are on circular Keplerian orbits, modified slightly to account for the self-gravity of the disc and gas forces. In stable discs without embedded planets, the surface density does not change markedly with time except at the inner and outer edges. The outer edge expands slightly, due to the initial sharp drop in density. At the inner edge, particles migrate inward due to viscous forces, participating in a complex balance between the softened potential of the star, the smoothed gas pressure forces and Keplerian shear. These particles are much closer to the central star than the inner Lindblad resonance of the planet, so we do not believe they play a large role in determining the evolution of the planet. We have not investigated schemes for halting planet migration, which might involve removing the inner reaches of the disc.

#### The embedded planet

The initial mass of the planet particle was a parameter. We used planet masses ranging from 10 $`\mathrm{M}_{}`$ to 2 $`\mathrm{M}_\mathrm{J}`$. The planet particle interacts via gravity only, and no gas drag is applied. Planets significantly more massive than a gas particle can capture gas particles, adding to the effective mass of the planet. These captured gas particles are not removed or treated specially; they remain as part of the envelope of the developing planet. The potential of the planet was spline-softened on a length-scale of one-fifth the initial Roche radius.

### 2.1  An initial gap

Sufficiently large planets can open a gap in the disc by exerting a tidal torque stronger than the local viscous angular momentum transfer. The balance between these forces determines the shape and width of the gap formed around such a planet, which can be calculated analytically by averaging over many orbits and assuming how and where the density waves are dampened. The formation of these gaps has been studied numerically before, primarily in light discs with heavy planets without allowing migration (Bryden et al., 1999), finding reasonable agreement with analytic models. For low-mass discs, the Type I migration time-scale may be sufficiently long for the standard explanation of gap formation to work. We found that planet migration and gap formation could not be decoupled. Therefore, we chose not to impose any initial gap, instead the formation of the gap is simulated in systems that produce one.

A Jupiter-mass planet embedded in an unperturbed disc without a gap is somewhat unphysical: how did the planet become so massive without affecting the disc? As seen in Fig. 3 below, the transition from Type I to Type II migration is very sudden, taking only a few orbital times. Therefore, the sudden appearance of a large planet in a calm disc is not entirely unreasonable, and we can imagine the following scenario: the planet forms, at a lower mass, further out in the disc, and undergoes Type I migration and accretion to obtain the location and mass we start with. To accurately model the capture of gas, the planet particle must be several times more massive than the gas particles. Therefore, future investigations using higher resolution should be able to push the initial planet mass lower.

### 2.2  Equation of state

We use an isothermal equation of state with a spatially fixed temperature profile. In doing so we are assuming that the disc temperature is controlled solely by the distance from the parent star, and that the energy of viscous heating is radiated away much faster than the dynamical time-scale of the gas. We do not account for shielding, so we miss the detailed thermal structure of the clumps that form. Furthermore, clumps on elliptical orbits will heat up and cool down along their orbits. Boss (2002) has found that using a more realistic equation of state suppresses clump formation in marginally stable discs.

### 2.3  Disc stability

The stability of the disc is characterized by the Toomre $`Q`$ criterion, defined as

$$Q=\frac{c_s\kappa }{\pi G\mathrm{\Sigma }}$$

where $`c_s`$ is the sound speed, $`\kappa `$ is the epicyclic frequency, and $`\mathrm{\Sigma }`$ is the surface density of the disc. This criterion measures the susceptibility to a local instability. Global spiral waves can cause the local density to increase, leading to a local collapse. The power laws we have chosen for the temperature and density make $`Qr^{1/2}`$. Hence we expect instabilities, if produced, to occur in the outer regions of the disc.

## 3 Numerical Realization

To model the disc–planet system we used Gasoline, a tree-based particle code originally developed for cosmological simulations (Stadel, 2001; Wadsley et al., 2004). It implements smoothed particle hydrodynamics, a Lagrangian treatment of hydrodynamics for gas. It has been used previously on solar system-scale problems by Richardson et al. (2000) and Mayer et al. (2002). The code was modified to apply a fixed temperature profile to the particles. The equations of motion include the gravity and hydrodynamics of the gas disc. We do not include effects due to magnetic fields.

### 3.1  Boundary conditions

SPH is a Lagrangian technique, so we do not need to impose any boundaries. All the particles are allowed to move where the forces dictate. Particles are not deleted or merged when they approach the central star or the embedded planet.

### 3.2  Simulation parameters

The tree–cell interaction opening angle we used was $`\theta =0.7`$. Using a multi-stepping technique, the particle timesteps are picked according to the local density (with $`\eta =0.1`$), and to satisfy the Courant condition (with $`\eta _c=0.4`$). The largest possible timestep was set to be 2.5 yr, which is 1/50 the orbital time at the outer edge of the disc. The mean molecular weight of the gas is 2 AMU, corresponding to molecular hydrogen.

### 3.3  Viscosity

The SPH formalism treats viscosity differently from the Shakura & Sunyaev $`\alpha `$ prescription popular in finite difference schemes. No ‘real’ viscosity is added to the equations of motion. Because so little is known concerning the sources of real viscosity, we choose to err on the side of caution by introducing as little viscosity as possible into the system, in an attempt to model the fluid with as high a Reynolds number as possible. However, some artificial viscosity is necessary to stabilize the integration of the equations of motion. Gasoline includes the standard Monaghan treatment of artificial viscosity (Monaghan, 1992), including a Balsara switch that reduces viscosity in regions of strong shear (Balsara, 1995). Using the SPH formulation of the equations of motion, we can estimate the kinematic viscosity $`\nu =\overline{\alpha }c_sh/8`$ (Nelson et al., 1998) associated with the artificial viscosity. This estimate is reasonable when there is little or no shocking (the presence of shocking adds artificial viscosity locally). There is strong shear in the disc, so the Balsara switch reduces this by a factor of approximately $`5`$, yielding a viscous time-scale $`\tau =r^2/\nu 10^4`$ yr everywhere for a run with $`10^5`$ particles in a stable disc. Since this is at least 5–10 times longer than our simulations, we believe that artificial viscosity will not seriously affect the dynamics of the disc. To confirm this, we performed two additional simulations where we halved and doubled the standard value of $`\overline{\alpha }=1`$. The change in migration and accretion rates was approximately 5 per cent, which is less than the difference due to the random seed used to determine the initial particle locations. This indicates that artificial viscosity only affects the simulation as a numerical stabilizer, and does not introduce spurious physical effects.

## 4 Results

Looking first for planet migration, we start with planets embedded in light, stable discs. Increasing the disc and planet masses increases the rate of Type I migration, and brings us to gap formation and Type II migration. Continuing to marginally stable discs, we find that the density waves excited by an already formed planet can trigger gravitational instability in discs that would otherwise remain stable. Finally, we push the disc mass high enough so that no external trigger is needed to drive the gravitational collapse mode of planet formation.

Some of the factors that determine which category a system will end up in include the disc density profile and total mass, disc temperature, initial planet mass and orbital radius. We made choices for the density and temperature profiles, and varied the total disc mass, the initial planet mass and the orbital radius. We ran a simulation for each variation of the parameters for a given amount of time (326 yr), and classified the final result. We found the disc stability to be the most influential factor in determining the behaviour of a disc–planet system. The qualitative results of our suite of simulations can be seen in Fig. 1. We have arranged our results in a map of parameter space, illustrating the relative positions of the four interactions.

In stable discs (the upper region), planets of all masses exhibit Type I migration. In more massive discs, planets can open a gap, with a bias toward more massive planets. A precise measurement of the interface between these two regions can show us the criteria for gap opening. To differentiate between the regions labelled ‘Type I’ and ‘Gap Formation’ we look for the halting of the migration of the embedded planet within the time we simulated. See Sec. 5.4 for our definition of a gap and the reasoning behind it.

Some planets can trigger the formation of more planets via gravitational instability. This requires a marginally stable disc and a massive planet. The precise shape we have given this region is a guess based on our experience. In the bottom region, the disc is gravitationally unstable, even in the absence of embedded planets.

Note that the location of the boundaries between regions in this map are dependent on the time period of interest. We are deliberately simulating for only a few hundred years, to match the time-scale of triggered and wholly unstable planet formation. As such, many of the systems we classify as Type I do satisfy the standard thermal and viscous criteria for gap formation, and would form deeper depressions in the density profile and eventually shift from Type I to Type II migration if we evolved them several times longer. In addition, classification over such a short period of time is important, as we measure very high migration and accretion rates.

We have shown a two-dimensional slice through the three-dimensional parameter space we explored. While convenient and, we believe, largely accurate, this does hide some of the complexity we found. For example, we looked for triggered instability by placing planets at different orbital radii in different mass discs such that the value of $`Q`$ was the same at the location of each planet. Since $`Q`$ is dimensionless and we used power laws for the density and temperature, we might expect the same result in each system. However, the system with closer planet and more massive disc did not fragment, while the further planet and less massive disc system did. This may be due to edge effects or a non-linear dependence on the effectiveness of $`Q`$ as a predictor of collapse.

### 4.1  Type I migration

Linear theory suggests that a seed planet in a light gas disc will undergo Type I migration with a rate that scales linearly with the planet mass and disc surface density. The planet excites density waves in the disc which torque back on the planet, leading to a change in the orbital angular momentum of the planet. As a result of the power-law index we chose for the density profile, we always see inward migration. During Type I migration the planet accretes gas from the disc with the help of the spiral density waves it excites. The internal waves take angular momentum from the disc, and the external waves give angular momentum back to the disc. On both sides, mass on nearby orbits is shepherded on to the planet.

Table 1 lists migration and accretion rates for the systems we classified as Type I migration. For the lighter discs, the migration rate is a constant function of time, whereas a more massive disc causes the migration rate to decrease as the simulation progresses and a gap is formed. We find that the accretion rate scales linearly with the planet mass and disc surface density. Further, the migration rate scales linearly with the disc surface density. However, changing the planet mass has almost no effect on the migration rate, confirming the findings of Nelson & Benz (2003). If we compare our results with the analytic expression (eq. 70 from Tanaka et al. (2002)), we find our migration time-scales ($`\tau _M=a/|\dot{a}|`$) to be approximately 6 times shorter. Most previous simulations of Type I migration have kept the planet on a fixed circular orbit, and estimated migration time-scales by measuring the torque on the planet from the gaseous disc. Our migration timescales are actual measurements of the rate of change of the orbit of the planet.

### 4.2  Gap formation

More massive seed planets exert strong enough tidal torques to open a gap in the disc. By clearing the gas from the location of the Lindblad resonances, the torque exerted on the planet by the disc vanishes. Previous theoretical work and simulations have suggested that formation of a gap drastically reduces the accretion rate on to the planet. We have not found this to be true. Even when the gap has cleared, the planet still perturbs the disc. These perturbations act to pull mass off the edges of the gap and on to the planet.

We found that gap formation was dependent more on the disc mass than the planet mass. This agrees with our claim that the self-gravity of the disc is important. At the start of the simulation the planet rapidly collects all the gas particles nearby, forming a gap within approximately 100 years. This accretion of gas at corotation is much faster than the accretion from close spiral arms or from gap edges, as described above. Since the migration time-scale is so short in the massive discs, there is significant migration during the gap formation. After the gap has formed, the migration stops or slows considerably. In 0.1 $`\mathrm{M}_{}`$ discs, the accretion rate slows to a steady value of approximately $`10^5`$ $`\mathrm{M}_{}`$/yr, independent of the mass of the planet. Fig. 2 shows the mass of a planet and its accreted envelope as a function of time, for 0.5, 1, and 2 Jupiter-mass seed planets. By comparing with Fig. 3, the semi-major axis of the planet as a function of time, we see the accretion slow as the planet clears its gap (changes from Type I to Type II migration). The planet peels gas off the edges of the gap, as it catches up with particles on the outer edge and is caught by particles on the inner edge. Our results confirm that the formation of ‘Super-Jupiters’ is possible in massive discs. As Type I migration ceases, the planet is left in the gap with a non-zero eccentricity. We saw eccentricities of 0.01 to 0.04, comparable to the current eccentricity of Jupiter.

The time-scale for Type II migration is the viscous time-scale of the disc. Since we are actively trying to reduce the viscosity of the disc in our simulation, we do not see significant migration after the gap has formed.

### 4.3  Planet–triggered instability

If we reduce the stability of the disc, the spiral arms generated by an embedded protoplanet can trigger the gravitational instability. Such triggered clumps are always exterior to the seed planet that drives the instability, because $`Q`$ decreases with distance from the star. The spiral arms of the planet must present a large enough density contrast to drive the local value of $`Q`$ below 1. Since the unperturbed disc stability decreases with distance from the central star, for a fixed disc stability, a planet farther out from the central star is more likely to trigger the collapse. If we place a Jupiter-mass planet at 12.5 AU, a disc with $`Q1.6`$ will fragment (see Fig. 4). Once the instability is triggered, the evolution of the system is very similar to the totally unstable case. The clumps form with non-negligible eccentricity, and strongly disturb the disc profile.

### 4.4  Unstable discs

If the disc is massive and/or cool enough, gravitational instabilities will cause it to collapse into bound objects, regardless of any embedded planets. Spiral waves propagate, and the density buildup in the arms causes the disc to fragment into gravitationally bound objects. Not one, but several of these bound objects are formed, and their interaction as effective point masses dominates the subsequent evolution. The previously smooth disc is ripped apart, as the clumps often have sizable eccentricities (see Fig. 5 for a progression of snapshots). We observe collisions, ejections, tidal stripping during close encounters, and violent disruption when plunging too close to the parent star. The eccentricities of these clumps can be high, and change as a result of dynamical interaction. As in Mayer et al. (2002), we find that the disc is unstable when Toomre’s $`Q`$ criterion is less than approximately $`1.4`$. The clumps form when a spiral arm with sufficient amplitude reaches the outer regions of the disc.

For this regime we present a simulation that was unstable, a 0.2 $`\mathrm{M}_{}`$ disc which has $`Q_{\mathrm{min}}=0.8`$. In this run, a spiral arm at approximately 20 AU fragments into four separate clumps after 153 years. Five more clumps form over the next 200 years. At this point, two clumps have developed strong leading arms, which fragment into approximately 10 more clumps over the next 30 years. Figs. 6 and 7 show the orbital and mass progression of some of the first clumps formed (the simulation forms many more clumps than detailed in the plots). We see large eccentricities and eccentricity changes, and mergers of several clumps. Once a significant fraction of the gas has been caught up in clumps, we expect that the furious production of new bound objects will cease, and the system will evolve as a set of point masses. After a long time (much longer than we have simulated) the system will reach dynamical stability, presumably with only a handful of objects remaining.

When a spiral arm fragments, it collapses to a rotating spheroid. All our planets formed initially with prograde rotation, however a later merger resulted in a counter-rotating clump. We are not modeling the cooling of the atmosphere of the planet accurately, so cannot comment on what the final spin of the planets will be. Because our simulations are three-dimensional, we can and do observe tilted planets, a result of interactions with other planets and vertical asymmetry in the spiral density waves.

Our work on wholly unstable discs extends the work of Mayer et al. (2002) to a larger region of parameter space, following the system for a longer period of time and spatially fixing the temperature profile. The critical value of $`Q`$ that heralds collapse is quite sensitive to the equation of state used (Boss, 2002). A more detailed equation of state and more physical treatment of radiative transfer are necessary to determine how likely planet formation via gravitational instability is in a given disc. Note that, though we have evolved the system well past the initial collapse, we do not believe we have accurately modeled the detailed structure of the planets that form because we are using a simple equation of state. However, we have accurately followed the gravitational interaction of the multiple planets that form. We can therefore assert that, in planetary systems that formed via gravitational collapse, the disc torque theory of migration does not play a role, as the two-body interactions completely dominate the orbital evolution of the planets.

## 5 Discussion

Here we describe how we calculated the mass of the planet, some tests we performed, further discussion of our migration results, and how we define and use the word ‘gap’ in this work.

### 5.1  Mass of the planet

The gas particles that are captured by the planet are not removed from the simulation or added to the planet particle. They continue to orbit in the disc along with the planet. When we give results for accretion rates we have to calculate a mass of the whole planet, including the initial point mass and the captured gas. We use a simple criterion to determine whether or not a gas particle counts toward the mass of the planet: a particle must be within the Hill sphere of the planet and have density greater than the Hill density (the planet mass divided by the volume of the Hill sphere). Since the Hill sphere depends on the planet mass, this process is iterated until the planet no longer grows, starting from the initial planet particle.

### 5.2  Resolution issues

Our standard run for Type I migration is a 0.01 $`\mathrm{M}_{}`$ disc with a Jupiter-mass planet starting on a circular orbit at 5.2 AU. We simulate the system for 326 years (approximately 28 orbits) and calculate accretion and migration rates as a function of time. Over this time period the migration and accretion rates are constant, and the planet migrates inward approximately 0.3 AU and grows in mass to approximately 1.8 $`\mathrm{M}_\mathrm{J}`$. We ran this model at several resolutions to check for convergence. At resolutions of 20000 to 400000 particles the accretion rates matched to within 5 per cent and the migration rates to better than 1 per cent. Given this agreement, we chose $`10^5`$ particles as our standard resolution.

The results of two $`10^5`$ particle runs can be seen in Fig. 8. Each simulation initially has $`\mathrm{\Sigma }(r)r^{3/2}`$, aspect ratio $`H/r=0.05`$, and an embedded Jupiter-mass planet initially on a circular orbit at 5.2 AU. The two density profiles shown are for disc masses of 0.1 and 0.01 $`\mathrm{M}_{}`$, after 326 years. The disc is three-dimensional, has free boundary conditions, and self-gravity is included. The density profiles can be compared to the results of previous simulations, for example fig. 3 of Nelson & Benz (2003) and fig. 10 of Kley (1999). Those simulations are both two-dimensional, treat viscosity differently, and apply boundary conditions at the inner and outer edges of the disc.

### 5.3  Migration rate differences

The Type I migration rates we found were higher than most previous numerical and semi-analytical studies. However they agree fairly well with a recent simulation of Nelson & Benz (2003), which does include the self-gravity of the disc.

We believe our treatment of the vertical dimension and the disc self-gravity is correct, and so the discrepancy should reveal invalid assumptions made in the linear theory of migration. Self-gravity always acts to increase the magnitude of density perturbations, so we would expect that the torques arising from these perturbations would be correspondingly increased, but we have not completed any simulations without self-gravity, so cannot comment further. To determine the effect of including the vertical dimension of the disc, we did several runs where all particles were constrained to the $`z=0`$ plane. For light discs, the 2D case consistently gave accretion rates lower than the 3D case by at least 20 per cent. For more massive discs, the 2D case appears to give a slightly higher accretion rate during gap formation but again a lower rate after the gap has formed (see Fig. 2). The difference in migration rates between 2D and 3D simulations was not clear.

### 5.4  Defining a gap

Determining if and when a planet forms a gap in a disc is a non-trivial problem. By what percentage must the density decrease around a planet to be called a gap, and how far must this density depression reach? These criteria necessarily invoke an arbitrary number. Fig. 8 shows an azimuthally and vertically averaged surface density of two simulations with different disc masses. In both cases the gas density near the planet has been markedly reduced, and these depressions could be called gaps. The planet in the heavier disc has stopped migrating (see Fig. 3). In the lighter disc, however, the planet continues to migrate inward at a rate independent of the depth of the density depression. Therefore, we choose to use the word ‘gap’ only when its dynamical effects are apparent, i.e. when the planet has stopped migrating. This is how we chose between the ‘Type I’ and ‘Gap Formation’ regions when classifying simulations for Fig. 1.

## 6 Conclusion

This work presents the first results from our ongoing investigation of planet formation and migration. It is a pioneering application of global SPH simulations to a field previously examined primarily with Eulerian grid-based simulations. Using the new technique, we have measured migration and accretion time-scales, witnessing Type I migration and gap formation. Our results for this region partially qualitatively agree with the linear theory of interaction with spiral density waves. We have observed planet formation via disc fragmentation, including a cascade effect from a single starter planet. Using massive discs, we find that it is easy to form large (tens of Jupiter masses) planets on a short time-scale (hundreds of years). If planets form via gravitational instability, the theory of migration via disc torques is not applicable.

To obtain an overview of gaseous disc–embedded planet interactions, we performed a suite of simulations, varying the planet mass, the initial location of the planet, and the disc stability. We have shown the orientation in parameter space of four broad categories of interaction. These interactions are: Type I migration; gap formation and Type II migration; planet-triggered instability; and unstable discs. The stability of the disc is the most effective predictor of the type of interaction between a gaseous disc and an embedded planet.

Planet formation and migration is a non-linear process that requires the reasoned application of computational resources to model. Future investigations will explore the details of each interaction we have described. In addition, we will refine our knowledge of the boundaries between regions of parameter space. Future observations will give us a range for $`Q`$ and a better understanding of the temperature profile in protoplanetary discs. Combined with the results presented here, we can say which migration and formation processes will be relevant.


text:
# Elemental Abundances in NGC 3516

## 1 Introduction

X-ray and UV spectra of the Seyfert 1 galaxy NGC 3516 (z=0.008836$`\pm 0.000023`$, Keel 1996) have shown evidence for a significant column of ionized gas along the line-of-sight to the nucleus (Kriss et al. 1996; Netzer et al. 2002; Kraemer et al. 2002). A recent Chandra observation of this source found it to exist at a historical low-state. The spectrum at that epoch, and comparison with previous X-ray data led Netzer et al. (2002) to suggest a model for NGC 3516 where a constant column of gas reacts to changes in the nuclear ionizing flux. That gas was constrained to have density $`>2.4\times 10^6`$cm-3 and exist at a distance $`<6\times 10^{17}h_{75}^2`$ cm from the nucleus. Netzer et al. (2002) also detected a strong Ovii 0.561 keV line and marginally detected a Nvi 0.419 keV line. Line measurements from the Chandra Low Energy Transmission Grating (LETG) data appeared consistent with emission from the X-ray absorber.

We present here XMM-Newton (hereafter XMM) RGS grating spectra of NGC 3516. These data allow us to rethink the origin of the X-ray line emission in NGC 3516.

## 2 The XMM RGS Observation

An XMM observation of the Seyfert 1 galaxy NGC 3516 was performed covering November 09 UT 23:12:51 – 11 UT 10:54:19. This observation was part of a multi-satellite campaign which including overlapping observations by RXTE and Chandra. The Fe K$`\alpha `$ line shows interesting structure and evolution during the observation, and the combined data in the hard X-ray regime are detailed by Turner et al. (2002). Here we present a detailed analysis of the XMM RGS grating data from 2001 November.

As noted by Turner et al. (2002) NGC 3516 had a flux $`F_{210keV}1.31.5\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, during the November 2001 observations. The source was in this flux state during the previous Chandra LETG observation (Netzer et al., 2002) and as observed by ASCA during 1999 (Figure 1). Unfortunately no useful RGS data were obtained from an earlier epoch observation with XMM in 2001 April, due to background flares during a period of high solar activity.

RGS data were processed using the SAS 5.3.3. version of rgsproc and spectra were extracted using standard regions (the source cell encompassing 97% of the cross-dispersed counts) and extraction criteria, resulting in exposure times $`114`$ ks in RGS-1 (R1) and $`108`$ ks in RGS-2 (R2). The full-band (0.34-2.0 keV) RGS background-subtracted count rates were $`0.062\pm 0.0008`$ (R1) and $`0.067\pm 0.0008`$ (R2). The background comprised 17% of the total count rate in R1 and R2, respectively. We note that due to problems with some of the CCD chips onto which the RGS spectrum is dispersed there are two prominent data gaps evident, R1 has a gap between $``$ 0.9-1.2 keV and R2 has a gap between $`0.51`$ \- 0.62 keV. We note that in this analysis we adopted cz=2649 km/s Keel (1996).

Ignoring emission and absorption lines for the moment, we find a good parameterization of the continuum shape in RGS data over the 0.3 - 2 keV band using a powerlaw of photon index $`\mathrm{\Gamma }=2.17_{0.06}^{+0.04}`$ attenuated by a Galactic column of cold gas $`N_H=2.9\times 10^{20}\mathrm{cm}^2`$ plus a column of ionized gas with $`N_{H(UV)}=7.7\times 10^{21}\mathrm{cm}^2`$. Such a column of ionized gas is expected from the UV and LETG analyses (Netzer et al., 2002; Kraemer et al., 2002) of data from an epoch close in time to these RGS data (and the ionized gas is modeled using CLOUDY90 for solar abundance material, Ferland et al. 1998). Inclusion of the ionized gas improves the fit from $`\chi ^2=1594/639dof`$ to $`1180/638dof`$ and provides a good model to the gross shape of the soft spectrum. Figure 2 shows the data compared to this model. However, our parameterization of the continuum shape is not the only possible model. A recent alternative applied to some Seyfert spectra parameterizes broad spectral wiggles as emission lines due to ionized material in the accretion disk (e.g. Branduardi-Raymont et al. 2001). However, we do not pursue such a parameterization in this paper. The absorption edges expected from the (previously observed) ionized absorber provide a good description of most of the broad features, and we believe this makes our interpretation of the spectrum a compelling way to proceed for this observation of NGC 3516.

The most immediate result evident from Figure 2 is that the RGS spectra show several prominent emission lines. The two strongest lines are those at 0.4198 keV and 0.5614 keV. We identify these with forbidden components of Nvi and Ovii, respectively. These two strong lines were measured carefully to determine bulk and turbulent velocities of the emitting gas. For the Nvi line we measured a rest-energy $`E=419.8\pm 0.2`$ eV and for the Ovii line $`E=561.4\pm 0.3`$ eV. Errors are 90% confidence. Intrinsic line widths were found to be $`\sigma =3.8_{3.8}^{+3.5}\times 10^4`$ keV and $`\sigma =5.6_{5.6}^{+4.4}\times 10^4`$ keV, respectively. The energy for Ovii f corresponds to blueshift of $`214_{80}^{+160}`$ km/s while for Nvi the blueshift is $`143_{143}^{+143}`$ km/s, i.e. consistent with no blueshift. For comparison, the absolute wavelength calibration for the RGS is $`8`$m Å, or $`100`$ km/s at 0.5610 keV. Thus the evidence for a blueshift for the Ovii line is marginal, the lower limit is close to the accuracy of the detector absolute wavelength calibration. The line widths correspond to velocities FWHM(Nvi) $`=646_{646}^{+581}`$ km/s and FWHM(Ovii) $`=704_{704}^{+554}`$ km/s. Line widths are consistent with zero, the upper limits constrain velocity broadening to $`1300`$ km/s. (For comparison the FWHM resolution of R1 is $`680`$ km/s at the observed energy of Nvi, 0.416 keV and $`850`$ km/s at the observed energy of Ovii, 0.556 keV). Changes in the detector effective area with energy cause some apparent asymmetries at first glance as the effective area changes significantly across some line profiles (e.g. Oviii at 0.654 keV). However we find no true asymmetry in any emission line profile.

To search for weak emission lines we slid a Gaussian template across the data, testing for an improvement to the fit at every resolution element compared to our model of the underlying continuum. The data were first binned by a factor of 8 to a binsize of 0.32 Å to ensure that there were greater than 20 photons in each spectral bin so $`\chi ^2`$ fitting could be employed. Examination of the RGS background spectra shows a broad bump centered around 0.39 keV (not fully understood but thought likely due to relatively high dark current in CCD2 for both RGS instruments; XMM Users Handbook, also see den Herder et al. 2001). Other known features in the RGS background include emission from Al K$`\alpha `$ close to 1.5 keV, due to Al in the detector housing. This means we need to subtract the RGS background to ensure we isolate the spectral properties of the source. Figure 3 shows the results of this sliding Gaussian test against the observed data in red. Plotted is the improvement in the fit-statistic. This the plot is useful in assessing the significance of each line, but not the strength of the line. It can be seen that the addition of an emission feature results in a large improvement in $`\chi ^2`$-statistic at several energies.

In order to assess the reality of these features, we have constructed a series of simulated datasets. For each RGS, 20 separate background spectra were generated by randomly adding the appropriate Poisson noise to each channel in the observed background spectrum. Similarly, 20 simulated source spectra were generated for each RGS, again with Poisson noise introduced, using our model of the underlying continuum convolved with the instrumental response. These simulated data sets were then combined and 20 R1+R2 pairs analysed in the same manner as the observed spectra above. Specifically, the sliding Gaussian test was performed for each of the 20 pairs, and the maximum reduction in fit-statistic at each energy noted. The result is plotted in blue in Figure 3, and therefore represents a 95% significance threshold at each energy (and takes into account all instrumental and background features).

Including the Nvi and Ovii lines discussed above, we find a total of 13 emission features (some of which are blends) in the observed spectrum above the 95% confidence threshold. Ten of these 13 occur at energies consistent with atomic transitions of abundant ions, and have understandable intensity ratios (below). Thus we consider these 10 feature to be robust detections. They are labeled on Figure 3 and listed in Table 1. Of the remaining three features, that at 0.39 keV is most likely due to the detector background described above. The final two features (at 0.697 and 0.706 keV; and apparent strengths corresponding to $`0.4\times 10^5`$ photons cm-2 s-1) do occur at energies close to atomic transitions. The former, is consistent with the Ovii ($`1s^21s4p`$) resonance line at 0.698 keV. However, the lack of any indication of the Ovii ($`1s^21s3p`$) line at 0.665 keV makes such an interpretation extremely unlikely. The line at 0.706 keV could be identified with Sxiv RRC at 0.707 keV. However, such an interpretation would lead us to expect a much stronger Oviii RRC at 0.871 keV than observed. In summary, only the 10 emission lines which we believe secure ($`\delta \chi ^2>10`$) and attributable to NGC 3516 appear in Table 1.

While some of the detected features have some significant width, e.g. the Cvi RRC, the features were too weak for us to fit for the width and thus the temperature of the emitting gas. However we did obtain an upper limit on width, $`kT0.01`$keV, i.e. $`T120,000`$K. This is consistent with the width of the RRC features in NGC 1068, based on which Kinkhabwala et al. (2002) argued that the X-ray emitting plasma was purely photoionized.

The emission lines were assumed to originate outside of the UV absorber (i.e. not absorbed by it). If the emitter were inside the UV absorber it would have to reside within $`10^{16}`$ cm of the nucleus with a density $`10^9`$ cm-3 or so (Kraemer et al., 2002). While such a location is possible, the velocity widths of the emission lines at FWHM$`<1300`$ km/s are rather narrow for such a small radial location. These widths are more suggestive of an origin on the pc-scale. Furthermore, some results from Seyfert 2 galaxies suggest the X-ray emitting gas in general has extent on the 10-100pc size-scale (e.g. Sako et al. 2000). All things considered, it seems most likely the emitting gas lies outside of the UV absorber in this case. The fact the absorption edges are imprinted on the soft spectrum must then be explained as the soft-band having a significant contribution from nuclear continuum emission, transmitted through the UV absorber (the fraction of the nuclear continuum scattered by the line-emitting gas is negligible, as shown in Section 3.2). Thus we suggest a full model of absorbed continuum component plus unabsorbed lines (with Galactic column absorbing all). The line fluxes in Table 1 are the inferred intrinsic fluxes at the source (i.e. corrected for Galactic absorption). In addition to tabulating the fitted line fluxes, we noted the total counts in the bin, those attributable to the continuum emission and thus those attributable to the emission line itself. Errors are calculated on the counts in the line and the same fractional errors are applicable to the fluxes listed.

The process was repeated to search for weak absorption lines. Following a similar process to that used to detect the emission lines, we found evidence for an absorption feature at 0.358 keV and broad features centered on 0.753 keV and 0.96 keV. The former has no obvious identification. The feature at 0.75 keV, while close to a jump in detector efficiency, appears real, and we identify this with the unresolved transition array due to inner-shell absorption by Feix - xi M-shell ions, as observed in some other AGN (e.g. IRAS 13349$`+`$2438, Sako et al. 2001). The feature at $`0.96`$ keV may be a blend of Fe lines from highly-ionized gas (Fe xx \- xxiv). Much weaker features at 0.54 keV and 0.61 keV do not have clear identifications.

Absorption lines and edges were expected to be detected in NGC 3516 from transmission of the nuclear continuum through ionized gas, previously observed in this source (e.g., Kriss et al. 1996; Netzer et al. 2002; Kraemer et al. 2002). The overall shape of the RGS spectra shows the imprint of the absorption edges expected from the X-ray absorption associated with the gas measured by STIS (Kraemer et al., 2002). The width of the UV absorption lines is $`2070`$ km/s, X-ray absorption lines with the same velocity structure would be too narrow to detect with the RGS (or Chandra grating data).

## 3 FUSE data

We reduced data from an archival Far Ultraviolet Spectroscopic Explorer (FUSE) spectrum, originally obtained on 2000 April 17, again, while NGC 3516 was in an extremely low flux state (see Kraemer et al. 2002). From the FUSE spectrum we measured an Ovi $`\lambda `$ 1038 flux of $`1.6\times 10^{14}`$ ergs s-1 cm-2, with a FWHM of $``$ 250 km s-1. Given the intrinsic 2:1 emission ratio of Ovi $`\lambda `$ 1032/ Ovi $`\lambda `$ 1038, line fluxes the total Ovi flux should be $``$ 4.0 x 10-14 ergs s-1 cm-2. However, as noted by Hutchings et al. (2001) in their analysis of this spectrum, the observed Ovi doublet ratio is $``$ 1:1, which they attributed to line transfer effects in optically thin gas. We think that it is more likely that the lines are absorbed by intervening gas, which may explain the asymmetric line profiles (see Hutchings et al. 2001) and the narrowness of the Ovi lines compared to the X-ray lines. However, it is not necessary that this absorption arises in the same material that produces the deep, variable X-ray absorption discussed in Netzer et al. (2002) and Kraemer et al. (2002), but could be a relatively small column of ionized gas (e.g., $``$ 1020 cm-2, $`U`$ $``$ 0.5) that covers the NLR of NGC 3516. This component may only be detectable when viewed against the NLR emission, much of which was missed in the small 0$`^{\prime \prime }.`$2 $`\times `$ 0$`^{\prime \prime }.`$2 STIS aperture, but is well-covered by the 20$`^{\prime \prime }.`$0 $`\times `$ 20$`^{\prime \prime }.`$0 FUSE aperture. (At a distance of 35.7 Mpc for NGC 3516, $`1^{\prime \prime }=173`$pc.) In fact, similar evidence for absorbed UV narrow-line profiles has been seen in STIS spectra of NGC 1068 (Kraemer & Crenshaw, 2000).

## 4 Physical Conditions in the Emission-Line Gas

### 4.1 Emission-Line Diagnostics

The dominance of X-ray emission lines in the soft spectrum of this Seyfert 1 galaxy is consistent with the low flux state of the source. If these lines have a $``$constant flux over years, then they would have been swamped by the nuclear emission at higher flux states.

This source is one of the few observed to date with detectable forbidden emission lines from 3 abundant elements in the He-like isoelectronic sequence (Nvi, Ovii and Neix). It is the first such case for a Seyfert 1 galaxy. Even at first glance, the emission line ratios look unusual. Compare for example the XMM RGS spectrum of NGC 3783 where the Nvi line was very weak compared to Ovii, as expected for solar abundance gas (Blustin et al., 2002). The details of these anomalous line ratios are discussed, in the context of the photoionization models, in the following section.

As suggested by various authors (Pradhan et al. 1985; Liedahl 1999; Mewe 1999; Porquet & Dubau 2000), the physical conditions of photoionized plasmas can be constrained via the relative strengths of the most intense lines from He-like species; the resonance line (1s$`^2`$1S0 – 1s2p 1P1), the two intercombination lines (1s2 1S0 – 1s2p 3P2,1), and the forbidden line (1s2 1S0 – 1s2s 3S1). The ratios of the forbidden line to the two intercombination lines (the $`R`$ ratio) is sensitive to density, while the ratio of the sum of the forbidden and intercombination lines to the resonance line (the $`G`$ ratio) is sensitive to electron temperature. If the resonance lines are relatively strong, the latter ratio can indicate temperatures in excess of those expected for a photoionized plasma. However, the resonance lines can be also be photo-excited (Sako et al., 2000), hence the $`G`$-ratio may not be a reliable indicator of either the temperature or means of excitation of the plasma. Based on their analysis of an XMM-Newton RGS spectrum of NGC 1068, Kinkhabwala et al. (2002) suggested that the excess emission in virtually all the detected resonance lines was the result of photo-excitation in a purely photoionized plasma.

In these data, the $`R`$ ratio of the Ovii f-line to the i-lines is $``$ 3.5 (see Table 1), which for a photoionized plasma implies an electron density ne $``$ several x 109 cm-3 (Porquet & Dubau, 2000). However, the weakness of the intercombination lines, plus any underlying contribution from the wings of the adjacent forbidden and resonance lines, can easily lead to an overestimate of the flux. Hence, it is likely that ne is lower. Furthermore, at sufficiently high density, the metastable 1s2s 3S1 level begins to collisionally depopulate, suppressing the forbidden line of the He-like triplet. For Nvi this occurs at ne $`>`$ 108 cm-3, while the depopulation of corresponding Ovii level occurs at densities $``$ factor of 10 greater (Porquet & Dubau, 2000). Given the relative strengths of the Nvi and Ovii forbidden lines, we suggest that the average density in the emission-line gas is ne $``$ 108 cm-3.

### 4.2 Photoionization Models

Photoionization models of the X-ray emission-line gas were generated using the Beta 5 version of CLOUDY (Gary Ferland, private communication). This version of the code includes new atomic data for the He-like ions. We modeled the emission-line gas as a single-zoned slab of atomic gas, irradiated by the central source, for which we used the same spectral energy distribution and luminosity ($`10^{43}`$ erg s-1) described in Kraemer et al. (2002). As usual, the models are parameterized in terms of the ionization parameter $`U`$, which is the ratio of ionizing photons per H atom at the ionized face of the slab. We initially assumed solar elemental abundances (Grevesse & Anders, 1989), which are, by number relative to H, as follows: He $`=`$ 0.1, C $`=`$ 3.4 x 10-4, N $`=`$ 1.2 x 10-4, O $`=`$ 6.8 x 10-4, Ne $`=`$ 1.1 x 10-4, Mg $`=`$ 3.3 x 10-5, Si $`=`$ 3.1 x 10-5, S $`=`$ 1.5 x 10-5, and Fe $`=`$ 4.0 x 10-5. The gas was assumed to be free of dust.

Prior to generating new models, we explored the possibility that the X-ray emission lines arise in the UV/X-ray absorber described in Kraemer et al. (2002). In the current low-flux state of NGC 3516, the models of the strongest UV components predicted ratios of ionic columns for Ovi:Ovii:Oviii to be roughly 1.00:0.67:0.03. Consequently, little Oviii Ly$`\alpha `$ emission is expected relative to Ovii emission for the UV gas seen in emission. This is not the case, we observe enough Oviii emission to know there must be an additional component of higher-ionization gas present. We also note the model for the production of Ovi $`\lambda \lambda `$ 1032,1038 emission lines from the UV gas yields fluxes $``$ 20 times as strong as the Ovii f-line. In the data they appear weaker than this, as noted further below.

Since the X-ray emission-lines cannot be formed in the UV/X-ray absorbers, we permitted $`U`$ and the total column density (NH $`=`$ Hi $`+`$ Hii) to vary, as a new model is required for the X-ray emitting gas. In Table 2, we compare the model predictions to the measured emission-line fluxes by scaling the predicted Ovii f-line flux to that observed. We find that line ratios can be roughly matched using a single-zoned model, with $`U`$ $`=`$ 1.1, NH $`=`$ 1.5 x 1021 cm-2, and solar abundances (Model 1). In order to boost the strength of the resonance lines via photo-excitation, we introduced a turbulent velocity of 50 km s-1, which is roughly equal to the average FHWM of the individual UV absorbers (Kraemer et al., 2002). If we assume ne $`=`$ 108 cm-3, our model results place the emission-line gas at a radial distance of 0.02 pc; although this is within the range determined for the UV absorbers, since the lower limit to the density is not well constrained. It is quite plausible that the emission-line gas lies further from the nucleus, as we discussed in Section 2. The model predicts a mean electron temperature of 7.4 x104K, in agreement with the constraints determined from the Cvi RRC.

The model predictions assuming solar abundance are within the measurement errors for the Ovii He-like lines, Neix f, and Oviii Ly$`\alpha `$, while the prediction for Cvi Ly$`\alpha `$ is slightly high. However, the strengths of the Nvi He-like lines and Nvii Ly$`\alpha `$ line are unpredicted by factors of $``$ 4 and $``$ 2, respectively, which is strong evidence against our initial assumption of solar abundances. The apparently anomalous line ratios in NGC 3516 suggests one of two possibilities; the other elements such as C and O are depleted onto dust grains (e.g. Shields & Kennicutt 1995), or the N abundance exceeds the Solar value.

Considering the case of depletion onto dust grains. UV data show many cases of enhanced nitrogen relative to carbon, e.g. NGC 5548 (Crenshaw et al., 2003) and Akn 564 (Crenshaw et al., 2001). In those cases C may be depleted onto dust grains. However, in NGC 3516 the data argue against dust depletion. While the ratio of N/C is high, the ratio N/O is also high and yet O is not heavily depleted onto dust grains (see Snow & Witt 1996). Also the ratio N/Ne is high and Ne cannot be depleted onto dust grains. Thus we rule out depletion.

Hence, it is most likely that nitrogen is over-abundant, perhaps by a factor of 2 – 3. Interestingly, Kinkhabwala et al. 2002 came to a similar conclusion for NGC 1068. As a result, we regenerated the model using a N/H ratio of 2.5 x solar (Model 2). As shown in Table 2, the fit for the Nvii line is quite good. While the Nvi lines are much closer to the observed strengths, they are still somewhat underpredicted. The latter may simply be a reflection of the difficulty in matching the emission lines with a single-zone model, since some additional Nvi could arise in a component of lower ionization gas. Our discovery of nitrogen overabundance led us to revisit the absorber used which fits the UV data and models the broad features in the soft X-ray band. First we note that the STIS data on which that absorber was based (Kraemer et al., 2002) showed saturated N and C lines, so limits from those lines did not provide tight constraints on the state of the gas for the UV model. Fine-tuning abundances for the absorber model does not make any significant difference to the fit to broad features and thus to the emission line analysis and conclusions drawn here.

Based on the model predictions, we can constrain the global covering factor of the emission-line gas. Our predicted flux for Ovii f emitted from the illuminated face of the photoionized slab is 1.7 x 105 ergs cm-2 s-1. Adopting a distance of 35.7 Mpc for NGC 3516, assuming Ho $`=`$ 75 km s-1 Mpc -1 (Ferruit, Wilson & Mulchaey, 1998), the total O VII f emission from NGC 3516 is $``$ 5.5 x 1039 ergs s-1. Hence, the emitting surface area must be $``$ 3.2 x 1034 cm2. Based on our assumed density, ne=$`10^8`$ cm-3, and luminosity in ionizing photons, the emission-line gas lies at a radial distance of 6.9 x 1016 cm, which yields a total surface area of 6.0 x 1034 cm2. Hence the covering factor of the emission line gas is $`F_c`$ $``$ 0.5. While we expect that the emission-line gas extends over a range in radial distance, and hence density, the total covering factor will be similar.

From the predicted covering factor and column density, we can determine the fraction of the nuclear continuum reflected via electron scattering within the emission-line gas. Assuming isotropic scattering, at small electron-scattering optical depths $`\tau _e`$ $`<`$ 1, the reflected fraction of continuum radiation $`f_r`$ $``$ Nelectron $`F_c`$ $`\sigma _T`$, where Nelectron is the electron column density ($``$ NH) and $`\sigma _T`$ is the Thomson cross-section. We find $`f_r`$ $``$ 5 x 10-4, hence a negligible fraction of continuum radiation will be scattered into our line-of-sight by this component.

Returning to the question of the strength of Ovi emission. The total Ovi flux derived from the FUSE data is $``$ 4.0 x 10-14 ergs s-1 cm-2. Our X-ray emission-line model predicts a total Ovi flux of 6.75 x 10-14 ergs s-1 cm-2, which, since it is somewhat higher than the FUSE value, may support the indication (from the doublet ratio) that these lines are absorbed. Under the assumption that the strengths of the emission lines are constant this Ovi measurement provides a tight constraint on the ionization state of the NLR gas in which the UV and optical emission-lines arise.

## 5 Discussion

The relation between abundance and AGN redshift, luminosity, and the inter-relation between the AGN and starburst regions are key to understanding AGN formation and evolution. Based on the ratios of the Nv $`\lambda `$ 1240 emission line to Civ $`\lambda `$ 1550 and Heii $`\lambda `$1640, it has been argued that QSOs show high metallicity (Z$``$ 1, where Z $`=`$ 1 for solar abundances) and, in particular, high N abundance at redshifts as high as $`z`$ $``$ 3, indicative of enrichment due to rapid star formation at epochs as early as $``$ 1 Gyr (Hamann & Ferland, 1993). Among relatively nearby ($`z`$ $``$ 0.05) AGN, Narrow-Line Seyfert 1s show evidence of large N abundances (e.g., Wills et al. 1999). It is, therefore, interesting that the X-ray spectra of NGC 3516, a low-redshift, Broad-Line Seyfert 1 galaxy, also shows evidence for an overabundance by a factor of 2 –3 of nitrogen in the central regions of the galaxy, with respect to other heavy elements. Furthermore, Kinkhabwala et al. (2002) presented evidence for a similar enhancement of nitrogen in the prototypical Seyfert 2 galaxy NGC 1068. Also, based on the strength of the N V $`\lambda `$1240 line, there is evidence of anomalous N abundances in the NLRs of the Seyfert 1 galaxies NGC 5548 (Kraemer et al., 1998) and NGC 4151 (Kraemer et al., 2000), although this could be at least partly due to enhanced emission by photo-excitation.

Nitrogen can be synthesized from carbon and oxygen produced within a star, referred to as primary production, or via the CNO process in intermediate mass stars (M $``$ 7 M; see Maeder & Meynet 1989) which already possess C and O, i.e., secondary production (e.g., Tinsley 1980). At certain temperatures N is enhanced at the expense of C (while at higher temperatures it would be at the expense of O). Thus it is possible to observe solar O, Ne etc yet obtain enhanced N with subsolar C. The enrichment of nitrogen via secondary production increases with overall metallicity such that N/H $``$ (O/H)2 $``$ Z2. Since the evidence for overabundant nitrogen is typically the high ratios of N lines to those of C or other heavy elements, one might expect that a factor of 3 enhancement of N/O would require Z $``$ 3 and N/H $`9`$ times solar! However, for Z $``$ 1, it is still possible to have N/H a few times solar, assuming the enhancement of nitrogen is accompanied by a loss of carbon and, possibly oxygen (e.g. Maeder & Meynet 1989). For example, assuming roughly solar initial abundances (Grevesse & Anders, 1989), N/H could be 2.5 times solar if approximately half of the carbon were used up in the production of nitrogen. Interestingly, our models with solar carbon abundance overpredicted the strength of the Cvi L$`\alpha `$ line by a factor of $``$ 1.5. To explore the possibility that some of the carbon was lost in the production of nitrogen, we generated a third model, with the carbon abundance set to 47% solar (to fully account for the enhancement of nitrogen) and all other model parameters kept identical to those of Model 2. The predicted line fluxes for Model 3, scaled to the Ovii f-line, are listed in Table 2. The Cvi L$`\alpha `$ line is now within the measurement errors, while the predictions for the other lines are essentially unchanged from those of Model 2. Hence, the emission-line fluxes are consistent with emission from photo-ionized gas in which the nitrogen is enhanced and the carbon depleted by the same amount. If the anomalous N/O and N/H ratios are consistent with Z $``$ 1 and are the result of nitrogen production in intermediate mass stars, these results provide tight constraints on the history of star formation in the nucleus of NGC 3516.

## 6 Conclusions

We have used XMM-Newton RGS grating spectral data to examine the physical conditions within the X-ray emission-line gas in the Seyfert 1 galaxy NGC 3516. The spectra show emission-lines from the He-like ions of N, O and Ne, and H-like ions of C, N, and O. Also, we have detected RRC from Cvi and Oviii. We have shown the following:

1. The RGS data show the soft X-ray absorber to be consistent with the UV absorbers detected in earlier HST/STIS observations. However, the UV absorbers cannot account for the X-ray line emission. Although the UV absorbers could account for some of the Nvi and Ovii emission, they are generally too low an ionization state to produce any of the higher ionization lines detected in the RGS spectra. Furthermore, if the emission from the UV absorbers were scaled to fit the Ovii f component, the models predict Ovi $`\lambda \lambda `$ 1032,1038 lines almost an order of magnitude stronger than observed in a recent FUSE spectrum. While some line variability may occur, an order of magnitude change in Ovi between the FUSE observation and the epoch reported here is highly unlikely.

2. From the Cvi RRC, we find $`kT`$ $``$ 0.01 keV, consistent with low temperature ($``$ 105 K), photoionized gas. Based on the ratios of the f- and i- lines from He-like Ovii and the relative strength of Nvi f, we suggest that the gas is in the low-density regime (ne $``$ 108 cm-2). We have been able to fit the observed emission-line ratios with a single zoned photoionization model, with NH $`=`$ 1.5 x 1021 cm-2 and $`U`$ = 1.1, however, in the case of solar abundances, the Nvi and Nvii lines are significantly underpredicted. Hence, we suggest that the N/H ratio is at least 2.5 times solar, which may be the result of secondary production of nitrogen in intermediate mass stars. Follow-up observations of the stellar population in the nucleus of NGC 3516 could help test this possibility.

3. In order to produce the observed Ovii line fluxes, the emission-line gas must have a global covering factor of $``$ 0.5. The scaled (predicted) Ovi $`\lambda \lambda `$ 1032, 1038 emission is slightly higher than that seen in the FUSE spectrum, although the Ovi line ratios and profiles show strong evidence for absorption, possibly from an additional component of UV absorption near systemic that cannot be deconvolved from the strong UV components seen in the STIS spectra. Hence, the X-ray model is generally consistent with the FUSE spectrum.

## 7 Acknowledgements

We are grateful to the XMM satellite operation team and to Martin Still for help with the RGS data. T.J. Turner acknowledges support from NASA grant NAG5-7538. We thank Tahir Yaqoob, Fred Hamann, Allen Sweigart, Fred Bruhweiler, Mike Crenshaw, Jose Ruiz and the anonymous referee for useful comments. We especially thank Gary Ferland and Ryan Porter for their generosity in allowing us to use the latest pre-release version of CLOUDY.