text:
# Modifications of the hydrogen bond network of liquid water in a cylindrical 𝐒𝐢𝐎_𝟐 pore

## 1 Introduction

The increasing interest in the study of structural and dynamical properties of confined water in different environments arises from the close connection with many relevant technological and biophysical problems . It is well known that the properties of liquid water at ambient conditions are mainly determined by the microscopic local tetrahedral order . Understanding how the connected random hydrogen bond network of bulk water is modified when water is confined in small cavities inside a substrate material, is very important, e. g., for studies of stability and enzymatic activity of proteins, oil recovery, or heterogeneous catalysis, where water-substrate interactions play a fundamental role. The modifications of the short range order in the liquid depend on the nature of the water-substrate interaction, i.e. hydrophilic or hydrophobic, as well as on its spatial range and on the geometry of the substrate. Many experimental studies of structural and dynamical properties of water in confined regimes have been performed . Water-in-Vycor is one of the most extensively studied systems. In particular, it has been found in a recent neutron scattering experiment  that the hydrogen bond network of water in Vycor is strongly distorted with respect to bulk water. In order to shed light on this behavior, which was not evident in previous experiments , a more extensive study of the microscopic details of the molecular arrangement is needed.

Computer simulation methods allow a more detailed exploration of the molecular arrangements inside the system than most experimental techniques. They can therefore be considered as an ideal microscopic tool to study the properties of water in confined geometries, subject, of course, to the usual limitations of the method and the interaction potential functions employed. In this paper we present results obtained by a computer simulation study of water confined in a model for a Vycor glass pore. The simulation cell is obtained by carving a single cylindrical cavity in a modeled silica glass, thus reproducing the symmetry, the average dimensions, and the structure of the rough surface of Vycor pores at the atomic level.

In the next section we will briefly review the method for building the simulation cell and the implementation of the molecular dynamics. In the following sections we summarize some of the key results and compare them with recent experiments. The discussion is more detailed than in earlier work by our groups .

## 2 Model for the Vycor pore and molecular dynamics calculations

A glass of $`\mathrm{SiO}_2`$ is obtained by computer simulation, following the method proposed by Brodka and Zerda , using potential parameters given by Feuston and Garofalini . The system is simulated as a simple ionic model for $`\mathrm{Si}^{4+}`$ and $`\mathrm{O}^2`$ ions . Inside a cell of vitreous SiO2 with a box length of approximately 71 Å a cylindrical cavity of 40 Å diameter, corresponding to the average pore size of the Vycor glass of interest, is created. The details for generating the cavity with a rough surface and microscopic structure representative for pores in Vycor glass are given in previous publications . On the surface of the pore two types of oxygen atoms can be distinguished, depending on the number of silicon atoms to which they are connected. Nonbridging oxygens (NBO’s) are bonded to only one silicon atom and are saturated with a hydrogen atom. The number density of the hydroxyl groups in our samples is 2.5 $`\mathrm{nm}^2`$, in good agreement with the experimentally determined value . The second species of surface oxygen atoms are bridging oxygens (BO’s), which are connected to two silicon atoms.

For the water-water interaction we assume the SPC/E model . The atoms of the substrate are allowed to interact with the water sites by means of the Coulomb potential, where different fractional charges are assigned to BO’s ($`0.629|e|`$), NBO’s ($`0.533|e|`$), silicon atoms ($`1.283|e|`$) and surface hydrogen atoms ($`0.206|e|`$). In addition both BO’s and NBO’s interact with the oxygen sites of water via a Lennard-Jones potential, whose parameters are $`\sigma =2.70`$ Å and $`\sigma =3.00`$ Å for BO and NBO atoms, respectively, and $`\epsilon =230`$ K in both cases . During the simulation the Vycor atoms are kept fixed, and periodic boundary conditions are applied along the $`z`$-axis. The shifted force method is used with a cut-off at 9 Å. As discussed in a previous paper , the use of larger cut-off or Ewald summation does not change the nature of the obtained results. The temperature is controlled by means of a Berendsen thermostat with relaxation times between 0.5 and 1 ps. This procedure simplifies running the simulations over long times but has negligible effects on the calculated properties. All calculations have been performed at room temperature for several hundred picoseconds, after an equilibration run lasting more than 50 ps. More details can be found in Ref. .

## 3 Structural properties of the confined water

The molecular dynamics calculations have been performed for different numbers of water molecules, corresponding to different levels of hydration. The density profiles of the oxygen and hydrogen atoms of water for 96% (2600 water) and 19% degree of hydration (500 water) are shown in Fig. 1. Density profiles and further results of several intermediate degrees of hydration have been published elsewhere . We observe already at low hydration the presence of a layer of water molecules wetting the substrate surface; at nearly full hydration two layers of water with higher than average density are evident. The same features are present in the density profiles for the hydrogen atoms. Few molecules are able to penetrate into the glass and are trapped inside small pockets close to the surface, which are a byproduct of the ’sample preparation process’. The molecules in these pockets can easily be identified in the snapshot along the pore axis in Fig. 2. At low hydrations the surface is not covered completely; adsorption occurs preferentially in several regions, in which small clusters are formed. We have made no attempt to identify specific surface structures that favor cluster formation.

In the comparison with experiment the definition of the hydration level deserves some consideration. Full hydration in the experimental sample preparation corresponds to an estimated water density which is $`11\%`$ less than the density of water at ambient conditions . If the simulated cavity would be a perfect cylinder of radius 20 Å, this would translate into $`N_W=2670`$ molecules. Due to the roughness of the glass surface and the existence of the small pockets inside the glass close to the surface, which trap several water molecules over the entire simulation length, the exact number of molecules can only be determined by a grand canonical simulation scheme. For the investigated number $`N_W=2600`$ we observe indeed that the density profile does not reach the bulk value in the center of the pore. If one assumes that full hydration corresponds to a bulk-like region in the center of the pore, we estimate that an additional 100 molecules would correspond to full hydration, not much different from the experimental value. Thus, $`N_W=2600`$ corresponds to a hydration level of approximately 96%.

We calculate the site-site radial distribution functions for $`N_W=2600`$, i. e., at approximately the experimental density. The calculation of these functions is not straightforward when dealing with a non-periodic system, since excluded volume effects must be carefully taken into account . Unlike a bulk liquid, a *uniform confined fluid*, i.e., a collection of non-interacting particles in a confining volume, where periodic boundary conditions are only along the $`z`$-axis, does not create a *uniform* radial distribution function $`g=1`$. Rather, the form of the radial distribution function depends on the geometry of the confining system.

If we disregard the confinement effects, the raw site-site pair correlation functions can be calculated in the usual way from the simulation by taking the average number $`n_{\alpha \beta }^{(2)}(r)`$ of sites of type $`\beta `$ lying in a spherical shell $`\mathrm{\Delta }v(r)`$ at distance $`r`$ from a site of type $`\alpha `$ and normalizing it to the average number of atoms in the same spherical shell in an ideal gas at the same density

$$g_{\alpha \beta }^{(MD)}\left(r\right)=\frac{n_{\alpha \beta }^{(2)}(r)}{\frac{N_\beta }{V_p}\mathrm{\Delta }v(r)},$$
(1)
where $`N_\beta `$ is the total number of $`\beta `$ sites, and $`V_p`$ is the volume of the cylindrical cell. The three resulting raw site-site correlation functions are shown in Fig. 3. It is evident that the three functions do not approach a constant value at large $`r`$. The normalization factor in (1) must also include a uniform profile,

$$g_u(r)=\frac{V_p}{\left(2\pi \right)^3}d^3QP_{cyl}\left(Q\right).$$
(2)
$`P_{cyl}(Q)`$ is the form factor of the cylindrical simulation cell of height $`L`$ and radius $`R`$

$$P_{cyl}\left(Q\right)=_0^1𝑑\mu \left[j_0\left(\frac{\mu QL}{2}\right)\right]^2\left[\frac{2j_1\left(QR\sqrt{1\mu ^2}\right)}{QR\sqrt{1\mu ^2}}\right]^2$$
(3)
where $`j_n(x)`$ are the spherical Bessel functions of the first kind of order $`n`$. The resulting $`g_u(r)`$ is the smooth dashed line shown in Fig. 3.

