text:
# Structural properties of silicon dioxide thin films densified by medium-energy particles

## I Introduction

In response to the rapidly increasing demand for highly-specialized applications in optical filtering devices, integrated electronics and photonics technology, considerable efforts have been devoted in the past few decades to the development of methods for growing thin silicon dioxide (SiO2) coatings having low stress and low defect concentrations, and exhibiting bulk-like optical properties. To this end, methods commonly referred to as ion-assisted deposition (IAD), in which growth is accompanied by a flux of energetic particles, are of particular interest: indeed, bombardment with low to medium-energy neutral or ionic particles — in the range 10–200 eV — is found to have beneficial effects on a number of physical properties of the deposited material. Besides fulfilling the above requirements, it also leads to reduced moisture absorption, higher density, lower chemical etch rate, better adhesion to the substrate, and therefore better overall long-term behavior. In addition, IAD methods allow processing at lower substrate temperatures, making them especially adapted to growth on polymers or graded-index coatings. While there evidently exists a threshold energy below which benefits are insignificant, too high an energy can have undesirable consequences, such as preferential sputtering or impurity incorporation.

In ion-beam assisted deposition (IBAD) — a special case of IAD — an external source provides a well collimated flux of ions with a narrow kinetic energy distribution. It has proven useful to discuss the dependance of structural changes of IBAD films upon bombardment in terms of three parameters: the mean kinetic energy of the ions, $`E_i`$, the ratio of ion-to-condensing particle fluxes at the surface, $`\varphi _i/\varphi _n`$ (where $`\varphi _n`$ is measured from the growth rate), and the energy delivered by each particle that condenses, $`E_p=E_i\varphi _i/\varphi _n`$. It has been observed that there exist “critical” values of these parameters below which the films are porous and above which they are densely packed. In the latter case, films are in general found to exhibit better physical properties, closer to their bulk, relaxed counterparts. For SiO2, $`E_p`$, determined from optical index, step coverage or optical transmittance measurements, is in the range 10–100 eV, depending on the particular set of parameters used in the growth process.

Clearly, bombardment profoundly affects the microstructure of the deposited material. For amorphous solids — which is the case of SiO2 films — static disorder hampers the measurement of detailed information beyond the second-neighbor shell. As a consequence, experimental results must be considered in the light of structural models. In this way, proper conclusions on the atomic mechanisms that mediate bombardment-induced densification can be drawn. Thus, for instance, models help understanding how changes in the medium-range structure of amorphous silica (a-SiO2) influence the Si–O–Si bond-angle, which in turn is related to such properties as the density, the optical index or the amount of built-in stress.

It is useful to discuss the structure of a-SiO2 in terms of short-range order (SRO) and medium-range order (MRO) correlations, at length scales $``$ 4 Å and in the range 4–10 Å, respectively. SRO is well defined in terms of such quantities as bond lengths and coordination numbers. It is mainly associated to the structural “building block” of the network, the Si(O1/2)4 tetrahedron. Both numerical simulations and experiment have shown the SRO correlations in permanently-densified (20%) a-SiO2 recovered from high pressures ($``$ 15-20 GPa) to be essentially the same as in the “normal” material.

In contrast, MRO correlations, which are at the origin of the so-called “first sharp diffraction peak” (FSDP), undergo important modifications upon compression: the height of the FSDP decreases and its position increases with density. Molecular-dynamics and Monte-Carlo simulations have attributed this effect to modifications in the topology of the network; in particular, the intensity of the FSDP is reported to be roughly proportional to the number of six-membered rings. (A $`n`$-membered ring is defined as a closed loop of $`n`$ Si–O bonds). Further, the x-ray diffraction spectra of porous silica thin films deposited at low substrate temperatures, using either electron-beam evaporation or low-pressure chemical vapor deposition (LPCVD), exhibit a strong low-angle scattering intensity and a FSDP height smaller than that of bulk a-SiO2. The latter effect, which diminishes upon annealing, has been attributed to differences at the MRO level. However, to the best of our knowledge, no systematic studies of the changes in MRO with the kinetic energy of the ions or the ion-current density have so far been reported; this constitutes one objective of the present work.

Several computational methods have been employed to assess the role of ion beams during thin film growth. Among them, molecular-dynamics (MD) simulations have provided useful insights into the understanding of ion-induced layer-by-layer growth, surface mobility, pore annihilation, stress and defect formation. These studies have focused on the role of collision-induced events for obtaining crystalline material. Although the precise values of the parameters required for this purpose are material-dependent, a beam energy of the order of a few tens of eV — in agreement with experiment — appears to be adequate, yielding ordered structures on the scale of many interatomic spacings. These simulations allow collision-induced events to be sorted out and relevant deposition parameters to be identified, thus providing useful information for better control of the growth process.

The situation is more complicated for the case of amorphous materials growth, as spatial order is short-ranged. Only very few MD studies of the structural properties of amorphous thin film growth have thus far been reported; notable exceptions are Refs.  for silicon and Refs.  for diamond-like carbon films. In the latter case, bonding is found to be partly $`sp^2`$ and partly $`sp^3`$, that is, between graphite and diamond. For the case of SiO2 — for which no such studies so far exist — considering the large number of polymorphs of the material, as-deposited films are expected to show a large diversity of structures as a function of density.

In view of this, and given the technological importance of the material, we have carried out detailed MD simulations of the growth of silicon dioxide on a-SiO2 substrates. More precisely, growth is modelled within a simplified IBAD scenario: we examine, at the atomic level, how the structure of the films evolves as a function of $`R`$, the ratio of medium- (30 eV) to low-energy (1 eV) SiO2 particle fluxes impinging on the surface. The rationale for the model is given in Sec. II. For the sake of simplicity, the energetic ions are taken to be SiO2 “molecules”; although this is a very crude approximation to the real growth process, our model is expected to provide much-needed information into the fundamental aspects of bombardment-induced densification. We find, on average, that densification is related to changes in the ring structure, in a way which is similar to the case of pressure-induced densification of “normal” (2.2 g/cm3) silica glass. The intensity of the FSDP and the number of six-membered rings are here found to increase with density. We also find that the crossover from low- to high-density films is associated with a tensile-compressive stress transition, occurring for $`E_p`$ of about 12–15 eV. It is observed, finally, that medium-energy bombardment promotes the annihilation of structural defects, namely three- and four-membered rings, which are at the origin of the $`D_2`$ and $`D_1`$ defect lines observed in the Raman spectrum of silica.

## II Computational details

MD consists in integrating the classical equations of motion for a set of interacting particles. The interactions can be derived from first principles, or expressed in some effective, empirical form whose parameters are fitted to experimental data. The first-principles approach is very accurate, but does not allow large systems to be dealt with over “long” times cales. These limitations are not so much of a problem in the empirical approach, but accuracy and “transferability”, that is, the ability for a model potential to adequately describe a phase to which it was not fitted, is an issue. For the problem considered here, empirical potentials are unavoidable. The model developed by Nakano et al. has been successfully employed to describe the bulk and the surface of disordered phases of SiO2, viz. amorphous and porous silica, which are evidently relevant to the present problem. In this model, silicon is assumed to be coordinated to four oxygen atoms, that is, the network is chemically ordered. The total potential energy of the system is the sum of two- and three-body contributions:

$`V={\displaystyle \underset{i<j}{}}V_{ij}^{(2)}(r)+{\displaystyle \underset{i,j<k}{}}V_{ijk}^{(3)}(\stackrel{}{r}_{ij},\stackrel{}{r}_{jk}),`$ (1)
where $`i`$, $`j`$, and $`k`$ run over all particles. The two-body term is written as

$`V_{ij}^{(2)}(r)`$ $`=`$ $`A_{ij}({\displaystyle \frac{\sigma _i+\sigma _j}{r}})^{\eta _{ij}}+{\displaystyle \frac{Z_iZ_j}{r}}e^{\frac{r}{\lambda }}`$ (3)
$`{\displaystyle \frac{\alpha _jZ_i^2+\alpha _iZ_j^2}{2r^4}}e^{\frac{r}{\xi }},`$
where the various terms represent the steric repulsion, the Coulomb interaction due to charge transfer, and a charge-dipole interaction due to atomic polarizability, respectively. The three-body term, which describes the covalent bonding, is expressed as

$`V_{ijk}^{(3)}(\stackrel{}{r}_{ij},\stackrel{}{r}_{jk})`$ $`=`$ $`B_i\mathrm{exp}({\displaystyle \frac{\mu }{r_{ij}r_0}}+{\displaystyle \frac{\mu }{r_{jk}r_0}})`$ (6)
$`\times [\mathrm{cos}(\theta _{ijk})\mathrm{cos}(\theta _{ijk}^0)]^2`$
$`\times \mathrm{\Theta }(r_0r_{jk})\mathrm{\Theta }(r_0r_{ij}),`$
where $`\mathrm{\Theta }`$ is the Heaviside step function, and $`\theta _{ijk}^0`$ is the equilibrium value of the $`(\stackrel{}{r}_{ji},\stackrel{}{r}_{jk})`$ bond angle. This term is calculated for nearest-neighbor Si–O–Si and O–Si–O triplets only. The model was fitted to various physical properties of bulk SiO2, such as cohesive energy, elastic constants and the phonon density of states; the values of the parameters can be found in Ref. .

As demonstrated by surface extended x-ray absorption fine structure and x-ray photoelectron spectroscopy measurements, the atomic structure of silica coatings strongly depends on the particular deposition process used: both silicon and oxygen atoms can be found in different environments, allowing the presence of homoatomic bonds. This is particularly true when sub-oxides are formed. These effects are not considered in our model; thus, the role of structural defects — such as peroxide radicals and, more generally, homopolar bonds — on relaxation mechanisms are not taken into account. For the case of (nearly) stoichiometric situations, experimental data on SiO2 produced by plasma enhanced chemical vapor deposition (PECVD), RF sputtering, secondary ion deposition, and ion-beam deposition indicate that the basic “building-block” of the network is the Si(O1/2)4 tetrahedron. Thus, chemical ordering is expected to dominate in these cases, and the model potential of Nakano et al. is assumed to correctly describe atomic correlations.