The properly normalized oxygen-oxygen and oxygen-hydrogen site-site distribution functions are then obtained by dividing $`g_{\alpha \beta }^{MD}(r)`$ by $`f_cg_u(r)`$. The correction factor $`f_c`$, which lies between 1.0 and 1.04, takes into account the roughness of the surface and has been adjusted in such a way that the corrected pair correlation functions are close to one in the region $`r>8`$ Å. The corrected pair correlation functions are shown in Fig. 4 and are compared to the corresponding functions of SPC/E water at ambient conditions (dashed lines). We note that the modification of the oxygen-oxygen function relative to bulk water shows a similar trend as the experimental one (see Fig. 7a in Ref. ). The first minimum becomes shallower and fills in. Such behavior was also observed in experiments and computer simulations of water under pressure, which was explained as resulting from a collapsing hydrogen-bond network . In the oxygen-hydrogen function we notice that the first peak is lower in amplitude and the first minimum shifts towards shorter distances, as was observed in the experiment (Fig. 7b in Ref. ). Agreement with experimental results is less satisfactory in the case of the hydrogen-hydrogen pair correlation function. The experimental function is changed much more than the simulated one. The differences may be due to the fact that, in the experiment, the hydroxyl hydrogen atoms of Vycor are not distinguished from water hydrogen atoms.

We conclude that our results are in qualitative agreement with the experimental trend for the confinement-induced changes in water structure. According to Ref. , those changes indicate a distortion of the H-bond ordering in the confined water, induced by the interaction with the substrate. The distortion of the hydrogen bond network deserves further investigation. Since there are several approximations made in extracting the water-water contribution from the measured structure factors, computer simulations can contribute to our understanding by first providing a more detailed analysis of the simulated hydrogen bond network than is possible in experiment, and second by investigating the differences between water molecules located at various distances from the substrate surface.

## 4 Hydrogen bond network of confined water

In a first attempt to characterize the hydrogen bond network, we calculated the number of hydrogen bonds between water molecules and between water molecules and surface atoms. A geometric criterion is applied. A hydrogen bond exists between two oxygen atoms, if the angle between the intramolecular O-H vector and the intermolecular O$`\mathrm{}`$O vector is less than $`30^{}`$, provided that the O$`\mathrm{}`$O separation is less than $`3.35`$ Å. In Fig. 5 we show the results for $`N_W=2600`$. Close to the surface the number of water-water hydrogen bonds is reduced; a large number of interfacial water molecules is engaged in hydrogen bonds between water molecules and oxygen atoms of the pore surface. Taking the sum of both contributions, one can see that the total number of hydrogen bonds remains constant and the loss of water-water hydrogen bonds is compensated by water-Vycor bonds. Only bridging oxygen atoms of the Vycor surface contribute significantly to the number of hydrogen bonds; the amount of bonds between hydroxyl groups on the surface and water molecules is small and therefore not shown in Fig. 5. It is worth noticing that at this high level of hydration the average number of H-bonds is always lower than the number of bonds formed by bulk SPC/E water.

In Fig. 6 we have investigated the distribution of hydrogen bond angles in more detail for the simulation with $`N_W=2600`$. The figure depicts the distribution of $`\mathrm{cos}\varphi `$, where $`\varphi `$ is the angle between intramolecular O-H and intermolecular H$`\mathrm{}`$O vector (see inset). Only that part of the distribution which corresponds to almost linear hydrogen bonds is shown. The top frame shows the distribution of water-water hydrogen bond angles close to the center of the cylinder (solid line) and close to the pore surface (dashed). While there is a reduction in the number of water-water bonds, their nature is not changed. By and large, linear water-water hydrogen bonds are favored. There is however, a tendency to have slightly more distorted hydrogen bonds close to the pore surface, as can be inferred from the reduced height of the distribution (dashed line) at $`\mathrm{cos}\varphi =1`$ and the increase for $`\mathrm{cos}\varphi <0.9`$.

The center frame of Fig. 6 contains the corresponding distribution for the hydrogen bonds between water molecules and the bridging oxygen atoms. The distribution is almost identical to that between water molecules. Hence, the loss in water-water bonds is quantitatively (see Fig. 5) and structurally compensated by hydrogen bonds to Vycor atoms. The linearity of hydrogen bonds between water molecules and surface oxygen atoms on hydrophilic substrates was already observed by Lee and Rossky . The angle distribution between water molecules and hydroxyl groups on the surface (bottom frame of Fig. 6) is uniform; there is no tendency to form linear bonds between water and NBO’s as in the case of water and BO’s. The probable cause is the fact that, contrary to Lee and Rossky’s work, all substrate atoms including the hydrogen atoms are immobile during the simulation and that, consequently, the hydroxyl groups are not able to reorientate in an optimal way to form hydrogen bonds with water molecules. Almost linear hydrogen bonds can, however, form between water molecules and BO’s due to the rotational mobility of the water molecule. Since the total number of hydrogen bonds between water and NBO’s is very small, we did not attempt to modify our model of the Vycor surface.

We have further analyzed the local structure around a water molecule. We explored the arrangement of nearest neighbor molecules by means of the distribution of the angle $`\gamma `$, which is the angle between the two vectors from an oxygen atom of a water molecule with the oxygen atoms of the two closest water neighbors (see inset of Fig. 7). Figure 7 shows the angular distribution function $`P(\mathrm{cos}\gamma )`$ in the fully hydrated system with $`N_W=2600`$ for different radial layers (full lines, as indicated). For comparison, the corresponding distribution for bulk water is given as the dashed lines. The arrows point to the average value of $`\mathrm{cos}\gamma `$. Both in experiment and in simulation, the distribution function for bulk water at ambient conditions shows a well defined peak at $`\gamma 105^{}`$ ($`\mathrm{cos}\gamma 0.26`$), representing the tetrahedral order, and a secondary peak at $`\gamma 54^o`$ ($`\mathrm{cos}\gamma 0.6`$), attributed to neighbors located in the cavities of the hydrogen bond network .

In Vycor glass the distributions $`P(\mathrm{cos}\gamma )`$ near the center of the pore are similar to the bulk distribution. In the distribution in the innermost region with density less than ambient density, a more ordered tetrahedral structure than in the bulk phase is quite obvious. Closer to the boundary of the pore, the minimum between the main peak of the tetrahedral coordination and the peak at $`\mathrm{cos}\gamma =0.6`$ begins to fill in and for $`\rho >16`$ Å the distribution is monotonic, indicating a departure from the preferential tetrahedral arrangement. In the layer between 18 and 19 Å, there is an obvious preference for small angles. The average value, $`\mathrm{cos}\gamma `$ (arrows in Fig. 7), shifts to larger values from the center to the surface of the pore. The shift to larger values is a consequence both of the substrate-induced distortion of individual water-water hydrogen bonds and of the higher local density, which leads to a collapse of the tetrahedral structure due to packing effects.

## 5 Summary and Conclusions

We have performed molecular dynamics simulations of water in porous Vycor glass for different degrees of hydration, ranging from 19 to 96 %. In the present work, we analyzed the modifications of the hydrogen bond network of water at the highest degree of hydration. At nearly full hydration the structural properties of water in the center of the pore are similar to the bulk phase. We did not compare the chemical potential of water in the pore to that of bulk water.

Close to the surface the structure changes in several ways. Due to the hydrophilic properties of the surface the density distribution of hydrogen and oxygen atoms shows the formation of two layers of water molecules which interact with the Vycor surface; the water-water hydrogen bond network in the boundary region is distorted, as hydrogen bonds between water molecules are partially substituted by bonds from water to bridging oxygens of the Vycor glass. Hydrogen bonds close to the surface are still preferentially linear but show a larger deviation from linearity than hydrogen bonds in the center of the pore or in bulk water. The local tetrahedral arrangement of water molecules, however, is destroyed close to the surface due to the geometric confinement and the increased local density.

The presence of the pore surface also changes the radial pair distribution functions relative to those of bulk water. The minimum of the oxygen-oxygen radial distribution function fills in with respect to the bulk; this observation is familiar from water under pressure, where it also indicates a decrease in hydrogen bond order. The comparison of the pair correlation functions with experiment is encouraging, and we believe that further studies at different temperatures would be helpful for a deeper understanding of the phenomena concerned with the modification of hydrogen bonding in confined water.

## Acknowledgments