In IBAD, silicon (di-)oxide is evaporated in a low-pressure oxygen atmosphere; part of the condensing particles therefore have an energy smaller than 0.1 eV. Simultaneously, a beam of O+/O$`{}_{2}{}^{}{}_{}{}^{+}`$ ions with energy in the range 30–500 eV is directed onto the substrate, in order to promote surface relaxation as well as to hamper the formation of suboxides due to preferential sputtering of oxygen atoms. Although the experimental procedure is relatively simple, an atomistic description of it is a formidable task, as it requires an environment-dependent potential for oxygen ions, neutral gaseous particles, and for surface and bulk atoms. For simplicity, because large systems have to be dealt with, we describe the growth of silica in terms of a single type of condensing particle, namely SiO2, for both the deposited species and the incident energetic ions. Since the clusters are free to dissociate on the substrate (and actually do so), we expect this approximation to be of little consequence on the final structure of the deposited film. However, because oxygen atoms are preferentially sputtered as the deposited energy is raised (leading to a 10% deficit in the worst cases), SiO3 clusters were introduced when necessary so as to maintain the stoichiometry at 2:1. Also, the energy transferred during a collision varies with the mass of the incoming particle, which, in turn, may influence the value for which dense films are obtained. Thus, care must be taken when comparing these results with experiment. We also neglect atomic ionisation and electronic loss during collisions, a reasonable approximation at low energies, especially for the case of insulators. The present model is evidently crude, but we expect that the generic features of the structure will be relevant to real materials, even though the dynamics of the deposition process is not described precisely.

The substrate on which growth takes place consists of a $`27.59\times 27.59\times 10`$ Å3 amorphous slab containing 747 atoms, extracted from a bulk ($`27.59^3`$ Å3) sample prepared using the activation-relaxation technique (ART) of Barkema and Mousseau. The substrate is first subjected to an extensive thermal annealing cycle in order to allow the surface to relax. A 3 Å-thick layer at the bottom is held fixed, in order to mimic a semi-infinite system, but also to prevent the substrate from moving because of momentum transfer from incoming atoms. The remaining 7 Å are coupled to a heat bath at room temperature, consistent with experiment. Velocity renormalization was used to maintain the system at constant temperature.

SiO2 particles of roughly thermal kinetic energy ($`E=0.020.06`$ eV) were used in an attempt to first grow samples without bombardment. However, this led to unstable structures and to clusters of a dozen atoms desorbing from the growing film. Much better results were obtained with $`E=1`$ eV, an energy similar to that used by Hensel et al. (2 eV) in their MD study of bombardment-induced epitaxial growth of elemental silicon. For the high-energy particles, we used 30 eV, as suggested by the experimental results of McNeil et al. As noted previously, important process parameters are expressed in term of the ion-to-condensing particle fluxes ratio, $`\varphi _i/\varphi _n`$. However, this quantity is not known a priori ($`\varphi _n`$ is deduced from the measurement of the growth rate); thus, it is not a adequate parameter for monitoring the simulation procedure. Instead, we used, as other authors, the ratio, $`R`$, of low-to-medium energy particle fluxes. Low- and high-energy particles were chosen at random according to the desired $`R`$ value, and deposited at normal incidence onto the substrate. The equations of motion were integrated using the Verlet algorithm with an adaptative time step, so as to adjust to abrupt changes of the velocities during collisions. The time interval between the deposition of two successive particles was chosen in such a way as to provide a reasonable compromise between computational effort and proper relaxation of the surface. In practice, we found it adequate to relax the system during 4 and 8 ps after depositing a low- and a high-energy particle, respectively, before a new particle was introduced.

As the film grows, “normal” thermal diffusion into the substrate is not sufficient to drain away the excess energy brought about by incoming particles. Therefore, the entire sample was submitted to a 1 ps thermalisation phase each time before introducing a new particle. Prior to this, the temperature was around 350 K and 320 K for 30 eV and 1 eV SiO2 particles, respectively. Since the constants of diffusion are small at these temperatures, the equilibration procedure is expected not to interfere with the growth process in any significant manner.

Simulations were carried out for several values of the parameter $`R`$, namely 0, 0.08, 0.22, 0.38, 0.55, 0.73, 1.07, 1.53, 2 and 2.73. It can be determined a posteriori that these correspond to $`\varphi _i/\varphi _n`$ values of 0, 0.07, 0.18, 0.32, 0.41, 0.49, 0.62, 0.75, 0.80 and 0.9, respectively. In each case, about 1300 atoms ($``$430 clusters) were deposited. Following deposition, the samples were equilibrated during 20 ps at 300 K, after which a 40-ps microcanonical run was carried out in order to let the sample reach structural equilibrium. The corresponding “true” duration of simulated growth is about 3 ns. All calculations were performed using the program groF, a general-purpose MD code for bulk and surfaces developed by one of the authors (LJL).

## III Results and discussion

### A Structure of samples grown without bombardment

In order to assess the role and importance of bombardment on the properties of the grown material, we first discuss the case $`R=0`$, where all incoming particles have the same energy, viz. 1 eV. Figure 1 shows a snapshot of the system at the end of the run. The material is evidently highly disordered and porous. In Fig. 2, we plot, for this system, the corresponding integrated number of atoms from the bottom of the simulation cell; we note that it increases approximately linearly with $`z`$ in the region 13-53 Å. This portion will henceforth be referred to as the “bulk-like” part of the deposit, where structural analysis will be performed in order to limit artifacts due to the surface and the substrate. The same procedure is used for samples grown with bombardment, which leaves us with 20 to 40 Å of bulk material, about 1000 atoms in all cases.

The overall mass density in the bulk-like region of the slab, 1.3$`\pm `$0.1 g/cm3, is about 40% less than the “normal” density of a-SiO2 (2.2 g/cm3) and 12% less than the density of SiO2 grown by electron-beam evaporation on a silicon substrate at room temperature. The discrepancy between experimental and calculated densities is possibly due to the short duration of the simulation runs, which may not be long enough to allow complete relaxation, and thus full densification. (The computational effort required to achieve complete relaxation is formidable, and it constitues one of the outstanding problems of disordered-system simulations). For this model of silica, it is worth noting that the density is close to the critical value of 1.4 g/cm3, below which pore percolation occurs and leads to fracture. This was not observed here; however, as will be shown below, the film contains high tensile stress ($``$ 2 GPa). Thus, relaxation is expected to proceed by densification and/or fracture formation. This is consistent with the observation of reduced crack toughness and adhesion to the substrate for the case of coatings grown at low substrate temperature, $`T_S`$, with insufficient bombardment.

We now discuss the structure and topology of the deposited material, in order to provide a reference for samples grown in the presence of bombardment (cf. Sec. III B). Figure 3 shows the total pair distribution function for the bulk-like part of the sample, as well as for bulk silica; also shown are the partial distributions $`g_{\mathrm{SiSi}}`$, $`g_{\mathrm{SiO}}`$ and $`g_{\mathrm{OO}}`$. Here again, the amorphous character of the material is clearly evident: order beyond the nearest-neighbour shell rapidly fades out. The nearest-neighbor Si–O, Si–Si and O–O distances are listed in Table I, along with those of bulk a-SiO2, as well as other material produced by secondary ion deposition (SID) at room temperature and LPCVD at $`T_S=430^{}`$ C;  the calculated and experimental values are seen to be in satisfactory agreement. The partial pair correlations reveal a broadening of the first peak of the Si–Si distribution compared to the “normal” glass, presumably a consequence of the more disordered structure of the deposited film. Also, small pre-peaks to the left of the nearest-neighbour peak are seen to develop in the Si–Si and O–O correlation functions, at 2.45 and 2.32 Å respectively. As will be demonstrated below, these are due to the presence of two- and three-membered rings, leading to Si–O–Si and O–Si–O bond angles of 90 and 130 degrees, respectively.

In spite of the disordered nature of the network, the structural building block remains, to a good approximation, the Si(O1/2)4 tetrahedron, even though it might be slightly distorted. Indeed, the coordination of silicon atoms, calculated from the area under the first peak of gSiO, is about 3.9. Further, as displayed in Fig. 4(a), the O–Si–O bond angle distribution has a major peak centred at 109.7, the usual tetrahedral bond angle.

The topology of the network can be described more precisely in terms of rings of nearest-neighbor bonds, where an $`n`$-membered ring is defined as the shortest closed path of $`n`$ Si–O bonds originating on a given silicon atom.  In normal silica glass, the statistics are dominated by five-, six-, and seven-fold rings. Smaller rings are unlikely to occur because they are too expensive in terms of elastic bond-bending energy. The existence of “odd-numbered” rings results in a broad distribution of Si–O–Si bond angles. In the normal material, this distribution is centered around 143 and has a width of about 25, as can be seen in Fig. 4(b). In contrast, the mean Si–O–Si bond angle for evaporated SiO2 is known from infrared spectroscopy measurements to be 7–10 smaller than in a-SiO2. This is attributed to the presence of three- and four-membered rings, manifested by the presence of two sharp peaks at 606 and 495 cm-1 in the Raman spectrum of chemical vapor deposited (CVD) films, produced at low substrate temperature.

Such small rings are a common feature of disordered SiO2 networks, and they are a consequence of rapid quenching rates; thus, they are also expected to be present in SiO2 produced using evaporation methods. In the present model, the fraction of three- and four-membered rings is about 0.1 and 0.15 per deposited silicon atom, respectively, roughly ten and two times that of normal a-SiO2. The mean value of the Si–O–Si bond angle is 139. These results are consistent with the above-mentioned experimental results for evaporated films. In addition, two-membered rings (edge-sharing tetrahedra) are present here, but absent in the normal material. As can be seen in Fig. 4(b), the presence of such highly-strained rings gives rise to a sharp peak at low angles, $`92^{}`$, in the Si–O–Si bond angle distribution, while the distributions at high angles agree within statistical accuracy. The $`92^{}`$ peak is directly related to two-membered rings, as can be appreciated from Fig. 4(c), while three-membered rings give rise to a shoulder at about 130 in the bond-angle distribution.

### B Structure of samples grown with bombardment

We now examine the structure of samples produced in the presence of bombardment. It is expected that the energy provided by energetic particles will have an “annealing” effect, thereby reducing the proportion of high-energy structures — mostly two- and three-fold rings.

Simulations were carried out for several values of the $`R`$ parameter. As an example, we show in Fig. 5 the final configuration for $`R=0.73`$. Comparing this with Fig. 1, the main role of bombardment — densification — is clearly evident: the mass density of this new sample is $`\rho =2.16\pm `$ 0.1 g/cm3, much higher than the value of 1.3 g/cm3 for the sample grown without bombardment, and close to 2.2 g/cm3 for normal a-SiO2. (The density of $`\alpha `$-quartz is 2.65 g/cm3). In Fig. 6(a), we plot $`\rho `$ as a function of $`R`$; it is seen to increase rapidly, reaching a plateau at about 2.3 g/cm3 for $`R0.6R_c`$. This “critical” $`R`$ value roughly marks the crossover from a state of tensile stress (at small values of $`R`$) to a state of compressive stress, as can be seen in Fig. 6(b) where the pressure is plotted vs $`R`$. Near $`R_c`$, the ratio of ion-to-condensing particle fluxes $`\varphi _i/\varphi _n0.40.5`$; this yields a critical energy, $`E_c`$, for obtaining high-density, low-stress films, $`E_c=E_i\varphi _i/\varphi _n`$, on the order of 12–15 eV per condensing SiO2 unit.

The calculated value of $`E_c`$ is in reasonable agreement with experiment, considering the relatively crude nature of our model and the evident difficulty in measuring this quantity in a reliable manner. Indeed, values have been reported between 10 and 100 eV. Souche et al.,  for instance, gave 35 eV, based on optical index measurements on SiO2 prepared by oxygen IBAD at $`T_S=250^{}`$C, for a beam energy $`E_i=150`$ eV. McNeil et al., using a similar deposition process but lower $`E_i`$, 30 eV, found the optical transmittance to deteriorate with high ion current density, indicating that $`E_c`$ is smaller than 30 eV; an estimate in the range 10–20 eV has been given. Finally, Al-Bayati et al., using 10 eV Si+/O+ ion beams to grow SiO2 on a silicon substrate at 350C, found the interface to be smooth and free of strain, and the film to be of excellent quality, with no evidence of suboxides; thus, 10 eV would appear to be an upper bound to $`E_c`$.

The structural changes which take place upon densification can be investigated in terms of short- and medium-range order modifications. At short range, Si–O and O–O nearest-neighbour distances are found to depend very little on $`R`$, and they are close to the values observed for the sample grown without bombardment (cf. Section III A). We note, however, a slight ($``$2%) increase of the Si–Si distance, consistent with an increase of the Si–O–Si mean bond angle from 139 for $`R=0`$ to 143 for $`R=2.73`$. The Si coordination, between 3.90 and 3.96, shows no systematic variation with rising density.

However, the picture is quite different at intermediate length scales. In Fig. 7 we plot, for both $`R=0`$ and $`R=2.73`$, the total pair distribution function, $`g(r)`$, and the running coordination number, $`N(r)`$; the latter, essentially the integral of the former, gives the mean number of atoms within a sphere of radius $`r`$ centered on an average atom, irrespective of the atomic species. It is clear that the $`N(r)`$ curves for the two samples are very similar for $`r3`$ Å, confirming the role of the Si(O1/2)4 tetrahedron as the structural building block (no homoatomic bonds are allowed in this model). Above this value, however, strong differences develop. Densification, therefore, involves structural modifications at or beyond second nearest-neighbors.

Amorphous silica is known to possess significant order at medium range, manifested in the so-called “first sharp diffraction peak” (FSDP) at 1.55 Å-1, and observed in neutron and x-ray diffraction experiments. Although the precise origin of this peak, in terms of stable entities, is still the matter of debate, MD simulations by Nakano et al. have shown that its position and shape arise from correlations between atoms at distances as large as 12 Å. The FSDP can therefore be viewed as a signature of the presence of MRO in a-SiO2. In order to assess this picture within the present approach, that is, to quantify possible modifications in MRO upon densification, we calculated the static structure factor, $`S(q)`$, which is also accessible experimentally. This is shown in Fig. 8 for the samples with $`R=0`$ and $`R=2.73`$, as well as for a-SiO2 prepared using the ART method; the curves were smoothed using Gaussian functions of width 0.1 Å-1 so as to mimic thermal as well as experimental broadening. There is evidently very little difference from one sample to another for $`q`$-values larger than 4 Å-1, consistent with the existence of a structural unit common to all three networks. In contrast, the low-$`q`$ portion of the spectra exhibit significant differences. In particular, for $`q1`$ Å-1, the intensity varies strongly with $`R`$, in a manner which is inversely related to the density.

While the finite size of the system, plus periodic boundary conditions, artificially contribute to the intensity at such low-$`q`$ values, they cannot be solely responsible for the observed variations. Since all three samples are treated in the same manner and are similar in size, these variations must be related to the presence of pores in the low-density samples, an effect which can be understood in term of a Babinet argument — the scattering by a pore is similar to that from a cluster having comparable dimensions.

Changes in the MRO upon densification are also apparent in the FSDP, which the present model places at about 1.60 Å-1, in satisfactory agreement with the experimental value, 1.55 Å-1. For $`R=0`$, that is, highly porous material, the FSDP shows up as a mere shoulder, whose intensity is 23% smaller than that of a-SiO2. For $`R=2.73`$, a well-defined peak is observed, which is only 7% smaller than that of the normal material. The complete picture is provided in Fig. 9, where we plot the intensity of the FSDP as a function of density. Although the statistics are not perfect, there is a clear tendency for the height of the FSDP to increase with increasing density. It is interesting to note that if the densification were the consequence of a homogenous elastic compression of the network, that is, without structural modifications, the mean atomic distances would decrease with increasing density, implying an overall shift of $`S(q)`$ towards higher $`q`$ values. Here, although the FSDP position is modified, no obvious trend is found.

The differences in heights of the FSDP among the various samples are consistent with the results of Himmel et al., who found a 28% smaller height in the x-ray static structure factor of electron-beam evaporated silica, compared with bulk fused silica. Our calculations are also in qualitative agreement with MD simulations of pressure-induced densification of silica nanophases, where the height of the FSDP was found to be 15 % smaller compared to a-SiO2 in the density range 1.67 to 2.03 g/cm3.

The ring structure of the various models can tell us more about changes in MRO upon densification: because the size of rings is typically in the range 4–10 Å, modifications in the MRO are related to those in the arrangements of neighboring Si(O1/2)4 tetrahedra. The occurrence of the various types of rings as a function of $`R`$, relative to $`R=0`$, is shown in Fig. 10. In spite of statistical uncertainties, the following clear trends can be identified. (i) The number of two-, three- and four-membered rings decreases with rising $`R`$; indeed, as noted in Sec. III A, these correspond to high-energy structures and therefore are expected to anneal out during low-energy bombardment. (ii) The number of five-, six-, and seven-membered rings increases with rising $`R`$; these rings being close to the most probable size in a-SiO2, they are thus strongly favored. (iii) The number of nine-membered rings decreases with density, indicating a reduction in the number of pores. In contrast, the relative occurrence of eight-membered rings exhibits no clear trend, most likely because they are relatively stable entities in dense, cristalline forms of silica, such as $`\alpha `$-quartz (2.65 g/cm3) and coesite (2.92 g/cm3).

MD simulations on elemental silicon growth have demonstrated, in agreement with experiments, that low-energy ion bombardment has an effect similar to annealing. For the case of SiO2 prepared by IBAD, our results indicate that the elimination of two-, three- and four-fold rings may be achieved with the use of a second particle beam with kinetic energy in the range 12-15 eV. The concentration of three- and four-fold rings was reported to decrease for CVD SiO2 submitted to annealing; however, to the best of our knowledge, no such investigation has been reported regarding their dependence upon IBAD ion-current density or ion kinetic energy. Such data would be of utmost interest.

The modifications in the ring statistics were pointed out in MD and Monte-Carlo simulations of pressure-induced densification of bulk silica. It was found that most topological modifications are mediated by defect annihilation and bond switching. The former mechanism proceeds through the formation of bonds between non-bridging oxygen and three-fold coordinated silicon atoms, while the latter involves the diffusion of coordination defects (hopping of a dangling bond from one atom to a neighbor, accompanied by the formation of a new bond). This has also been identified as an important relaxation mechanism in a-SiO2, and it appears to be characteristic of chemically-ordered networks. The latter point being an assumption of our model, it is expected that these mechanisms also play an important role during bombardment-induced densification. More precisely, since coordination numbers do not show systematic variations with increasing density, bond-switching is expected to dominate over the annihilation of coordination defects.

In pressure-densified material, the evolution of the ring populations is quite different from the one we observe: the number of four-, eight-, and nine-fold rings increase, while those of five- and six-fold rings, as well as the intensity of the FSDP, decrease with increasing density (in the range 2.2–3.2 g/cm3). Taken together with the present results, this may indicate that in its most stable form, a-SiO2 has maximum MRO at ambient temperature and pressure.

## IV Concluding remarks

In this article, we have presented the results of extensive and detailed molecular-dynamics simulations of the growth of SiO2 on a a-SiO2 substrate, with a view of understanding the process of densification induced by low-energy particle bombardment. The physical picture which emerges from our observations is the following: the network undergoes significant density variations through structural changes at an intermediate length scale; more precisely, ring statistics are affected, while the structural unit of the network — the Si(O1/2)4 tetrahedron — remains essentially unchanged. Upon bombardment, small rings, which are associated with defects in the network, rearrange to produce larger ones, thus reducing the energy cost associated with bond bending and increasing the mean Si–O–Si bond angle. Large rings, particularly nine-membered ones, “disintegrate” into five-, six-, and seven-membered rings, thus reducing the occupied volume and causing the density to increase. Simultaneously, this influences the intensity of the first sharp diffraction peak (FSDP), which is found to increase with rising density. By calculating the pressure inside the sample, we have found an “optimal” incident energy range for low-stress, high-density films of 12 to 15 eV per condensing SiO2 particle.