We gratefully acknowledge financial support by the Fonds der Chemischen Industrie and the Clothilde Eberhardt Foundation.


format: text

text:
# 1 On the left, the total (heavy line), potential (dashed line), kinetic (dot-dash line) and magnetic (full line) energies of three CMEs. On the right, their mass (diamonds) and velocity (asterisks).

LASCO and EIT Observations of Coronal Mass Ejections

K. P. Dere

Naval Research Laboratory, Code 7663, Washington DC 20375

A. Vourlidas, P. Subramanian

Center for Earth Observing and Space Research, Institute for Computational Sciences, George Mason University, Fairfax, VA 22030

Proceedings of the YOHKOH 8th Anniversary International Symposium, “Explosive Phenomena in Solar and Space Plasma”, Dec 6-8 1999, Sagamihara, Japan

The LASCO and EIT instruments on the SOHO spacecraft have provided an unprecedented set of observations for studying the physics of coronal mass ejections (CMEs). They provide the ability to view the pre-event corona, the initiation of the CME and its evolution from the surface of the Sun through 30 R. An example of the capability of these instruments is provided in a description of a single event (Dere et al., 1997). During the first 2 years of operation of LASCO and EIT on SOHO, a substantial fraction, on the order of 25 to 50%, of the CMEs observed exhibited structure consistent with the ejection of a helical magnetic flux rope. An examples of these has been reported by Chen et al. (1997) and Dere et al. (1999). These events may be the coronal counterpart of magnetic clouds discussed by Burlaga et al.(1981) and Klein and Burlaga (1982). They analyzed observations of magnetic fields behind interplanetary shocks and deduced that the field topology was that of a helical flux rope.

Recently, we have explored a number of the consequences of the helical flux rope description of these types of CMEs. Vourlidas et al. (1999) examined the energetics of CMEs with data from the LASCO coronagraphs on SOHO. The LASCO observations provide fairly direct measurements of the mass, velocity and dimensions of CMEs. Using these basic measurements, we determined the potential and kinetic energies and their evolution for several CMEs that exhibited a flux-rope morphology. Assuming magnetic flux conservation (’frozen-in’ fields), we used observations of the magnetic flux in a variety of magnetic clouds near the Earth to determine the magnetic flux and magnetic energy in CMEs near the Sun. Figure 1 shows these quantities for a few representative flux rope CMEs. In general, we find that the potential and kinetic energies increase at the expense of the magnetic energy as the CME moves out, keeping the total energy roughly constant. This demonstrates that flux rope CMEs are magnetically driven. Furthermore, since their total energy is constant, the flux rope parts of the CMEs can be considered to be a closed system above $``$ 2 $`R_{}`$.

Subramanian et al. (1999) examined images from LASCO to study the relationship of coronal mass ejections (CMEs) to coronal streamers. We wished to test the suggestion of Low (1996) that CMEs arise from flux ropes embedded in streamers near their base. It is expected that the CME eruption would lead to the disruption of the streamer. To date, this is the most extensive observational study of the relation between CMEs and streamers. The data span a period of 2 years near sunspot minimum through a period of increased activity as sunspot numbers increased. We have used LASCO C2 coronagraph data which records Thomson scattered white light from coronal electrons at heights between 1.5 and 6$`R_s`$. Synoptic maps of the coronal streamers have been constructed from C2 observations at a height of 2.5$`R_s`$ at the east and west limbs. We have superposed the corresponding positions of CMEs observed with the C2 coronagraph onto the synoptic maps. We identified the different kinds of signatures CMEs leave on the streamer structure at this height (2.5$`R_s`$). We find four categories of CMEs with respect to their effect on streamers:

1. CMEs that disrupt the streamer.
2. CMEs that have no effect on the streamer, even though they are related to it.
3. CMEs that create streamer-like structures.
4. CMEs that are latitudinally displaced from the streamer.

Figure 2 summarizes these results. CMEs in categories 3 and 4 are not related to the streamer structure. We therefore conclude that approximately 35% of the observed CMEs bear no relation to the pre-existing streamer, while 46% have no effect on the observed streamer, even though they appear to be related to it.

Previous studies using SMM data (Hundhausen 1993) have made the general statement that CMEs are mostly associated with streamers and that they frequently disrupt it. Our conclusions thus significantly alters the prevalant paradigm about the relationship of CMEs to streamers.

Subramanian and Dere (2000) have examined coronal transients observed on the solar disk in EIT 195 Å images that correspond to coronal mass ejections observed by LASCO during the solar minimum phase of January 1996 through May 1998. The objective of the study is to gain an understanding of the source regions from which the CMEs observed in LASCO images emanate. We compare the CME source regions as discerned from EIT 195 Å images with photospheric magnetograms from the MDI on SOHO and from NSO Kitt Peak, and also with BBSO H$`\alpha `$ images. The overall results of our study suggest that a majority of the CME related transients observed in EIT 195 Å images are associated with active regions. We have carried out detailed case studies of 5 especially well observed events. These case studies suggest that active region CMEs are often associated with the emergence of parasitic polarities into fairly rapidly evolving active regions. CMEs associated with prominence eruptions, on the other hand, are typically associated with long lived active regions. Figure 3 summarizes these results.

References

Chen, J., et al. 1997, ApJ, 338, L194

Dere, K. P., et al., 1997, Solar Phys., 175, 601.

Dere, K. P., Brueckner, G. E., Howard, R. A., Michels, D.J., Delaboudiniere, J.P., 1999, ApJ, 516, 465.

Burlaga, L., Sittler, E., Mariani, F., Schwenn, R., 1981, JGR, 86, 6673

Klein, L. W., Burlaga, L. F., 1982, JGR, 87, 613

Subramanian, P. and Dere, K.P., 2000, ApJ, in preparation

Subramanian, P., Dere, K.P., Rich, N.B., Howard, R.A., 1999, JGR, 104, 22331

Vourlidas, A., Subramanian, P., Dere, K.P., Howard, R.A., 2000, ApJ, in press


text:
# Pulse driven switching in one-dimensional nonlinear photonic band gap materials: a numerical study

## I INTRODUCTION

Nonlinear dielectric materials  exhibiting a bistable response to intense radiation are key elements for an all-optical digital technology. For certain input optical powers there may exist two distinct transmission branches forming a hysteresis loop, which incorporates a history dependence in the system’s response. Exciting applications involve optical switches, logic gates, set-reset fast memory elements etc . Much interest has been given lately to periodic nonlinear structures , in which because of the distributed feedback mechanism, the nonlinear effect is greatly enhanced. In the low intensity limit, these structures are just Bragg reflectors characterized by high transmission bands separated by photonic band gaps . For high intensities and frequencies inside the transmission band, bistability results from the modulation of transmission by an intensity-dependent phase shift. For frequencies inside the gap bistability originates from gap soliton formation , which can lead to much lower switching thresholds .

The response of nonlinear periodic structures illuminated by a constant wave (CW) with a frequency inside the photonic gap is generally separated into three regimes: i) steady state response via stationary gap soliton formation, ii) self-pulsing via excitation of solitary waves, and iii) chaotic. Much theoretical work has been done for systems with a weak sinusoidal refractive index modulation and uniform nonlinearity , or deep modulation multilayered systems , as well as experimental . One case of interest is when the system is illuminated by a CW bias and switching between different transmission states is achieved by means of external pulses. Such switching has already been demonstrated experimentally for various kinds of nonlinearities , but to our knowledge, a detailed study of the dynamics, the optimal pulse parameters and the stability under phase variations during injection, has yet to be performed.

In this paper we use the Finite-Difference-Time-Domain (FDTD)  method to study the time-dependent properties of CW propagation in multilayer structures with a Kerr type nonlinearity. We find our results generally in accord to those obtained for systems with weak linear index modulation , which were solved with approximate methods. We next examine the dynamics of driving the system from one transmission state to the other by injecting a pulse, and try to find the optimal pulse parameters for this switching. We also test how these parameters change for a different initial phase or frequency of the pulse. Finally, we will repeat all work for the case of a linear multilayer structure with a nonlinear impurity layer.

## II FORMULATION

Electromagnetic wave propagation in dielectric media is governed by Maxwell’s equations