The very brief ($``$4 ns) real-time equivalent duration of our simulations may not be sufficient to allow the samples to escape from local minima in the potential energy surface, leading to structures which are incompletely relaxed compared to the real material. Because relaxation processes, which are expected to increase MRO correlations, may occur within the time scale of an experiment, control over the film microstructure may be a difficult objective to attain when performing ion-assisted deposition experiments. In addition, most of MRO modifications arise from low-energy bombardment. For this reason, and as demonstrated by experimental results, we emphasize that any deposition process where the average ion kinetic energy is greater than about ten eV should be avoided. Besides ion-beam experiments, it was recently demonstrated that magnetron sputtering processes could produce low-energy high ion flux (up to 30 ions per neutral condensing particle), allowing amorphous silicon films with a tunable amount of MRO to be grown. We expect that such findings will stimulate futher experiments on SiO2 films, which will allow our results to be assessed. For example, a systematic study of the peaks associated with three- and four-membered rings in the Raman spectra of silica coatings, and the evolution of the FSDP, with the kinetic energy and flux of ions could help enhance our understanding of bombardment-induced densification.

Acknowledgments – A. L. is grateful to Gilles Dennler for useful discussions. Normand Mousseau is acknowledged for providing the ART model. This work is supported by grants from the Natural Sciences and Engineering Research Council of Canada (NSERC) and by the “Fonds pour la formation de chercheurs et l’aide à la recherche” (FCAR) of the Province of Québec. We are indebted to the “Réseau québécois de calcul de haute performance” (RQCHP) for generous allocations of computer resources.


format: text

text:
# The milliJansky 12𝜇m Population: First Followup Based in part on observations carried out at the European Southern Observatory, La Silla, Chile, and in part on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA

## 1 Introduction

The ISO satellite (Kessler et al., 1996) has vastly increased our observational capabilities in the mid- and far-IR, and has ushered in a new series of deep cosmological surveys at these wavelengths (eg. Oliver et al. 2000; Elbaz et al. 1999; Kawara et al. 1998; Puget et al., 1999; Clements et al. 1999). In the mid-IR these surveys reach flux limits more than 100 times fainter than our previous best surveys with IRAS (Hacking et al., 1987), whilst also attaining substantially improved spatial resolution. They thus represent an important new tool for observational cosmology. They are especially important since there is increasing evidence that much of the star-formation activity in the universe is obscured by dust, which re-radiates the energy absorbed at optical and UV wavelengths into the mid- and far-IR. While the discovery of luminous and ultraluminous infrared galaxies (LIGs and ULIRGs) in the local universe by IRAS was the first hint of the potential importance of obscured emission (see eg. Saunders & Mirabel, 1996, and references therein), it wasn’t until the discovery of the cosmic infrared background (CIB) (Puget et al., 1996) and large numbers of discrete sources in SCUBA submm surveys (eg. Eales et al., 1999) and ISOPHOT far-IR surveys (eg. Puget et al., 1999) that the cosmological significance became clear. It now seems that such obscured objects might contribute up to 2/3 of the integrated emission in the universe (Gispert et al., 2000) and might thus be the sites for the majority of star formation. Alternatively, a significant number of the obscured sources might be powered by active nuclei, and the sources uncovered by mid and far-IR surveys might thus tell us much about the formation of quasars, and thence the cosmic X-ray background (Almaini et al., 1999).

While much of the emission from obscured sources will be emitted in the far-IR, at wavelengths $`>`$50$`\mu `$m, there are substantial advantages in pursuing surveys in the mid-IR, at wavelengths between 10 and 20$`\mu `$m. These are primarily to do with the limitations of the current generation of instruments. Since ISO had only a 60cm mirror, the angular resolution attainable at long, far-IR, wavelengths is quite poor, leading to problems with confusion and with identifications for followup surveys. These difficulties are largely avoided in the mid-IR The mid-IR region was also the most sensitive on ISO, with improvements up to 100x over IRAS, while the far-IR typically only reached 10x greater depths, partly as a result of confusion noise (Puget et al., 1999). Despite operating at shorter wavelengths, the mid-IR is still capable of detecting obscured star formation regions. For example, the Antennae galaxy has optically prominent galactic nuclei, but mid-IR images, obtained with ISOCAM (Mirabel et al., 1998), show that the bulk of the star formation resides between these nuclei, in a region of high obscuration. The peak of the mid-IR flux also coincides with the peak of both the 850$`\mu `$m emission, as seen by SCUBA, and the peak of the molecular emission, as revealed by CO observations (Haas et al., 2000, and papers in Hughes & Lowenthal, 2000).

Mid-IR observations, however, do have some problems. These are principally the result of the complex spectral energy distributions (SEDs) in the mid-IR band. For starburst galaxies, which are likely to be the majority of the objects found in a deep mid-IR survey, the SED is dominated by the so-called Unidentified Infrared Bands (UIB), which produce prominent broad features at wavelengths of 6.2, 7.7, 8.6, 11.3 and 12.7 $`\mu `$m, and can vary in strength from object to object. Furthermore, there is also the possibility of emission and absorption line contributions, the strongest of which will be the Si absorption feature at 9.7$`\mu `$m. The effect of these features on cosmological surveys in the mid-IR has been extensively studied by Xu et al. (1998). They find that the effects of these features will be more pronounced the narrower the observational passband. While most deep cosmological mid-IR surveys have been conducted in the 15$`\mu `$m band, the present work uses the broader 12$`\mu `$m filter. The Xu et al. modelling clearly shows that this filter will be less susceptible to the effects of the complex SED.

The current results of deep mid-IR surveys with ISO are summarised by Elbaz et al. (1999; but see also Oliver et al., 2000; Serjeant et al., 1997; Clements et al., 1999; Serjeant et al., 2000). Substantial evolution in the number counts, beyond no-evolution extrapolations from the local population, has been found (eg. Elbaz et al., 1999). The meaning of these counts, though, is still unclear. Elbaz et al. (1998), for example, claim that a new population of objects becomes significant at flux levels of $`<`$ 1mJy at 15$`\mu `$m, producing a distinct bump in the number counts. They suggest that this is a new population of objects at z$`<`$1.5 contributing a substantial amount of the star formation history of the universe. Meanwhile, Xu (2000) suggests that a combination of simple luminosity evolution at a rate of L $`(1+z)^{4.5}`$ combined with the colour effects of the UIB features can produce the jump in number counts seen in the data. It is not possible to decide between these possibilities using the existing number count data. The only way to determine the true nature and significance of the sources uncovered by deep mid-IR surveys is to determine their redshifts, and thus to explicitly measure their luminosity function and evolution. Various models of this population can then be explicitly tested. We cannot yet directly obtain redshifts for mid-IR sources in the waveband of detection, but must instead use ground-based optical and/or near-IR observations. Identifications, and spectra, at these complementary wavelengths must thus be obtained.

We here report the first results from identification and redshift measurements for a deep 12$`\mu `$m ISO survey (Clements et al., 1999) conducted using the ISOCAM instrument (Cesarsky et al. 1996) on the ISO satellite. The survey covers a total area of 0.1 sq. degrees made up of four separate subregions separated by several degrees, and reaches a 3$`\sigma `$ flux limit of $``$300$`\mu `$Jy. A total of 148 sources were detected above the 3$`\sigma `$ limit, though a number of these may still be fragments of brighter sources. Where optical identifications in the US Naval Observatory Catalogue (USNO-A) are available, we can achieve star-galaxy separation using an optical-IR colour-colour method. This relies on the stellar spectra approximating a single temperature black-body, while the galaxies are assumed to have a substantial infrared excess. Very few of the galaxies found in the survey were previously known. Only one object, which was detected by IRAS, had a measured redshift. A full description of this 12$`\mu `$m survey is given in Clements et al. (1999) (hereinafter Paper 1). The goal of the present paper is to present the first results of the followup programme aimed at this survey.

The rest of this paper is organised as follows. In the next section we describe the sample of objects observed, and after that the observations themselves. The results of the spectroscopy are then discussed followed by brief comments on the identification of 12$`\mu `$m sources in two previously blank fields. Discussion of the luminosity function of the 12$`\mu `$m sources and then emission line diagnostics for their power sources follows. We then discuss the nature and significance of this population before drawing our conclusions. An H0 of 75kms-1Mpc-1 and $`q_0=0.5`$ are assumed throughout.

## 2 The Sample

The target list for the spectroscopic followup was constructed from the list of $`>5\sigma `$ non-stellar sources given in Paper 1. From these we then selected those objects with an optical ID brighter than R=19.6 magnitudes, as cross identified in the USNO-A catalogue. This provides a list of 25 objects. \[The object F3\_15, was found to be a star in paper 1 on the basis of its optical/IR colours. This stellar identification was inadvertently missed from Table 2 in that paper, but was included for the analysis work in that paper. It is thus not part of the basis sample in the current work.\]

As part of the INT Wide Field Survey (McMahon et al., 2000) and the UKIRT Mini-survey (Davies et al., private communication), much of Field 1 of Paper 1 was observed in the optical/near-IR, providing UBVRIJK images and photometry. The identifications and fluxes from USNO-A can thus be checked and are found to be broadly correct. The deeper optical data also allows us to go fainter in our identifications of ISO sources in Field 1. We find that all of the $`>5\sigma `$ ISO sources are identified in this field (the faintest is F1\_44 which is identified with an R band object with R=23.15). We also attempt to obtain identifications for fainter $`3`$ – 5$`\sigma `$ significance ISO sources using the new optical data. We find good identifications by eye for all except three ISO sources $`>3\sigma `$ within the area covered by the INT WFS image. These three unidentified sources are all of low signal-to-noise, (3.5, 3.0 and 3.2 $`\sigma `$) and so might be residual noise spikes or cosmic hits that have managed to escaped our ISOCAM reduction. The reliability of the source detection, discussed in paper 1, with only 4% of initially detected sources being rejected as spurious cosmic hits, suggests the that these may be real sources. We thus conclude that the optical imaging, to about R=23.5, which identifies 26 ISO sources and fails on these three, has a maximum identification incompleteness of $`11\%`$. This bodes well for obtaining complete identifications for the $`3`$ – 5$`\sigma `$ sources in the other three fields with deeper optical imaging, which has now been obtained (Clements et al., in preparation).