$`\mu {\displaystyle \frac{\stackrel{}{H}}{t}}=\stackrel{}{}\times \stackrel{}{E}`$ $`{\displaystyle \frac{\stackrel{}{D}}{t}}=\stackrel{}{}\times \stackrel{}{H}`$ (1)
Assuming here a Kerr type saturable nonlinearity and an isotropic medium, the electric flux density $`\stackrel{}{D}`$ is related to the electric field $`\stackrel{}{E}`$ by

$`\stackrel{}{D}=ϵ_0\stackrel{}{E}+\stackrel{}{P}_L+\stackrel{}{P}_{NL}=ϵ_0\left(ϵ_r+{\displaystyle \frac{\alpha |\stackrel{}{E}|^2}{1+\gamma |\stackrel{}{E}|^2}}\right)\stackrel{}{E}`$ (2)
where $`\gamma 0`$. $`\stackrel{}{P}_L`$ and $`\stackrel{}{P}_{NL}`$ are the induced linear and nonlinear electric polarizations respectively. Here we will assume zero linear dispersion and so a frequency independent $`ϵ_r`$. Inverting this to obtain $`\stackrel{}{E}`$ from $`\stackrel{}{D}`$ involves the solution of a cubic equation in $`|\stackrel{}{E}|`$. For $`\alpha 0`$ there is always only one real root, so there is no ambiguity. For $`\alpha <0`$ this is true only for $`\gamma >0`$. In our study we will use $`\alpha =1`$ and $`\gamma =(ϵ_r1)^1`$ so that for $`|\stackrel{}{E}|\mathrm{}`$, $`\stackrel{}{D}ϵ_0\stackrel{}{E}`$.

The structure we are considering consists of a periodic array of 21 nonlinear dielectric layers in vacuum, each 20 nm wide with $`ϵ_r=3.5`$, separated by a lattice constant $`a=`$200 nm. The linear, or low intensity, transmission coefficient as a function of frequency is shown in Fig. 1a. In the numerical setting each unit cell is divided into 256 grids, half of them defining the highly refractive nonlinear layer. For the midgap frequency, this corresponds to about 316 grids per wavelength in the vacuum area and 1520 grids per effective wavelength in the nonlinear dielectric, where of course the length scale is different in the two regions. Stability considerations only require more than 20 grids per wavelength . Varying the number of the grids used we found our results to be completely converged. On the two sides of the system we apply absorbing boundary conditions .

We first study the structure’s response to an incoming constant plane wave of frequency close to the gap edge $`\omega a/2\pi c=0.407`$. For each value of the amplitude, we wait until the system reaches a steady state and then calculate the corresponding transmission and reflection coefficients. If no steady state is achievable, we approximate them by averaging the energy transmitted and reflected over a certain period of time, always checking that energy conservation is satisfied. Then the incident amplitude is increased to its next value, which is done adiabatically over a time period of 20 wave cycles, and measurements are repeated. This procedure continues until a desired maximum value is reached, and then start decreasing the amplitude, repeating backwards the same routine. The form of the incident CW is

$`E_{CW}(t)=\left(A_{CW}+dA_{CW}{\displaystyle \frac{min\{(tt_0),20T\}}{20T}}\right)e^{i\omega t}`$ (3)
where $`t_0`$ is the time when the amplitude change started, $`A_{CW}`$ is the last amplitude value considered and $`dA_{CW}`$ the amplitude increment. One wave cycle $`T`$ involves about 2000 time steps.

## III RESPONSE TO A CW BIAS

The amplitude of the CW is varied from zero to a maximum of 0.7 with about 40 measurements in between. Results are shown in Fig. 1b, along with the corresponding one from a time-independent approximation. The agreement between the two methods is exact for small intensities, however, after a certain input the output waves are not constant any more but pulsative. This is in accord with the results obtained with the slowly varying envelope approximation for systems with a weak refractive index modulation . It is interesting that the averaged output power is still in agreement with the time-independent results, something not mentioned in earlier work. For higher input values, the solution will again reach a steady state just before going to the second nonlinear jump, after which it will again become pulsative. This time though, the averaged transmitted power is quantitatively different from the one predicted from time-independent calculations.

The nonlinear transmission jump originates from the excitation of a stationary gap soliton when the incident intensity exceeds a certain threshold value. Due to the nonlinear change of the dielectric constant, the photonic gap is shifted locally in the area underneath the soliton, which becomes effectively transparent, resembling a quantum well with the soliton being its bound state solution . The incident radiation coupled to that soliton tunnels through the structure and large transmission is achieved. We obtain a maximum switching time of the order of 100 $`T_r`$, or a frequency of 360 GHz, where $`T_r=2L/c`$ is roundtrip time in vacuum. The second transmission jump is related to the excitation of two gap solitons, which however, are not stable and so transmission is pulsative. The Fourier transform of the output shows that, after the second transmission jump, the system pulsates at a frequency $`\omega a/2\pi c=0.407\pm n\times 0.024`$, exactly three times the one of the first pulsating solution $`\omega a/2\pi c=0.407\pm n\times 0.008`$, where $`n`$ is an integer. For much higher input values the response eventually becomes chaotic. A more detailed description of the switching process as well as the soliton generation dynamics can be found in .

## IV PULSE DRIVEN SWITCHING

We next turn to the basic objective of this work. We assume a specific constant input amplitude $`|A_{CW}|=0.185`$ corresponding to the middle of the first bistable loop. Depending on the system’s history, it can be either in the low transmission state I, shown in Fig. 1c, or in the high transmission state II, shown in Fig. 1d, which are both steady states. We want to study the dynamics of a pulse injected into a system like that. More specifically, if it will drive the system to switch from one state to the other, how the fields change in the structure during switching, for which pulse parameters this will happen and if these parameters change for small phase and frequency fluctuations. We assume Gaussian envelope pulses

$`E_P(0,t)=A_Pe^{(tt_05t_w)^2/t_w^2}e^{i\omega t}`$ (4)
where $`A_P`$ is the pulse amplitude, $`t_0`$ the time when injection starts, and $`W=2t_wc`$ is the pulse’s full width at $`1/e`$ of maximum amplitude. The beginning of time $`t`$ is the same as for the CW, so there is no phase difference between them. After injection we wait until the system reaches a steady state again and then measure the transmission and reflection coefficients to determine the final state. During this time we save the field values inside the structure every few time steps, as well as the transmitted and reflected waves. This procedure is repeated for various values of $`A_P`$ and $`t_w`$, for both possible initial states. Our results are summarized in Fig. 2. White areas indicate the pulse parameters for which the intended switch was successful while black are for which it failed. In Fig. 2a, or the “Switch” graph, the intended switching scheme is for the same pulse to be able to drive the system from state I to state II and vise versa. Fig. 2b, or “Switch All Up”, is for a pulse able to drive the system from I to II, but fails to do the opposite, ie. the final state is always II independently from which the initial state was. Similarly, Fig. 2c or “Switch All Down”, is for the pulse whose final state is always I, and Fig. 2d or ”No Switch”, for the pulse that does not induce any switch for any initial state.

We find a rich structure on these parameter planes. Note also that there is a specific cyclic order as one crosses the curves moving to higher pulse energies: $``$ d $``$ b $``$ a $``$ c $``$ etc. This indicates that there must be some kind of energy requirements for each desired switching scheme. After analyzing the curves it was found that only the first one in the “Switch” graph could be assigned to a simple constant energy curve $``$$`W|A_p|^2`$. Since any switching involves the creation or destraction of a stationary soliton, then this should be its energy. In order to put some numbers, if we would assume a nonlinearity $`|\alpha |=10^9`$ cm2/W, then we would need a CW of energy $`34`$ MW/cm2 and a pulse of width $`W/c`$ of a few tenths of femtoseconds and energy $``$$`2.5\mu `$J/cm2. These energies may seem large, but they can be sufficiently lowered by increasing the number of layers and using an incident frequency closer the the gap edge.

In order to find more about how the switching occurs, we plotted in Fig. 3 the effective transparent areas and the output fields as a function of time, for the first three curves of Fig. 2a. As expected, for the pulse from the first curve, the energy for the soliton excitation is just right, and the output fields are small compared to the input. For the other curves however, there is an excess of energy. The system has to radiate this energy away before a stable gap soliton can be created. It is interesting to note that this energy goes only in the transmitted wave, not the reflected, and it consists of a series of pulses . For the second curve in Fig. 2a there is one pulse, for the third there are two etc. The width and frequency of the pulses are independent of the incident pulse, they are the known pulsating solutions we found in the CW case. So the system temporally goes into a pulsating state to radiate away the energy excess before settling down into a stable state. If this energy excess is approximately equal to an integer number of pulses (the solitary waves from the unstable solutions), then we will have a successful switch, otherwise it will fail. A similar behavior is found in the system’s response during switch down for the first three curves in Fig. 2a, using exactly the same pulses as before for the switch up. So the same pulse is capable of switching the system up, and if reused, switching the system back down. Using the numbers assumed before for the nonlinearity $`\alpha `$, the pulses used in Fig. 3 are (a) $`W/c=14`$ fs, $``$=2.5 $`\mu `$J/cm2, (b) $`W/c=28`$ fs, $``$=12 $`\mu `$J/cm2, (c) $`W/c=42`$ fs, $``$=32 $`\mu `$J/cm2.

Up to now, the injected pulse has been treated only as an amplitude modulation of the CW source, ie. they had the same exactly frequency and there was no phase difference between them. The naturally rising question is how an initial random phase between the CW and the pulse, or a slightly different frequency affect our results. We repeated the simulations for various values of an initial phase difference, first keeping them with the same frequency. We find that although the results show qualitatively the same stripped structure as in Fig. 2, there are quantitative differences. The main result is that there is not a set of pulse parameters that would perform the desired switching successfully for any initial phase difference. Thus the pulse can not be incoherent with the CW, ie. generated at different sources, if a controlled and reproducible switching mechanism is desired, but rather it should be introduced as an amplitude modulation of the CW. However, if this phase could be controlled, then the switching operation would be controlled, and a single pulse would be able to perform all different operations.

The picture does not change if we use pulses of slightly different frequency from the source. We used various pulses with frequencies both higher and lower than the CW, and we found a sensitive, rather chaotic, dependence on the initial phase at injection time. The origin of this complex response, if it is an artifact of the simple Kerr-type nonlinearity model that we used, and if it should appear for other kinds of nonlinearities, is not yet clear to as. More work is also needed on how these results would change if one used a different $`|A_{CW}|`$ not in the middle of the bistable loop, a wider or narrower bistable loop etc., but these would go more into the scope of engineering.

## V LINEAR LATTICE WITH A NONLINEAR IMPURITY LAYER

Besides increasing the number of layers to achieve lower switching thresholds, one can use a periodic array of linear layers $`ϵ=ϵ_0ϵ_r`$ with a nonlinear impurity layer  $`ϵ=ϵ_0(ϵ_r^{}+\alpha |\stackrel{}{E}|^2)`$ where $`ϵ_rϵ_r^{}`$ and we will use $`\alpha =+1`$ and $`\gamma =0`$. This system is effectively a Fabry-Perot cavity with the impurity (cavity) mode inside the photonic gap, as shown in Fig. 4a. The bistable response originates from its nonlinear modulation with light intensity. The deeper this mode is in the gap, the stronger the linear dispersion for frequencies close to it. Because of the high $`Q`$ of the mode, we can use frequencies extremely close to it achieving very low switching thresholds . Here however we only want to study the switching mechanism, so we will use a shallow impurity mode.

The bistable input-output diagram, the output fields during switching and the field distributions in the two transmission branches are also shown in Fig. 4. We observe a smaller relaxation time and of course the absence of pulsating solutions. The parameters used are $`ϵ_r^{}=1`$ and $`\omega a/2\pi c=0.407`$ which corresponds to a frequency between the mode and the gap edge. We want to test if a pulse can drive this system in switching between the two different transmission states, and again test our results against phase and frequency perturbations. The two states shown in Fig. 4 are for an input CW amplitude of $`|A_{CW}|=0.16`$. The results for coherent, pulse and CW, are shown in Fig.5. We see that any desired form of switching can still be achieved, but the parameter plane graphs no more bare any simple explanations as the ones obtained for the nonlinear superlattice. Repeating the simulations for incoherent beams and different frequencies we obtain the same exactly results as before. Only phase-locked beams can produce controlled and reproducible switching.

## VI CONCLUSIONS

We have studied the time-dependent switching properties of nonlinear dielectric multilayer systems for frequencies inside the photonic band gap of the corresponding linear structure. The system’s response is characterized by both stable and self-pulsing solutions. We examined the dynamics of driving the system between different transmission states by pulse injection, and found correlations between the pulse, the stationary gap soliton and the unstable solitary waves. A small dependence on the phase difference between the pulse and the CW is also found, requiring coherent beams for fully controlled and reproducible switching. Similar results are also found for the case of a linear periodic structure with a nonlinear impurity.

###### Acknowledgements.

Ames Laboratory is operated for the U. S. Department of Energy by Iowa State University under contract No. W-7405-ENG-82. This work was supported by the Director of Energy Research office of Basic Energy Science and Advanced Energy Projects, the Army Research office, and a PENED grand.


text:
# Silicate emission in Orion Based 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 Orion nebula is one of the most studied star–forming regions in the Galaxy. The ionizing stars of the Orion nebula (the Trapezium stars, the hottest of which is $`\theta ^1`$ Ori C, O6) have eroded a bowl–shaped H ii region into the surface of the Orion molecular cloud. The Orion bar is the limb–brightened edge of this bowl where an ionization front is progressing into the molecular cloud. It is seen as an elongated structure at a position angle of approximately 60°. The Orion nebula extends to the North. The Trapezium stars are located at an angular distance of approximately 2.3 arc minutes from the bar, corresponding to 0.35 pc at a distance of 500 pc. The molecular cloud extends to the other side of the bar, but also to the back of the Orion nebula. The bright star $`\theta ^2`$ Ori A (O9.5Vpe) lies near the bar, and is clearly in front of the molecular cloud since its color excess is only E(B–V) $``$ 0.2 mag.

Figs. 1 and 2 illustrate the geometry of the region observed. Fig. 1 shows six representative images of the region of the Orion bar (see the figure caption for details). Fig. 2 shows the contours of the the \[Ne iii\] 15.5$`\mu `$m fine-structure line emission which delineates the H ii region. The emission in one of the mid-IR bands (hereafter called the Aromatic Infrared Bands, AIBs) at 6.2$`\mu `$m traces the Orion bar (an edge-on PhotoDissociation Region or PDR). The AIBs are usually strongly emitted by PDRs. The Trapezium region was avoided because of possible detector saturation.

Pioneering infrared observations by Stein & Gillet (Stein (1969)) and Ney et al. (Ney (1973)) discovered interstellar silicate emission near 10$`\mu `$m in the direction of the Trapezium. This was confirmed by Becklin et al. (Becklin (1976)) who also noticed extended silicate emission around $`\theta ^2`$ Ori A. Since that time, interstellar silicate emission has been found by the Infrared Space Observatory (ISO) in the H ii region N 66 of the Small Magellanic Cloud (Contursi et al. in preparation) and in a few Galactic compact H ii regions (Cox et al. in preparation) and Photodissociation Regions (PDRs, Jones et al. in preparation). The emission consists of two broad bands centered at 9.7 and 18$`\mu `$m, which show little structure and are clearly dominated by amorphous silicates.

We report in the present article ISO observations of the Orion bar and of a part of the Orion nebula made with the Circular Variable Filter of the ISO camera (CAM-CVF) which allowed imaging spectrophotometry of a field $`3\mathrm{}\times 3\mathrm{}`$ at low wavelength resolution (R $``$ 40). We also use an ISO Short–Wavelength Spectrometer (SWS) observation which provides higher–resolution ($`R1000`$) spectroscopy at a position within the H ii region (see Fig. 2). This spectrum was taken as part of the MPEWARM guaranteed-time program. We show here that these new ISO data confirm and extend previous observations of the amorphous silicate emission and also give evidence for emission by crystalline silicates.

Sect. 2 of this paper describes the observations and data reduction. In Sect. 3, we discuss the emission of dust and gas. The silicate emission is characterized in section 4 through modelling of the observed continuum IR emission. Finally, conclusions are presented in Sect. 5. Our observations also give information on the fine–structure lines and on the AIBs. This will be presented in Appendices A and B respectively.

## 2 Observations and data reduction