The basic properties of the sources observed in this project are given in Table 1.

## 3 Observations

Most of the observations were carried out at the ESO 3.6m telescope at La Silla from Nov. 13 to 15 1999. We used the EFOSC2 instrument. This provides both imaging and long slit spectroscopic capabilities, using filters and grisms (Patat. 1999). The detector is a 2060x2060 thinned Loral CCD which has good sensitivity throughout the optical spectrum. The pixel scale is 0.157”/pixel. Most of the spectra discussed here were taken using a 2” wide slit, though a few used a 1” slit, and Grism 12, which covers the spectral region from $``$ 6000 to 10000 Å at 2.12Å/pixel dispersion. A few objects were also observed in the blue using Grism 11 which covers $``$3400 to 7500Å at the same resolution. Wavelength calibration was achieved using a He-Ar arc lamp attached to the spectrograph, while flat fielding used a dome flat. The contribution of sky lines to the spectra was calculated by fitting and subtracting a cubic spline function to the spatial direction for each detector row in the dispersion direction. Standard IRAF tasks were used for these procedures. The CCD response function was calibrated out of the spectra using observations of the spectroscopic standard stars LTT9239 and LTT1020. Spectroscopy integration times ranged from 900 to 1800s. Typical observing conditions were for seeing between 1” and 2”, and with reasonable transparency. The conditions, however, were not sufficiently stable to be photometric.

Two objects, F1\_0 and F2\_0 were observed on a previous run, from 24 to 25 December 1997, which was severely affected by poor conditions and equipment problems. This was also conducted at the ESO 3.6m using EFOSC2. The observational setup was broadly similar, though a different grism, Grism 1, was used. This covers $``$3200 – 10000 Å at 6.3Å/pixel. The observations and data reduction for these objects followed the same scheme described above. The reduced resolution, however, means that these observations cannot be used for emission line diagnostics.

In addition to the spectroscopic observations, R band imaging was also obtained for two objects during the 1999 run. These were F4\_2, and F4\_11, chosen to be bright 12$`\mu `$m sources, with strong detections in the ISOCAM observations, but without a bright (R$`<`$19.6) optical identification. The images were taken as a series of four separate 200s integrations each shifted by -30” in both RA and DEC. The images were flat fielded using an imaging dome flat, and them median stacked with the appropriate shifts taken into account. Residuals left after the dome flat were them corrected be median stacking the images for both objects with each other without the shifts corrected, to provide a secondary, sky flat image.

## 4 Results

### 4.1 Spectroscopic Results

The reduced spectra were all visually inspected and initial classifications, as stars, emission line galaxies and, in one case, a possible absorption line galaxy, were made on this basis. The vast majority of objects turned out to be emission line galaxies. These spectra were further analysed to extract redshifts, determine emission line strengths, and to search for broad emission lines. Redshifts and derived properties are summarised in Table 2. The following method was used to extract spectral properties: first the brightest pixel in a given wavelength range is assumed to be H$`\alpha `$. Then a composite spectrum, including a background flux constant wrt. wavelength, H$`\alpha `$, \[NII\] and \[SII\] emission lines, is matched to the spectrum. The lines are assumed to have Gaussian shapes, with \[NII\] and \[SII\] lines having the same width but different from H$`\alpha `$. The strengths of the H$`\alpha `$, \[NII\] and \[SII\] lines are allowed to vary, but the ratios of the \[SII\] and \[NII\] doublets are fixed at their physical values (Osterbrock, 1989). The assumed redshift is also allowed to vary. The $`\chi ^2`$ between the model and the observed spectrum is minimised using a downhill simplex method (AMOEBA in IDL). Where a broad line is suspected, an additional model component consisting of an extra H$`\alpha `$ with larger initial Gaussian width is added to see if the $`\chi ^2`$ is significantly reduced. This is found to be the case in only one object. The results of these fits are given in table 3, as are the results of the spectroscopy more generally. Spectra are shown in Fig 1 and, for the one quasar in the sample, F4\_12, in Fig.2.

Given the redshifts we then calculate the 12$`\mu `$m luminosities of the objects in the usual manner. We include a K-correction to account for the effects of the complex spectral energy distribution (SED) in the mid-infrared. The K-correction adopted is a simple linear fit to the models of Xu et al. (1998). Over the redshift range of interest (0 – 1.2) the Xu K-corrections can be approximated as $`K=1.7z`$, using the definition of K-correction from Xu et al. (1998). However, this is the mean K-correction and does not take into account the effects of SED differences from one object to another. These differences, though, are smaller in the ISO 12$`\mu `$m band than in any of the other ISO filters used for deep ISOCAM surveys (Xu et al. 1998). The uncertainties produced by this scatter should be less than $``$0.2 mag, and thus should not have a major effect on the luminosity function we derive from these luminosities.

#### 4.1.1 Notes on Specific Objects

F3\_3: This object was classified as a galaxy based on its location on the B-R-F12 colour/colour diagram (Paper 1). However, spectroscopy clearly shows that the object is an M star. The previous misclassification is probably due to the large molecular absorption bands in the M-star spectrum suppressing optical emission below what would be expected from a simple black-body spectrum.
F1\_34: The H$`\alpha `$ emission line in this object appears to be broad, $``$1000 kms-1 (cf. 400kms-1 for the other, narrow lines) indicating the likely presence of a Seyfert 1 active nucleus.
F4\_12: At z=1.2 this is the highest redshift object in the sample. It is also the only object in the sample characterised by a clear broad line spectrum, containing emission lines of MgII and AlII/CIII\]. There also appears to be a broad emission feature to the blue of MgII which could be due to FeII. The MgII FWHM is found to be 4300$`kms^1`$ (observed frame). The spectrum is shown in Figure 2.

## 5 Blank Fields

Ten sources from Paper 1 are as yet optically unidentified. At the same time as our optical spectroscopy, we thus also obtained moderately deep images for two 12$`\mu `$m optical blank field sources, F4\_11 and F4\_2. These are shown in Figure 3. In both cases we find plausible faint optical identifications within 7” of the ISO source position. In the case of F4\_11 the identification is isolated and unambiguous. In the case of F4\_2, however, there would appear to be a group of similar faint objects close to the primary identification. There is thus the potential that the ISO source may be a combination of two, possibly interacting, objects, or that the identification might be one of the other optical sources in the field.

We can obtain an estimate for their magnitude by comparison with magnitudes of known objects imaged at the same time. We find that F4\_11 has R$``$20 and F4\_2 R$``$22. These results are consistent with the results from the INT WF survey data discussed above.

With identifications obtained for these two sources, there remain 8 sources with ISO detections $`>5\sigma `$ that still lack optical identifications. Since these are the objects with the highest 12$`\mu `$m to optical ratios, they are clearly interesting targets. Their optical faintness suggests that they are at high redshifts, and they are thus likely to be highly luminous. Further observations will be proposed to obtain identifications and redshifts for these, and fainter ISO sources.

## 6 The Luminosity Function

The data presented in this paper provides redshift information for a complete subsample of the 12$`\mu `$m survey. We are complete to the 5$`\sigma `$ flux limit of the survey, and to R=19.6 magnitudes for the optical identifications and redshift information. We can thus use the accessible volume method (Avni & Bahcall, 1980) to determine the 12$`\mu `$m luminosity function in this sample, whose 12$`\mu `$m flux limit is a factor of $`>`$ 20 deeper than existing determinations of the mid-IR luminosity function (Hacking et al., 1987). We use the standard 1/Vmax method to determine the luminosity function, with the volume determined from the maximum distance at which a source would be detected ie. the accessible volume. The luminosity function contribution of each object is 1/Vmax where Vmax is given by the redshift at which an object would drop out of either the survey data or the followup programme, whichever is smaller.

Thus:

$$\mathrm{\Phi }(L)dL=\underset{i}{}\frac{1}{V_i}$$
(1)
where $`\mathrm{\Phi }(L)`$ is the number of objects $`Mpc^1`$ in the luminosity range $`L`$ to $`L+dL`$ and

$$V_i=min(V_{opt,i},V_{12,i})$$
(2)
where $`V_{opt,i}`$ is the volume in which the $`i`$th object could have been detected down to the R=19.6 magnitude limit of the optical followup, and $`V_{12,i}`$ is the volume in which the $`i`$th object could have been detected down to the 500$`\mu `$Jy flux limit in the original survey. The volumes are calculated taking into account the variable flux depths of the original 12$`\mu `$m survey (see Paper 1).

The luminosity function produced from these considerations is shown in figure 4. As can be seen the luminosity function is a fairly close match to existing 12$`\mu `$m luminosity function determinations based on low redshift IRAS data (see eg. Fang et al 1998). The only substantial excursion from the local luminosity function occurs at the highest luminosities. Since these high luminosity bins are also dominated by the highest redshift objects this might be taken as evidence for evolution. Combining the two highest luminosity bins to provide better statistics (3 objects) at the highest luminosities suggests a $`3\sigma `$ significance difference between the current luminosity function and that derived from IRAS (Fang et al. 1998; Rush et al., 1993). The mean redshift for objects in this bin is 0.7, implying density evolution at a rate of $`(1+z)^8`$ or luminosity evolution at a rate of $`(1+z)^5`$ with large uncertainties on these exponents. These values are broadly consistent with the evolution in the ultraluminous IRAS galaxies suggested by Kim & Sanders (1998) and with that suggested by Xu (2000) for the 15$`\mu `$m population. However, we as yet have only a small number of objects at the highest luminosities. Indeed, the highest redshift bin contains only one object, F4\_12. This is one of only two broad-line objects in the survey as well as being at the highest redshift, so might be an unusual object biasing the statistics. There may also be additional biases entering the data, eg. Malmquist or Eddington bias, since the sample size is so small. These issues will be addressed when the identification programme is complete and we have a larger number of sources.