Imaging spectrophotometry was performed with the 32$`\times `$32 element mid-IR camera (CAM) on board the ISO satellite, using the Circular Variable Filters (CVFs) (see Cesarsky et al. 1996a  for a complete description). The observations employed the 6″ per pixel field-of-view of CAM. Full scans of the two CVFs in the long-wave channel of the camera were performed with both increasing and decreasing wavelength. The results of these two scans are almost identical, showing that the transient response of the detector was only a minor problem for these observations. The total wavelength range covered is 5.15 to 16.5$`\mu `$m and the wavelength resolution $`\lambda /\mathrm{\Delta }\lambda 40`$. 10$`\times `$0.28 s exposures were added for each step of the CVF, and 7 more at the first step in order to limit the effect of the transient response of the detectors. The total observing time was about 1 hour. The raw data were processed as described in Cesarsky et al. (1996b ), with improvements described by Starck et al. (Starck (1998)) using the CIA software111CIA is a joint development by the ESA Astrophysics Division and the ISOCAM Consortium led by the ISOCAM PI, C. Cesarsky, Direction des Sciences de la Matière, C.E.A., France. The new transient correction described by Coulais & Abergel (Coulais (1998)) has been applied but the corrections introduced are minimal, as mentioned above. The bright star $`\theta ^2`$ Ori A is visible in the maps of several spectral components and has been used to re–position the data cube. This involved a shift of only 2″. The final positions are likely to be good to 3″ (half a pixel).

All the maps presented here were obtained from the CVF data cube and have approximately the same resolution: namely, 6″ pixels at the short wavelengths increasing to about 8″ at 15$`\mu `$m; see Appendix C for more details.

In several of these maps a faint emission can be seen on the south–east part of the ISOCAM field of view. This feature does not correspond to anything conspicuous in published images of the region, in particular in the near–IR images of Marconi et al. (Marconi (1998)). It is a spurious feature due to multiple reflections of the strong Trapezium between the detector and the CVF filter wheel, as shown by the ISOCAM ray tracing studies of Okumura (koko (1999)).

The complete SWS scan (2.4-46$`\mu `$m) was reduced with the latest version of SWS-IA running at the Institut d’Astrophysique Spatiale. Calibration files version CAL-030 were used.

Fig. 3 presents the SWS spectrum ($`\lambda /\mathrm{\Delta }\lambda 1000`$) obtained inside the H ii region at the position indicated on Fig. 2. Fig. 3 also shows the comparison of the SWS spectrum with that of the CAM-CVF pixels averaged in the SWS aperture. The agreement between these spectra is excellent, well within 20 percent for the continuum.

## 3 Gas and dust emission

The spectra towards the H ii region of the Orion nebula are shown in Figs. 3 and 4. The CAM-CVF spectrum of Fig. 3 is representative of the whole field because CVF spectra obtained at different positions in the H ii region and around $`\theta ^2`$ Ori A look qualitatively similar (compare Figs. 3 and 10 which show the CVF spectra of different pixels; note particularly the rising long wavelength portion of the spectra).

In the SWS spectrum, a large number of unresolved lines from atoms, ions and molecules are visible. We note the Pf$`\alpha `$ recombination line of hydrogen (emitted by the warm, ionized gas of the H ii region) and the molecular hydrogen pure rotation lines S(2) and very faintly S(3) and S(5) (stemming from the cooler, molecular PDR gas). The simultaneous presence of these lines reflects the variety of physical conditions present along the line of sight. Clearly, we are looking at emission from the H ii region mixed with some emission from the background PDR. These unresolved lines are briefly discussed in Appendix A.

The other striking fact of the SWS spectrum is the strong continuum peaking at about 25$`\mu `$m. It is emitted by warm dust in the H ii region, but dust from the background PDR probably also contributes. The broad emission bands of amorphous silicates centered at 10 and 18$`\mu `$m are visible. The classical AIBs at 6.2, 7.7, 8.6, 11.3 and 12.7$`\mu `$m dominate the mid-IR part of the spectrum. As discussed by Boulanger et al. bbcr (1998), the mid-IR spectrum can be decomposed into Lorentz profiles (the AIBs) and an underlying polynomial continuum. Maps of the various AIBs constructed in this way all show the same morphology originating mainly from the PDR gas in the Orion bar (see Appendix B). We will hereafter use the 6.2$`\mu `$m-band as representative of the behaviour of the AIBs.

In Fig. 5 we compare the behaviour of the mid-IR continuum emission and of the AIBs. Clearly, the AIB emission is concentrated in the Orion bar whereas the 15.5$`\mu `$m-continuum emission extends throughout the whole CAM field and shows a local peak around $`\theta ^2`$ Ori A (note that the mid-IR emission around this star is foreground because the star lies in front of the nebula). The continuum emission, however, appears to peak towards $`\theta ^1`$ Ori C, outside the region observed with ISOCAM.

The contrast in the emission morphology between the bands and continuum can be interpreted in terms of the photodestruction of the AIB carriers in the hard UV-radiation field of the H ii region. The AIB carriers must be efficiently destroyed while the larger grains are much more resistant (e.g. Allain et al. 1996). We detail the modelling of the dust thermal emission in the next section.

### 3.1 Modelling the dust emission

To account for the observed SWS spectra, we have calculated the thermal equilibrium temperature of dust in the Orion H ii region as a function of distance of the Orion bar from the Trapezium stars, assuming that $`\theta ^1`$ Ori C (an O6 star) dominates the local radiation field. We use the optical constants of the amorphous astronomical silicate of Draine (Draine85 (1985)) and of the amorphous carbon AC1 of Rouleau & Martin (Rouleau (1991)). Assuming typical interstellar grain sizes (e.g. Draine & Lee DandL84 (1984)), we find a temperature range of 85–145 K for amorphous silicates and a range of 110–200 K for amorphous carbon, corresponding to grains of radius 1500 and 100 Å respectively, at a distance of $`0.35`$ pc from $`\theta ^1`$ Ori C (the distance of the Orion Bar to the Trapezium stars).

Using, for simplicity, discrete dust temperatures consistent with those calculated above (Tsilicate = 80 K and 130 K, Tcarbon = 85 K and 155 K) we are able to satisfactorily model the continuum emission spectrum from the dust in the Orion H ii region at the position of the ISO-SWS spectrum. In Fig. 4 we show the calculated emission spectrum from our model where we adopt the carbon/silicate dust mass ratios of Draine & Lee (DandL84 (1984)). In the calculated spectrum we have included the emission from carbon grains at 300 K, containing $``$ 1 percent. of the total carbon dust mass, in order to fit the short wavelength continuum emission. The hot carbon grain emission mimics that of the stochastically-heated Very Small Grains (VSGs, Désert et al. 1990). The 300 K temperature represents a mean of the temperature fluctuations for these small particles in the radiation field of $`\theta ^1`$ Ori C, and therefore indicates a lower mass limit of $``$ 1 percent for the mass of the available carbon in VSGs.

The results of our model show that the emission feature in the 10$`\mu `$m region is dominated by amorphous silicates at temperatures of the order of 130 K, but that there may also be a small contribution from amorphous carbon grains in the 12$`\mu `$m region (Fig. 4). We also note broad “features” in the SWS spectrum, above the modelled continuum in Fig. 4, at $`1520`$$`\mu `$m, $`2028`$$`\mu `$m and longward of 32$`\mu `$m, that are not explained by our model. These features bear a resemblance to the major bands at 19.5, 23.7 and 33.6$`\mu `$m seen in the crystalline forsterite spectra of Koike et al. (Koike (1993)) and of Jaeger et al. (Jaeger (1998)). Bands in these same wavelength regions were noted by Jones et al. (Jones (1998)) in the SWS spectra of the M 17 H ii region and were linked with the possible existence of crystalline Mg-rich olivines in this object. Thus, similar broad emission bands are now observed in the 15–40$`\mu `$m wavelength region of the SWS spectra of two H ii regions (Orion and M 17). These bands resemble those of the crystalline Mg-rich silicate forsterite. Another band at 9.6$`\mu `$m is probably due to some sort of crystalline silicate, and will be discussed in more details in the next section.

This dust model is simple–minded but emphasizes dust spectral signatures in the mid-IR continuum which was the main aim here. More detailed modelling treating temperature fluctuations and taking into account the grain size distribution is underway (Jones et al. in preparation).