These hints at evolution in the high luminosity/redshift end of the population can only be confirmed by determining redshifts for the rest of the sample. Since the remaining galaxies are optically faint (R$`>`$19.6) there is a strong possibility that substantially more high redshift galaxies will be found, allowing a much better examination of the high redshift and high luminosity end of the luminosity function.

## 7 Emission Line Diagnostics

The nature of the ionising spectrum in a galaxy can be determined by examining emission line ratios (Baldwin et al., 1981). Standard ratios include H$`\alpha `$/\[NII\] and H$`\alpha `$/\[SII\], which can be determined for most of the objects in the current sample, as well as \[OI\]/H$`\alpha `$, \[OIII\]/H$`\beta `$ and \[OII\]/H$`\beta `$. These latter ratios are not available for most of the objects in the sample since our spectral range did not extend far enough into the blue, or go deep enough to detect the weak \[OI\]6300 line. Consideration of the standard emission line diagnostic plots shows that we can attempt a classification using just the H$`\alpha `$, \[NII\] and \[SII\] lines we have for the majority of the objects. Our classification diagram is shown in Figure 5. On this basis we find that the majority of our 12$`\mu `$m sources (11/16 where appropriate data is available) have HII region-like spectra ie. are powered by starbursts, while 4/16 have Seyfert2-like spectra ie. contain and may be powered by active nuclei. Two further objects have broad line Seyfert1-like spectra, producing a total of 6/18 objects that contain AGN. Three objects are not clearly classified by this technique and could be composite objects or LINERs.

We can extend this classification scheme for the four objects at high enough redshift and with sufficient signal-to-noise in their spectra to detect H$`\beta `$ and \[OIII\]. Of these objects only one is fully classified using just H$`\alpha `$, \[NII\] and \[SII\]. This object, F1\_12, is classified as a Sy2 by both methods. Of the other three objects where we can use H$`\beta `$ and \[OIII\], F1\_5 was previously unclassified due to a poor \[SII\] detection, and the others, F1\_11 and F4\_5, have mixed classifications. Once H$`\beta `$ and \[OIII\] ratios are included we find that F1\_5 and F1\_11 appear to be starbursts (ie. HII-region-like). F4\_5 appears to have a confused classification, with \[OIII\]/H$`\beta `$ vs. \[SII\]/H$`\alpha `$ indicating an HII region spectrum, and \[OIII\]/H$`\beta `$ vs. \[NII\]/H$`\alpha `$ indicating a Sy2 spectrum. On further examination it appears that the \[SII\] line strength may be underestimated due to a sky absorption feature at the extreme red of the optical band ($``$9300Å). We thus tentatively conclude that it is in fact a Sy2.

The final tally is thus 14 starbursts, 7 AGN (Sy1 and Sy2), 1 mixed and one unclassifiable. Of the objects where we have a firm classification, 1/3 (7/21) are AGN-type and 2/3 (14/21) are starbursts. We find no significant difference between 12$`\mu `$m luminosities of the two classes, HII and Seyfert, though the most luminous objects are from the Seyfert class.

The Rush et al. (1993) sample of local 12$`\mu `$m selected objects had an AGN fraction of 13%, so the larger fraction we have found here, 33%, may represent some evolution in the population. Alternative explainations, such as different K-corrections between AGN and non-AGN containing objects, greater prominance in the optical for AGN objects, are also possible. These possibilities will be investigated when the identification programme for the complete 12$`\mu `$m survey is finished.

## 8 Discussion

### 8.1 Physical Associations

If we examine the redshift distribution, we find several clumps of objects that appear to group in redshift. These include an apparent group around z=0.11 in field 1, and around z=0.17 in field 4. We can test the reality or otherwise of these associations by examining the positions of the objects on the sky. In field 1 we do indeed find that the objects F1\_0, F1\_7 and F1\_9 appear to have a close physical association. They form a small group covering a region $``$ 1 arcminute across, corresponding to $``$ 100kpc at z=0.11, using our assumed cosmology (see Fig. 6). These objects may thus be part of a small group or weak cluster. One other object in field 1 is at roughly the same redshift, F1\_10 at z=0.113. It is much further away on the sky, 7.5 arcminutes, corresponding to 750kpc. It might thus be part of some weak larger scale structure, but this is difficult to assess with the small numbers of redshifts we have. A simple analysis, using the luminosity function derived in section 6, suggests that the association of F1\_10 with the other three objects is only marginally significant, at $``$90% level.

In field 4 we find that F4\_9 and F4\_3 are close in both redshift (difference $`<`$0.0005) and on the sky (39 arcsecond), with a separation corresponding to 94kpc at z=0.17. They are thus likely to be a close pair of galaxies. F4\_1, in contrast, is significantly further away in both redshift (0.004 in z) and angular separation (4.8 arcminutes). It is thus unlikely that there is a strong physical association between F4\_1 and the other two objects.

It is thus possible that 5 of our 12$`\mu `$m sources appear in two groups, though better statistics are needed to determine if the correlation function of 12$`\mu `$m selected objects is any different from optically selected, normal galaxies.

### 8.2 Bolometric Luminosities

Spinoglio et al. (1995) found a correlation between 12$`\mu `$m luminosity and bolometric luminosity. The basis for this correlation was argued to be that the wavelength is a ‘pivot point’ between the optical emission, which rises more slowly than linearly with bolometric luminosity, and the mid- to far-IR, which rises more rapidly than linearly with bolometric luminosity. They also found that the constant of proportionality relating 12$`\mu `$m and bolometric luminosity is different for starburst and AGN sources — thus $`L_{bol}=14L_{12}`$ for starbursts and $`L_{bol}=5L_{12}`$ for AGNe, where $`L_{12}=\nu (12\mu m)L_\nu (12\mu m)`$ in solar units. If we assume that our optical classifications of these objects genuinely identify the source of their luminosity, then we can calculate the bolometric luminosity of the sources. These are given in Table 2. We find that the typical bolometric luminosity for our 12$`\mu `$m sources is around 10$`{}_{}{}^{10.5}L_{}^{}`$. The population of mid-IR sources so far uncovered would thus appear to have relatively moderate luminosities, and to not represent any substantial new population of dust-enshrouded starbursts or AGN. Our highest luminosity source, F4\_12, however, turns out to have a bolometric luminosity of 10$`{}_{}{}^{12.7}L_{}^{}`$, approaching that of Hyperluminous infrared galaxies. Since this is the highest redshift object in the sample, and our luminosity function analysis seems to indicate strong evolution at high redshifts or high luminosities, then one can speculate that this might be the first member of a lurking population of obscured high luminosity systems that we have yet to uncover. Completion of our identification programme can test this idea.

### 8.3 The Nature and Evolution of the 12$`\mu `$m MilliJansky Population

Low luminosity 12$`\mu `$m selected sources appear to be similar to the local samples. We find that about 1/3 of them contain AGNe while the remainder are powered by starbursts. These fractions are somewhat higher than those discussed by Rush et al. (1993). At the highest luminosities there are hints of strong evolution being detected, but these are currently dominated by small number statistics. With deeper observations we will not only be able to confirm the presence of evolution but also test claims based upon the 15$`\mu `$m surveys that a substantial new population of objects becomes significant at the $``$1mJy flux level (eg. Elbaz et al., 1998; Aussel et al., 1999). It should be noted that selection of objects at 12$`\mu `$m is much less affected by K-corrections due to the complex mid-IR spectral energy distributions of star-forming galaxies and AGN (Xu et al., 1998) than selection at 15$`\mu `$m. It has been suggested (Xu, 2000) that the bump in 15$`\mu `$m number counts at $``$ 0.4mJy is entirely due to this K-correction effect and not due to any new population. While the current survey is working at somewhat brighter fluxes we find no clear evidence for such a population. However, if such sources are optically faint, whether intrinsically or as a result of their redshift, they may well have been missed by the optical selection criterion. Completion of the optical identifications and redshift measurements for all the objects in the survey is needed to test these claims and is thus a high priority.

### 8.4 Future Prospects

The present paper is only the first step in following up the galaxies found in this 12$`\mu `$m ISO survey. We have here discussed the results of redshift measurements on 26 of the 36 galaxies with fluxes above the 5$`\sigma `$ flux limit of our survey, forming a complete subsample with R$`<`$19.6. While we are able to obtain interesting conclusions, they are clearly limited by small number statistics, and we are as yet unable to study objects with a large F12/Fopt ratio (eg. F4\_2 which has F12 of 3.6mJy, but an R mag. of $``$22). Completing the redshift observations for the 10 $`>5\sigma `$ sources so far unobserved is clearly a high priority. In field 1, where we have deeper optical data, nearly all of the 12$`\mu `$m sources are identified with objects brighter than R=23.2. If this is representative of the rest of the sample, then complete identifications and spectroscopy should be fairly easy to obtain. Perhaps more interesting is the prospect of extending our redshift observations to fainter 12$`\mu `$m fluxes. If we go to a 4$`\sigma `$ detection limit, rather than 5$`\sigma `$ we would include $``$ 45 more sources. All of these objects, again, are detected in our deep optical observations of Field 1, so such a study is eminently feasible, given deeper optical imaging and a moderate amount of time for spectroscopy. We will thus be able to substantially increase the number of redshifts we have for faint 12$`\mu `$m sources, and will be able to place the tentative conclusions of the present paper on a much firmer footing. Such observations will make an excellent parallel study to the 15$`\mu `$m surveys also underway (eg. Oliver et al. (2000); Elbaz et al. (1999)).

## 9 Conclusions

We have obtained spectroscopy for a complete R$`<`$19.6 subsample of our 12$`\mu `$m survey with the ISO satellite. The redshifts of these objects range from 0.032 to 1.2, with derived bolometric luminosities ranging from 109.3 to 1012.7 L. There is tentative evidence for strong evolution at the high luminosity and/or high redshift end of the luminosity function. We find that 1/3 of these sources contain an AGN, while the rest seem to be powered by starbursts. At this stage our conclusions are limited by small number statistics, but there are very good prospects for extending the identification programme for this survey, and thereby allowing stronger, more definitive conclusions to be drawn.
Acknowledgements

It is a pleasure to thank Scott Croom and Steve Warren for the provision of the INT Wide Field Survey optical data for Field 1, and the UKIRT minisurvey team for providing similar infrared data. The Digitised Sky Survey was produced at the STSCI under US Government Grant NAG W–2166. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The USNO-A survey was of considerable help, and we would like to express our thanks to all those who helped to produce it. We would also like to thank Herve Aussel, Dave Alexander and Matt Malkan for useful discussions. DLC acknowledges support from a PPARC postdoctoral position and, during the early stages of this work, an EU TMR Network post (FMRX-CT96-0068). We would also like to thank the anonymous referee for many useful comments.


text:
# Creating vortex rings and three-dimensional skyrmions in Bose-Einstein condensates

## Abstract

We propose a method of generating a vortex ring in a Bose-Einstein condensate by means of electromagnetically-induced atomic transitions. Some remnant population of atoms in a second internal state remains within the toroidal trap formed by the mean field repulsion of the vortex ring. This population can be removed, or it can be made to flow around the torus (i.e. within the vortex ring). If this flow has unit topological winding number, the entire structure formed by the two condensates is an example of a three-dimensional skyrmion texture.

03.75.Fi,05.30.Jp

With the recent controlled creation of vortices  in trapped atomic Bose-Einstein condensates (BECs), and evidence for vortex nucleation above a critical flow velocity , the dynamics of topological structures in weakly interacting superfluids has become an active subject of experimental research. In this paper we propose methods of preparing a vortex ring within a BEC, allowing the controlled study of a vortex structure away from surface effects . The vortex ring is created by means of driving electromagnetic (em) fields inducing transitions between internal atomic levels. The em field couples to the relative phases between the levels, in such a way as to transfer angular momentum to the atoms. With an appropriately chosen superposition of em fields the local axes of this effective rotation form a closed circular loop, producing a vortex ring in the BEC. We also consider the stability of the vortex ring, and identify a way to use a second atomic component of the BEC to increase the stability of the ring against collapsing under its own string tension. Conversely, the vortex ring acts as a trap for the second component, providing purely atom-optical confinement of cold atoms. Finally, the combined two-component structure may be identified as a three-dimensional (3D) skyrmion , shown in Fig. 1. This nonsingular texture is of interest in its own right, as a topological structure beyond the simple vortex.

Topological defects such as vortices and monopoles possess cores at which the order parameter is singular, and indeed such defects are characterized by winding numbers defined on paths or surfaces that enclose the singular cores. The simplest topological defect relevant to superfluid physics is that of a current flowing around a closed 1D path parametrized by the angular coordinate $`\varphi `$. The simple mapping $`\alpha (\varphi )=\varphi `$, from physical space into order parameter space, where $`\varphi `$ is defined on a 1D circle, has ‘winding number’ one, which will be unchanged if $`\alpha (\varphi )`$ is continuously deformed. Topologically nontrivial configurations are also possible without singular cores: Textures are defined by the way in which the compact order parameter space is ‘stretched over’ physical space. An SU(2) order parameter, for instance, takes values on a sphere. One can ‘puncture’ this sphere, stretch the resulting hole to infinity, and so ‘spread’ the order parameter space over a 2D plane. Such a mapping, or any continuous deformation of it, is known as a 2D (or ‘baby’) skyrmion. Interpreting the order parameter sphere as the spin orientation of a spinor condensate, one can realize a 2D skyrmion as a coreless point vortex (or in 3D, a coreless vortex line) , such as has been created experimentally at JILA .

The full, 3D skyrmion results from taking an $`S^3`$ order parameter space and stretching it over $`R^3`$ in a similar way . An $`S^3`$ order parameter space is afforded most simply by a two-component BEC, whose interactions effectively fix the total density $`|\psi _1|^2+|\psi _2|^2=\rho `$ of the two complex macroscopic wave functions $`\psi _1,\psi _2`$. The low energy degrees of freedom therefore reside on a 3-sphere of angles $`0<\alpha ,\beta \pi `$, $`0<\gamma 2\pi `$ (a 3-sphere is a 3D subset in 4D, such that $`_{n=1}^4x_n^2=R^2`$). For spherical polar coordinates in physical space, a 3D skyrmion is given by any continuous deformation of the mapping $`\gamma (r,\theta ,\varphi )=\varphi ,\beta (r,\theta ,\varphi )=\theta ,\alpha (r,\theta ,\varphi )=\lambda (r)`$ for monotonic $`\lambda `$ with $`\lambda (0)=0,\lambda (\mathrm{})=\pi `$. The topologically invariant winding number

$$W=\frac{1}{2\pi ^2}d^3x\mathrm{sin}^2\alpha \mathrm{sin}\beta det\left(\frac{\alpha ^i}{x^j}\right),$$
(1)
which obviously vanishes for the ground state of the two-species BEC, equals one for the skyrmion. (The topological invariance of $`W`$ is not hard to understand: because the determinant in the integrand is precisely the Jacobian for the change of variables from $`x^i=(x,y,z)`$ to $`\alpha ^i=(\alpha ,\beta ,\gamma )`$, the integral is simply the volume of the unit 3-sphere, regardless of the precise form of $`\alpha ^i(x^j)`$.)

We can obtain a physical picture of the 3D skyrmion by writing out the macroscopic wave functions:

$`\left(\begin{array}{c}\psi _1(𝐫)\\ \psi _2(𝐫)\end{array}\right)`$ $``$ $`\sqrt{\rho (r)}\left(\begin{array}{c}i\mathrm{sin}\alpha \mathrm{sin}\beta \mathrm{exp}(i\gamma )\\ \mathrm{cos}\alpha i\mathrm{sin}\alpha \mathrm{cos}\beta \end{array}\right)`$ (2)
$`=`$ $`\sqrt{\rho (r)}\left(\begin{array}{c}i\mathrm{sin}[\lambda (r)]\mathrm{sin}\theta \mathrm{exp}(i\varphi )\\ \mathrm{cos}[\lambda (r)]i\mathrm{sin}[\lambda (r)]\mathrm{cos}\theta \end{array}\right).`$ (3)
With this realization of the skyrmion, $`|\psi _1|^2`$ vanishes on the $`z`$ axis and at infinity, and is concentrated in a toroidal region. On the other hand, $`|\psi _2|^2`$ vanishes on the circle $`\theta =\frac{\pi }{2},r=\lambda ^1(\frac{\pi }{2})`$. Expanding $`\psi _2`$ about any point on this circle, we find $`\psi _2\delta r+i\delta \theta `$, indicating that the nodal circle of $`\psi _2`$ is a vortex line. Hence component $`\psi _2`$ forms a vortex ring, within whose core $`\psi _1`$ resides and flows azimuthally with winding number one. Of course, any continuous deformations of Eq. (2) will still be a skyrmion, and the two-fluid interpretation of the skyrmion will always provide a useful intuitive picture. Moreover, unlike in the case of wave function monopoles , changing the atomic basis that defines the two fluids does not lead to a different picture of the skyrmion, but simply produces the same combination of vortex ring and current, with a different spatial orientation.

We propose a method of creating a 3D skyrmion in a trapped BEC using em fields, in a scheme motivated by viewing the skyrmion as this combination of vortex structures. Since several good methods have been proposed for generating simple vortices, and some of these have already been successfully implemented, engineering the state of $`\psi _1`$ presents no new difficulties of principle. Creating the vortex ring state in $`\psi _2`$ is more of a challenge. Although vortex rings may be formed through instabilities, as of a dark soliton in an insufficiently narrow trap , a more controlled method of preparing them may also be desirable, in order to study their rich range of motion and excitation.

We consider em-induced coherent scattering of atoms between different internal levels. By superposing different driving em fields we can construct a Rabi amplitude with a spatially-dependent phase profile. In particular, nodes of the em field amplitude may be topological singularities of the em phase. The coupling of the em field to the phase of the coherent matter wave allows the topological singularities of the em fields to be imprinted on the matter field. We demonstrate the preparation process of topological objects by numerically studying this coherent matter-wave dynamics.

The dynamics of the BECs with the Rabi coupling between the levels $`|i`$ and $`|j`$ follows from the coupled Gross-Pitaevskii equation (GPE)

$$i\mathrm{}\dot{\psi }_i=(H_i^0+\delta _i+\underset{k}{}\kappa _{ik}|\psi _k|^2)\psi _i+\mathrm{}\mathrm{\Omega }_{ij}^{}\psi _j.$$
(4)
Here the kinetic energy and the trapping potential are introduced in $`H_i^0`$:

$$H_i^0\frac{\mathrm{}^2}{2m}\mathbf{}^2+\frac{1}{2}m\omega _i^2(x^2+\alpha _i^2y^2+\beta _i^2z^2).$$
(5)
We have also defined the interaction coefficients $`\kappa _{ij}4\pi \mathrm{}^2a_{ij}N/m`$. Here $`a_{ii}`$ denotes the intraspecies scattering length in internal level $`|i`$ and $`a_{ij}`$ ($`ij`$) stands for the interspecies scattering length. The Rabi frequency, $`\mathrm{\Omega }(𝐫)`$, describes the strength of the coupling between the two internal levels. The total number of BEC atoms, the atomic mass, and the detuning of the em fields from the resonance are denoted by $`N`$, $`m`$, and $`\delta _j`$, respectively.

Several papers  have proposed methods of injecting vorticity into BECs by means of the interaction between em and matter fields. In our adaptation of these methods to create a 3D skyrmion, the simplest part is the step of producing the ‘inside’ component $`\psi _1`$ with a quantized circulation about the $`z`$ axis. This can be achieved by transferring population from an initial internal state $`|0`$ using the Rabi frequency $`\mathrm{\Omega }_{01}(𝐫)`$,