The broad continua that lie above the model fit (i.e. $`2028`$$`\mu `$m and $`>`$ 32$`\mu `$m, Fig. 4) can be associated with crystalline silicate emission bands. This seems to be a robust conclusion of this study. The features are too narrow to be explained by single-temperature blackbody emission and are therefore likely to be due to blended emission features from different materials. Unfortunately, having only one full SWS spectrum and CVF spectra that do not extend beyong 18$`\mu `$m, we are unable to say anything about the spatial variation of these broad bands in the Orion region.

Interestingly, broad plateaux in the $`1520`$$`\mu `$m region have been associated with large aromatic hydrocarbon species containing of the order of a thousand carbon atoms (van Kerckhoven et al. vanKerckhoven (2000)). However, in this study the integrated intensity of the $`1520`$$`\mu `$m plateaux do vary by a factor of up to 10 relative to the aromatic carbon features shortward of 13$`\mu `$m. Thus, the origin of these broad emission features does remain something of an open question at this time.

## 4 Tracing the silicate emission

To delineate the spatial extent of the 10$`\mu `$m-silicate emission conspicuously visible in Fig. 3 and 4, we proceed as follows. We start with the spectrum towards $`\theta ^2`$ Ori A, which shows the most conspicuous silicate emission and we represent the AIBs by Lorentz profiles, see Fig. 6 (top). Next we subtract them from the CVF spectra. The remaining continuum has the generic shape of a blackbody on top of which we see the broad bands corresponding to the silicate emission, Fig. 6 (middle). Finally, we subtract a second order polynomial from the continuum thus obtaining the well known silicate emission profile at that position, Fig. 6 (bottom). The profile thus obtained is then used as a scalable template to estimate the emission elsewhere, see Appendix C for more details.

On top of the broad band of amorphous silicate centered near 9.7$`\mu `$m we see a band centered at nearly 9.6$`\mu `$m, which we ascribe to crystalline silicates (Jaeger et al. Jaeger (1998)). This band was also used as a scalable template as explained above and in Appendix C. Finally, the S(5) rotation line of H2 at 6.91$`\mu `$m is present and is probably blended with the \[Ar ii\] line at 6.99$`\mu `$m.

In Fig. 7, we see that the spatial distribution of the 9.7$`\mu `$m-feature of amorphous silicate is quite similar to that of the 15.5$`\mu `$m-continuum. The 15.5$`\mu `$m continuum emission includes a strong contribution from silicates (see Fig. 4), but a peak in the silicate emission around $`\theta ^2`$ Ori A is also evident. The silicate emission is thus predominantly due to larger grains. The narrower 9.6$`\mu `$m feature is mapped in Fig. 8. We note its similarity to the distribution of the 9.7$`\mu `$m broad band: this fact lends support to our assignation of this band to crystalline silicate.

Due to the low spectral resolution of the CAM-CVF, however, the 9.6$`\mu `$m feature will certainly blend with the S(3) pure rotational line of molecular hydrogen - if present. To check this we have compared our 9.6$`\mu `$m map to that of molecular hydrogen in its fluorescent vibrational line 1$`0`$ S(1) (2.12$`\mu `$m). Courtesy of P.P. van der Werf (van der Werf et al. vanderwerf (1996)), we reproduce in Fig. 9 the map of the fluorescent molecular hydrogen emission. This latter correlates better with the AIB emission as traced by the 6.2$`\mu `$m-feature (bottom figure) than it does with the tentative crystalline silicate emission (top), namely they both peak along the bar. This is not surprising because the H2 and AIB emitters require shielding from far-UV radiation to survive. Conversely, the 9.6$`\mu `$m silicate feature is stronger where H2 is weak as can be seen around $`\theta ^2`$ Ori A. In addition, the H2 S(3) rotational line at 9.66$`\mu `$m is detected in the ISO-SWS spectrum of the Orion bar presented in Verstraete et al. (1999, in preparation) with an intensity of $`6\times 10^7`$ W m-2 sr-1. This value is a factor of 16 below the median flux of the 9.6$`\mu `$m feature in our map, namely $`10^5`$ W m-2 sr-1. We can thus safely conclude that our 9.6$`\mu `$m-emission predominantly originate from silicates. A confirmation of the identification of the 9.6$`\mu `$m band with a crystalline silicate dust component would be possible if a second signature band were seen in our spectra. The SWS spectrum (Fig. 4) shows only broad emission bands that are difficult to characterise, and additionally, the chacteristic crystalline olivine band in the $`11.211.4`$$`\mu `$m region (e.g. Jaeger et al. Jaeger (1998)), if present, is blended with the 11.2$`\mu `$m aromatic hydrocarbon feature. Additionally, most of the chacteristic crystalline bands fall longward of the CVF spectra. Thus, it is difficult to self-consistently confirm the 9.6$`\mu `$m band identification with the presented data.

In summary, emission in the 9.7$`\mu `$m band of amorphous silicate emission exists everywhere inside the Orion H ii region. Previously, amorphous silicate emission had only been seen in the direction of the Trapezium (Stein & Gillett Stein (1969); Forrest et al. Forrest (1975); Gehrz et al. Gehrz (1975)). We may assume that the 18$`\mu `$m band is also widely present in the region, as witnessed by the single SWS spectrum (Fig. 4) and by the generally rising long wavelength end of ISOCAM spectra; the two spectra shown, Figs. 3 and 10 are quite representative of the steeply rising continuum longward of 15$`\mu `$m.

### 4.1  The interstellar silicate and H2 emission around $`\theta ^2`$ Ori A

The case of $`\theta ^2`$ Ori A is particularly interesting because the geometry is simple and therefore allows quantitative calculations. Moreover, the thermal radio continuum, the recombination lines and the fine–structure lines are faint in the neighbourhood of this star (Felli et al. Felli (1993); Pogge et al. Pogge (1992); Marconi et al. Marconi (1998), and the present paper, Fig. 2). $`\theta ^2`$ Ori A is classified as an O9.5Vpe star and shows emission lines (see e.g. Weaver & Torres–Dodgen Weaver (1997)). It is a spectroscopic binary and an X-ray source. There is little gas left around the star and the observed silicate dust (Fig. 8) is almost all that is visible of the interstellar material left over after its formation. Indeed, O stars are not known to produce dust in their winds which are probably much too hot, so that the silicates we see here must be of interstellar origin.

The mid–IR continuum observed towards $`\theta ^2`$ Ori A can be accounted for by combining emission of warm silicate and carbon grains (see Fig. 10). The model continuum was obtained in the same way as for the SWS observation (see Fig. 4) and with the same assumptions. The grain temperatures are consistent with the heating of interstellar grains by the strong radiation field of the star.

As discussed above, the band near 9.6 $`\mu `$m (Fig. 6 bottom and Fig. 8) may be due to crystalline silicates, any contribution of the S(3) H2 line to this band is minor. Another band at 14 $`\mu `$m (see Figs. 4 and 10) might also be due to crystalline silicates. Amongst the crystalline silicates whose mid–IR absorption spectra are shown by Jaeger et al. (Jaeger (1998)), synthetic enstatite (a form of pyroxene) might perhaps match the $`\theta ^2`$ Ori A spectrum. The interest in the possible presence of crystalline silicates around this star is that they would almost certainly be of interstellar origin, pre–dating the formation of the star. Observations at longer wavelengths are needed for a definitive check of the existence of crystalline silicates and for confirming their nature. Such observations do not exist in the ISO archives and should be obtained by a future space telescope facility.

## 5 Conclusions

We obtained a rather complete view of the infrared emission of the Orion nebula and its interface with the adjacent molecular cloud. The most interesting results are the observation of amorphous, and possibly crystalline, silicates in emission over the entire H ii region and in an extended region around the bright O9.5Vpe star $`\theta ^2`$ Ori A. We have fitted the mid–IR continuum of the H ii region and around $`\theta ^2`$ Ori A with the emission from amorphous silicate and amorphous carbon grains at the equilibrium temperatures predicted for the grains in the given radiation field. This shows that both types of grains can survive in the harsh conditions of the H ii region. A number of bands (the 9.6$`\mu `$m bump seen in Fig. 6; the excess 14$`\mu `$m emission indicated in Figs. 4 and 10) suggest emission from crystalline silicates (essentially forsterite) in the H ii region. Crystalline silicates may also exist around $`\theta ^2`$ Ori A, but further, longer wavelength observations are required to confirm their presence.