$$\mathrm{\Omega }_{01}(𝐫)=\mathrm{\Omega }_0[\mathrm{sin}(k_1x)i\mathrm{sin}(k_1y)]\mathrm{cos}(k_2z),$$
(6)
with a similar set of parameters as in Ref. . In particular, $`k_1^1`$ and $`k_2^1`$ are both larger than the size of the sample. The purpose of the $`\mathrm{cos}(k_2z)`$ is to confine $`\psi _1`$ along the $`z`$ direction. Note that in this scheme the Rabi frequency could represent a two-photon transition, with one photon from two orthogonal standing waves in the $`xy`$ plane and the other from a standing wave along the $`z`$ direction. Alternatively, the $`\mathrm{cos}(k_2z)`$ could simply represent collimation of the $`xy`$ beams, which would drive a single photon transition by themselves. The topological singularity in the em field at $`x=y=0`$ results in a single quantum of circulation about this axis, in the component $`\psi _1`$, as displayed in Fig. 2. One may also generate additional terms in Eq. (6). For instance, the standing wave along the $`z`$ axis could be further collimated in the $`xy`$ plane, limiting the spatial extent of $`\psi _1`$.

We propose two different methods of creating the more complicated part of the skyrmion, the closed vortex ring in $`\psi _2`$. The first is an adaptation of the method of Ref. , and makes use of two internal levels and the circulating state $`\psi _1`$ once it has been generated. To add the vortex ring requires a Rabi field coupling from state $`|1`$ to state $`|2`$, which is a phase-coherent superposition of a standing em field and a Gaussian beam, both along the $`z`$ axis:

$$\mathrm{\Omega }_{12}(𝐫)=\mathrm{\Omega }_0\left[1\eta ^2\mathrm{exp}[\rho ^2/\xi ^2]+iA\mathrm{sin}(k_3z)\right],$$
(7)
with $`A=2\eta /(k_3\xi )`$ and $`\rho (x^2+y^2)^{1/2}`$. We are interested in $`k_3|z|1`$, $`\rho \xi `$, and $`\eta 1`$ in such a way that $`|k_3(\rho \eta \xi )|1`$. Then we obtain

$`\mathrm{\Omega }_{12}(𝐫)/\mathrm{\Omega }_0`$ $``$ $`\rho ^2/\xi ^2\eta ^2+iAk_3z`$ (8)
$``$ $`2\eta /\xi (\rho \eta \xi )+iAk_3z`$ (9)
$`=`$ $`Ak_3[(\rho \eta \xi )+iz].`$ (10)
The phase of $`\mathrm{\Omega }_{12}(𝐫)`$ has a form of the quantized circulation with unit winding number around the closed ring $`\rho =\eta \xi `$ and $`z=0`$. Alternatively, the Gaussian beam resulting in the exponential term in Eq. (7) may also be replaced by two standing em fields: $`[\mathrm{cos}(2x/\xi )+\mathrm{cos}(2y/\xi )]/2`$.

For this first method, the field of Eq. (7) is tuned away from resonance. A localized shifting field is then applied as in Ref., to ‘turn on’ the transition in a small region, which is then swept around the trap to create a full vortex ring, and also undo the azimuthal phase winding of the $`|1`$ atoms, so that the population shifted into $`|2`$ will have no angular momentum about the $`z`$ axis. While this scheme takes advantage of an already demonstrated technology, it has the disadvantage that it tends to produce a highly excited skyrmion state, which may (as discussed below) then evolve self-destructively.

We therefore propose, in our second method, to couple atoms with em fields from level $`|0`$ to level $`|2`$, at the same time as we generate the simple vortex state in level $`|1`$. We need a Rabi field coupling from $`|0`$ to $`|2`$, of the same spatial form as $`\mathrm{\Omega }_{12}`$ of Eq. (7), as well as the $`\mathrm{\Omega }_{01}`$ field of Eq. (6). This scheme requires control of an additional internal state $`|0`$, but with a weak coupling it allows adiabatic transfer of the atomic populations , producing a ground state skyrmion.

For computational simplicity we demonstrate the three-level preparation process by considering two short nonadiabatic Rabi pulses: The first Rabi pulse of the form (7) generates the vortex ring in level $`|2`$ from the initial ground state population in level $`|0`$, and then the second pulse with the configuration (6) generates the trapped current in level $`|1`$ from the remaining atoms in level $`|0`$. Our numerical approach uses three states and can still generate an excited skyrmion, but it is amenable to fully 3D simulation, and our simulation indicates that the method can work. A careful study of an adiabatic population transfer will require a realistic model for the dissipation of the BECs, as well significant computational resources, and must be left for future investigation.

In our numerical simulations we assume that the traps are isotropic, $`\alpha _i=\beta _i=1`$, and that the trapping frequences are equal for all components. We choose the wave number of the em fields $`k_1=k_2=k_3=2\pi /(8l)`$, with $`l[\mathrm{}/(m\omega )]^{1/2}`$, $`\delta _i=0`$, $`\eta =0.04`$, $`\xi =30l`$, and $`\kappa _{ij}=100\mathrm{}\omega l^3`$. The number of atoms is determined by the previous relations as $`N=25l/(\pi a)`$.

We integrate the GPE on a spatial grid of $`128^3`$ points, using a split-operator method. The stability of the nonlinear evolution is improved by simultaneously propagating multiple copies of the wave function. Relatively small number of time steps (less than one thousand) provide a sufficient accuracy making the integration executable. Even with up-to-date workstations the optimized simulations can take several hours. Figures 1-3 represent the 3D skyrmion in the case where we have first created the vortex ring in level $`|2`$ by applying the Rabi frequency $`\mathrm{\Omega }_{12}(𝐫)`$ \[Eq. (7)\], with $`\mathrm{\Omega }_0=5000\omega `$ for $`0<t\omega 0.04`$. This is followed by creating a vortex line in level $`|1`$ by the driving $`\mathrm{\Omega }_{01}(𝐫)`$ \[Eq. (6)\], with $`\mathrm{\Omega }_0=80\omega `$ for $`0.04<t\omega 0.05`$. With the chosen set of parameters a large population remains in $`|0`$. This remnant population may be coupled out of the trap immediately after the em pulses as in the recent vortex experiments . Due to the shortness of the pulses the interactions of the atoms in level $`|0`$ do not have time to affect the structure of the skyrmion. On the other hand, for an efficient adiabatic population transfer the nonlinearity can be a problem. However, the techniques to overcome this were developed in Ref. .

The effort of creating such an intricate structure as a vortex ring, let alone a 3D skyrmion, is obviously only justified if its structure can then be observed. As with the filled vortex lines of JILA, the interior BEC in state $`|1`$ will enable detection of the density and phase profile of the skyrmion. An additional advantage of the skyrmion, as opposed to an empty vortex ring or to a vortex ring filled with a nonrotating interior BEC, is that it is less unstable against shrinking to zero radius. The kinetic and interaction energy of a vortex line gives it an effective string tension tending to shorten the vortex line, and in the case of a vortex ring, the ring’s radius $`R`$ (not its core thickness!) is energetically unstable to shrinking. In an infinite, homogeneous background, and in the absence of dissipation, the first-order dynamics of a vortex line means that a vortex ring will not actually shrink, but will instead move along its axis, at a velocity determined by its radius. But if its radius approaches its core thickness, a vortex ring can be annihilated (it is not topologically stable), and an excited vortex ring may also possibly twist or pinch, and break up into smaller rings. Moreover, motion at a constant speed cannot continue in a finite, inhomogeneous BEC. Dissipation is also present in experimental conditions, and although straight vortices have been found to have a long lifetime before floating out of the BEC cloud, the energy slope due to string tension may be much steeper than that due to vortex buoyancy (on the healing length scale rather than trap scale). Filling the vortex core lowers the string tension, and hence extends the lifetime of the filled vortex structure in comparison with that of an empty ring.

Although in other field theories skyrmions are often unstable against shrinking to zero size without the topological change of defect nucleation, the separate conservation of both atomic species in our case forbids this decay channel. This does not prevent the vortex ring from collapsing, but it does mean that the component $`|1`$ filling the ring cannot be eliminated. Furthermore, since the filling has one unit of angular momentum, there is a centrifugal barrier inhibiting the shrinking. As it tends to pin the vortex ring down, the ring may only collapse via the tunneling of atoms through $`\psi _1`$. Although the vortex rings are dynamically stable , numerical evolution in imaginary time shows that filled vortex rings, and even skyrmions, are still energetically unstable, on long enough (imaginary) time scales. Skyrmions are topologically stable as long as the order parameter is well defined everywhere; but they can be destroyed by the nucleation of a defect. Our simulations indicate that, while the energy decreases slowly but monotonically, the vortex ring in $`\psi _2`$ can pull through the state $`|1`$ if its inner radius approaches the healing length, so that a healing-length radius circle of zero total density nucleates. However, the decay via the tunneling of atoms through $`\psi _1`$ becomes increasingly slow with the larger occupation in $`|1`$ . A vortex ring with a nonrotating filling can, if its radius becomes too small, annihilate itself without significantly lowering the total density at any point.

Without dissipation, numerical evolution in real time shows that skyrmions as well as empty vortex rings do move through the background cloud, but instead of simply exiting the cloud as they approach its surface, they generate a background flow in the rest of the cloud which is able to draw them back towards the center. Much more extensive numerical studies will be needed to determine the stability of the resulting oscillation. Further analytical and numerical studies of skyrmion and vortex ring dynamics are clearly warranted, but in this paper we have shown that engineering and study of these advanced BEC structures should be feasible with current techniques.

This research was supported by the NSF and by EPSRC. We thank C.S. Adams, J. Dziarmaga, H.T.C. Stoof, and W.H. Zurek for valuable discussions, and gratefully acknowledge the hospitality of ITAMP at the Harvard-Smithsonian Center for Astrophysics.