Do the observed crystalline silicates result from processing of amorphous silicates in the H ii region or in the environment of $`\theta ^2`$ Ori A? Silicate annealing into a crystalline form requires temperatures of the order of 1000 K for extended periods (Hallenbeck et al. Hallenbeck (1998)). The dust temperatures observed in the H ii region and around $`\theta ^2`$ Ori A are considerably lower than this annealing temperature. One might however invoke grain heating following grain–grain collisions in the shock waves that are likely to be present in the H ii region. However, grain fragmentation rather than melting is the more likely outcome of such collisions (Jones et al. JTH (1996)). It is probable that the crystalline silicates observed here were already present in the parent molecular cloud, and probably originate from oxygen–rich red giants.

Emission by both amorphous and crystalline silicates has been observed with ISO around evolved stars (Waters et al. Waters (1996); Voors et al. Voors (1998)). The crystalline silicates there must have been produced locally by annealing of amorphous silicates. Gail & Sedlmayr (Gail (1999)) have shown that this is possible, and that both amorphous and crystalline forms can be released into the interstellar medium. However, there is no evidence for absorption by crystalline silicates in the general interstellar medium in front of the deeply embedded objects for which amorphous silicate absorption is very strong (Demyk et al. 1999, Dartois et al. Dartois (1998)). Consequently, crystalline silicates represent only a minor fraction compared to amorphous silicates. It would be difficult to detect the emission from a small crystalline component of dust in the diffuse interstellar medium because the dust is too cool (T$`20`$K) to emit strongly in the $`1540\mu `$m wavelength region. Observations of H ii regions and bright stars provide the opportunity of observing this emission due to the strong heating of dust. Emission from amorphous and crystalline silicates is seen around young stars (Waelkens et al. Waelkens (1996); Malfait et al. Malfait (1998)) as well as in comets (Crovisier et al. Crovisier (1998)). There are also silicates in meteorites, but their origin is difficult to determine because of secondary processing in the solar system. Crystalline silicates in comets, and perhaps in interplanetary dust particles believed to come from comets (Bradley et al. Bradley (1992)), must be interstellar since the material in comets never reached high temperatures. However, the silicates probably experienced changes during their time in the interstellar medium. It is interesting to note that while very small grains of carbonaceous material exist, there seem to be no very small silicate grains in the interstellar medium (Désert et al. 1986).

###### Acknowledgements.

A.P. Jones is grateful to the Société de Secours des Amis des Sciences for funding during the course of this work. We are grateful to P.P. van der Werf for providing us with his map of molecular hydrogen emission.

Appendix A: the mid–IR line emission from the Orion nebula

We presented in Fig. 2 a map of the studied region in the \[Ne iii\] line at 15.5 $`\mu `$m. Fig. 11 displays the map of the \[Ar ii\] line at 7.0 $`\mu `$m superimposed on the map of the \[Ar iii\] 9.0$`\mu `$m line. These maps illustrate the ionization structure of the Orion nebula. The spectral resolution of the CVF does not allow a separation of the the \[Ar ii\] line at 6.99 $`\mu `$m from the S(5) pure rotation line of H2 at 6.91 $`\mu `$m. However, the bulk of the H2 emission come from deeper in the molecular cloud than that of \[Ar ii\], ie. more to the south-west (see Fig. 9) and the contamination by the S(5) line is probably minor. The SWS spectrum shown here and that taken towards the bar (Verstraete et al. 1999, in preparation) in which the \[Ar ii\] and the H2 S(5) line are well separated from each other, show that the H2 line is a factor 4 or 5 weaker and hence cannot seriously contaminate the \[Ar ii\] map.

The emission by the singly–charged ion \[Ar ii\] is concentrated near the ionization front on the inner side of the bar. This is very similar to what is seen in the visible lines of \[N ii\] $`\lambda `$6578 and \[S ii\] $`\lambda `$6731 (Pogge et al. Pogge (1992) Fig. $`1c`$ and $`1d`$). The detailed correspondence between the maps in these three ions is excellent: note that the optical maps are not much affected by extinction. The ionization potentials for the formation of these ions are 15.8, 14.5 and 10.4 eV for Ar ii, N ii , S ii respectively, and are thus not too different from each other.

The emission from the doubly–charged ions \[Ne iii\] and \[Ar iii\] shows a very different spatial distribution, with little concentration near the bar but increasing towards the Trapezium. The \[Ne iii\] map (Fig. 2) is very similar to the \[O iii\]$`\lambda `$5007 line map (Pogge et al. Pogge (1992) Fig. $`1e`$), as expected from the similarity of the ionization potentials of \[Ne ii\] and \[O ii\], respectively 41.1 and 35.1 eV. However, the distribution of the \[Ar iii\] line (Fig. 11) is somewhat different, with a trough where the \[Ne iii\] and the \[O iii\] lines exhibit maxima. Ar iii is ionized to Ar iv at 40.9 eV, almost the same ionization potential as that of Ne ii, so that Ar iv (not observable) should co–exist with Ne iii and Ar iii with Ne ii. A map (not displayed) in the 12.7 $`\mu `$m feature, which is a blend of the 12.7 $`\mu `$m AIB and of the \[Ne ii\] line at 12.8 $`\mu `$m, is indeed qualitatively similar to the \[Ar iii\] line map in the H ii region. It differs in this region from the maps in the other AIBs, showing that it is dominated by the \[Ne ii\] line.

As expected, the dereddened distribution of the H$`\alpha `$ line (Pogge et al. Pogge (1992) Fig. $`3b`$), an indicator of density, is intermediate between that of the singly–ionized and doubly–ionized lines.

Appendix B: the AIB emission

Maps of the emission of the 6.2 and 11.3$`\mu `$m AIBs are shown in Fig. 12. We do not display the distribution of the other AIBs because they are very similar. All the spectra of Figs. 3, 4 and 6 show the classical UIBs at 6.2, 7.7, 8.6, 11.3 and 12.7 $`\mu `$m (in the CAM-CVF data the latter is blended with the \[Ne ii\] line at 12.8 $`\mu `$m). There are fainter bands at 5.2, 5.6, 11.0,, 13.5 and 14.2$`\mu `$mvisible in the SWS spectrum of Fig. 3: they may be AIBs as well. All the main bands visible in the CVF spectra are strongly concentrated near the bar. Emission is observed everywhere, because of the extension of the PDR behind the Orion nebula and the presence of fainter interfaces to the South–East of the bar. We confirm the general similarity between the distributions of the different AIBs through the Orion bar observed by Bregman et al. (Bregman (1989)).

We thus conclude that, although the excitation conditions vary greatly from the Trapezium region towards the South–West of the bar, the mixing of fore– and background material along the line of sight does not allow us to observe spectroscopical changes in the AIB emission features (due e.g. to ionization or dehydrogenation as in M17-SW, Verstraete et al. verstraete (1996)).

Appendix C: Estimates of emission strengths

Spectral emission maps have been obtained using one or another of three different methods. The emission from well defined and rather narrow spectral features, viz. AIBs and ions, can be estimated either by numerical integration of the energy within the line and an ad–hoc baseline (method 1), or by simultaneous fit of Lorentz (Boulanger et al. , bbcr (1998)) and/or gauss profiles, including a baseline, determined by a least square fitting algorithm (method 2). The strength of features not amenable to an analytical expression, like the suspected amorphous silicate emission (see Fig. 6) has been estimated using the following method (method 3). We have constructed an emission template consisting of all the observed emission features, each one arbitrarily normalised to unit peak intensity, see Fig. 13. A least square computer code was then used to obtain, for each of the $`32\times 32`$ lines of sight, a set of multiplying coefficients for each feature present in the template plus a global parabolic baseline so as to minimize the distance between the model and the data points. The number of free parameters is then eleven “line intensities” and three polynomial coefficients, for a total of 14 free parameters to be determined from 130 observed spectral points per line of sight. The main drawback of this method is that it does not allow for varying line widths or line centres; however, given the low resolution of ISOCAM’s CVF this is not a serious drawback. We have found that integrated line emission estimated from methods 2 and 3 give results that agree to within 20 percent; numerical integration of Lorentzian line strengths, on the other hand, badly underestimates the energy carried in the extended line widths and hence this method has not been used.