text:
# Twin model for orthorhombic perovskites

## I Introduction

Transition metal oxides with the perovskite structure exhibit a large variety of interesting physical phenomena, including high $`T_c`$ superconductivity, ferroelectricity, colossal magnetoresistance, and a variety of spin, charge, and orbital orderings. Doping can change the magnetic and electronic properties, which can be reflected in the details of the crystal structure. We are particularly interested in the structural response of phenomena of electronic origin, like Jahn-Teller distortions or orbital order.

Since the metal-oxygen bond lengths and angles determine the exchange interactions, a thorough knowledge of the structure is necessary to understand the physics behind these phenomena. Typically crystal structure information is generated by neutron powder diffraction (NPD). NPD gives direct information about the distances of lattice planes and is very accurate in determining the changes of the lattice parameters with temperature and pressure. Using a full pattern Rietveld Analysis the atomic positions can be determined with an accuracy of $`10^4`$-$`10^3`$.

The alternative is single crystal X-ray diffraction (SXD). SXD is a very powerful experimental application to study the crystal structure of a large variety of materials. In contrast with powder diffraction, which maps all crystal planes to one dimension ($`d`$ spacing, or angle $`2\theta `$), SXD observes all crystal plane reflections separately, even if they are related by symmetry. The relative co-ordinates of the atoms in the unit cell can be determined from the intensity distribution with the same accuracy as NPD. The major difference is that NPD is sensitive to the nucleus, while SXD measures the electron density, which is of interest to understand the electronic properties.

However, for many perovskites SXD is not used because of unconventional twinning of the crystals. Twinning can complicate the structure determination by creating extra reflection spots and/or superimposing reflections on top of each other. Twinning of the perovskites is complex as the cubic parent structure allows not only simple twin axis and planes, but also 3D twin relations. In this paper we will show a method to analyse twinned single crystals. Moreover, this twin model allows a full structure determination including detailed information on distortions of both structural and electronic origin. Besides, we also show that one does not need a larger data set, as is usual in the case of twinning.

Using SXD, we observe a $`2a_p`$x$`2a_p`$x$`2a_p`$ unit cell. It is commonly accepted that the unit cell is $`\sqrt{2}a_p`$x$`2a_p`$x$`\sqrt{2}a_p`$. The observed unit cell originates from a three-dimensional type of twinning that is not restricted to manganites. The model is most likely of general application for a large variety of perovskite $`Pnma`$ crystals. This twinning is unique and it involves a distribution of the $`b`$ axis over three perpendicular cubic axes.

## II The twin structure for manganites

The manganites have generated considerable attention because of the colossal magnetoresistance effect. The role of magnetic order has been widely discussed as the double exchange, inducing the ferromagnetic order, is required to generate a metallic ground state. The role of orbital order is much less understood. While local Jahn-Teller distortions are crucial in explaining the localization of the charge carriers in the paramagnetic state, the long range Jahn-Teller ordering is not well studied, except for undoped LaMnO3. We will show that the destruction of long range orbital order is a second prerequisite for the metallic ground state. In fact, we will show in a separate paper that with increasing hole doping of LaMnO3 not the exceeding of a critical concentration but the suppression of the orbital order generates the metallic ground state. The orbital ordering in perovskites with degenerate $`e_g`$ electrons can be easily measured, whereas for degenerate $`t_{2g}`$ electrons the Jahn-Teller distortions are much smaller. Furthermore, the $`Pnma`$ symmetry can accommodate not only the Jahn-Teller ordering, but also a 3D rotation of the octahedra, known as the GdFeO3 distortion. In this paper we will focus on the main reason why SXD has not been widely used for these perovskites, namely twinning. Twinning in doped LaMnO3 originates from the transition of the highly symmetric cubic parent structure to the orthorhombic symmetry, that accommodates both the GdFeO3 and the Jahn-Teller distortion. In thin films $`ac`$ twinning has been observed. The solution of the twin relations in the crystals allow us to study in detail the ordering of these compounds as influenced by temperature, magnetic state and doping concentration. In this paper we will focus on the twin relations.

At high temperatures most AMnO3 are cubic with a unit cell of $`a=3.9`$ Å. Due to both the small average A-site radius and the JT effect the MnO6 octahedra are rotated and distorted at lower temperatures. We will show that our crystals are twinned. We present here an accurate description of the origin of the twinning, the detwinning process and show that we can still determine the temperature dependence of the structure by measuring only the reflections of the main fraction in one octant of the $`hkl`$ space.

Twinning is reported before in AMnO3 perovskites and related structures. For instance, neutron powder experiments on LaMnO3 showed reflections, at temperatures above the Jahn-Teller transition, that could be indexed on a double cubic unit cell. But full pattern refinement was only possible in the orthorhombic space group $`Pnma`$. Single-crystal electron diffraction showed a double cubic unit cell for SrSnO3, even though CaSnO3 is orthorhombic. Electron microscopy, however, showed the existence of coherent domains, with perpendicular orientations of the doubled $`b`$ axis. The structure was deduced via the O’Keefe-Hyde relations.

The single crystal was mounted on an ENRAF-NONIUS CAD4 single crystal diffractometer. The temperature of the crystal was controlled by, heating, a constant nitrogen flow. Initial measurements were done at 180 K. Temperature dependent measurements were performed on a diffractometer with an adjustable temperature set-up, between 130 K and 300 K.

## III Crystal structure

Most structure research of perovskites focuses on the Mn-O distances and the Mn-O-Mn angles, as these parameters determine the super- and double exchange interactions. Here, we like to stress the importance of a complete structure refinement, including the La-position and the rotation/distortion of the MnO6 octahedra. The basic deformations that determine the deviation from the cubic structure are the GdFeO3 distortion and the Jahn-Teller distortion.

The basic building block of the perovskite structure is a 3.9 Å cube with Mn in the centre and O at the face-centers. The oxygen ions co-ordinate the Mn to form MnO6 octahedra. The A atoms are located at the corners of the cube. The undistorted ’parent’ cubic structure rarely exists, but distorts conventionally to an orthorhombic or rhombohedral symmetry. For an ideal perovskite the ratio between the radii of the A site ion and the transition metal ion is such that the tolerance factor

$$t=\frac{r_{A^{2/3+}}+r_{O^2}}{\sqrt{2}(r_{Mn^{4/3+}}+r_{O^2})}$$
(1)
is equal to one. In the La-Ca system the tolerance factor varies from 0.943 to 0.903 going from LaMnO3 to CaMnO3.

### A GdFeO3 rotation

Due to the small radius of the A-site ion, with respect to its surrounding cage, the MnO6 octahedra tilt and buckle to accommodate the lanthanide. This is known as the GdFeO3 distortion. The cubic state allows one unique oxygen position. Due to the GdFeO3 distortion we need two inequivalent positions to describe the structure. O1 is the in-plane oxygen, on a general position, $`x,y,z`$. Two opposite Mn-O1 bonds have the same length, but the perpendicular bonds need not to be equal. O2 is the apical oxygen, located on a fourfold $`x,\frac{1}{4},z`$ position on the mirror plane. Mn-O2 bonds are always of the same length. Both in the undistorted and the distorted perovskite, the O-Mn-O bond angles are $`180^{}(90^{})`$, but due to the buckling Mn-O-Mn bond angles are no longer $`180^{}`$. A pure GdFeO3 distortion can be obtained with equal Mn-O bond lengths.

### B Jahn-Teller distortion

The Jahn-Teller effect originates from the degenerate Mn $`d^4`$ ion. Two, possible, distortions are associated with the Jahn-Teller effect. Q2 is a orthorhombic distortion, with the in-plane bonds differentiating in a long and a short one. Q3 is the tetragonal distortion with the in-plane bond lengths shortening and the out-of-plane bonds extending, or vice versa. The main result of the JT effect, Q2 and Q3, is that the Mn-O distances become different, which lifts the degeneracy of the $`t_{2g}`$ and $`e_g`$ levels.

### C Glazer’s view on the octahedra

There are three possible rotations for a rigid MnO6 octahedron. One, $`\eta _y`$, changes the atomic positions of the O2 atom only. The others have an effect on the atomic positions of both O1 and O2. These movements will also affect the position of the La atom. Glazer identified all possible sequences of rotations for the perovskite system. Two of these sequences create the symmetry elements that make up the standard perovskites. In Glazer’s notation: $`a^{}b^+a^{}`$ ***) $`p^{}q^+r^0`$ means that along  a rotation of size $`p`$ is alternated positive and negative, along  a rotation of size $`q`$ is always (arbitrarily) positive, and along  there is no rotation.) yields orthorhombic $`Pnma`$ and $`a^{}a^+a^{}`$ will result in rhombohedral $`R\overline{3}c`$. As Glazer worked in a pseudocubic $`2a_p`$x$`2a_p`$x$`2a_p`$ system we have to translate these rotations to the orthorhombic unit cell, which gives $`a^{}b^+0`$. This means that the rotation around the $`c`$ axis, $`\eta _z`$, is zero. This is in agreement with the observation that the four positions of O2 near the Mn at 0,0,0 have the same distance to the $`y=0`$ plane, albeit two are positive and two are negative. By rotating the MnO6 octahedra around the $`x`$ axis, two oxygen are raised from their cubic position and two are lowered. Here we neglect the influence of a rotation $`\eta _y`$ on these coordinates. As a result of $`\eta _z=0`$ the fractional $`x`$ co-ordinate of O1 is also zero. In practice small deviations are found, indicative of the non-rigid behaviour of the MnO6 octahedra. Typically in AMnO3 two rotations and the Q2 JT distortion are observed.

One effect of the GdFeO3 distortion is that the number of formula units per unit cell is enlarged. The unit cell is doubled in the $`b`$ direction, with respect to the original cubic cell. The $`ac`$ plane is also doubled, resulting in $`b2a_p`$, and $`ac\sqrt{2}a_p`$. Despite the rotations and distortion, the structure remains in origin cubic. The shifts of the atoms in the unit cell are small. Furthermore, Mn stays roughly octahedrally surrounded although the Mn-O distances tend to differ and some O-Mn-O angles may be no longer perfectly $`90^{}`$ or $`180^{}`$. The small deviations from the cubic symmetry are reflected in the intensities of the reflections. To illustrate this point, we transform an arbitrary, cubic, reflection $`hkl`$ to the orthorhombic setting by $`h^{}k^{}l^{}`$ = $`h`$+$`l`$,$`2k`$,$`h`$-$`l`$. Thus orthorhombic reflections with $`k^{}`$ even and $`h^{}`$ and $`l^{}`$ both even or both odd stem from planes that already existed in the original cubic unit cell. We observe in diffraction patterns of orthorhombic perovskites that reflections that originate from cubic crystal planes have higher intensities than those that do not.

## IV Twin model

Reflections in reciprocal space are experimentally observed at a regular distance in three orthogonal directions, corresponding to a cubic lattice spacing of $`7.8`$ Å in real space. Although the three orthonormal axes have equal lengths, we could not observe the threefold rotation axis along $`111`$, required for cubic symmetry. However, studying the intensity distribution of planes in $`hkl`$ space with constant $`h`$, $`k`$ or $`l`$ showed much regularity as is shown in Fig. 1. We propose a twinned structure consisting of coherent $`Pnma`$ domains. Due to partial overlap this results in a metric cubic system with $`a2a_p7.8`$ Å.

Twinning is often observed in crystals with a reduced symmetry. For instance, the orthorhombic perovskite LaMnO3 has lower symmetry than cubic SrTiO3. Conventional twin models keep one characteristic axis unchanged and form the domains by rotation around that axis. This is commonly observed in constrained epitaxial thin films and pseudo-two-dimensional crystals like YBa2Cu3O7-δ. Another standard twin is inversion twinning, found in non-centrosymmetric systems. Due to the inversion twin, pseudo-centrosymmetry is obtained or a net polarisation in ferroelectric compounds is reduced to zero. Both twin models have usually the property that either all twin domains have reflections that lie on top of each other, merohedral twins, or there is a non-commensurate set of reflections, e.g. in monoclinic unit cells. Here we propose a more extensive form of twinning, where the different reflections of different domains partially coincide, giving rise to the observed metrically cubic system.

The transformation of cubic to orthorhombic symmetry requires a designation of $`a`$, $`b`$, and $`c`$ with respect to the degenerate cubic axes. There are three possibilities to position the doubled $`b`$ axis along the three original cubic axes. Thus we propose that the three fractions’ $`b`$ axes are oriented along the three original axes of the cubic unit cell, as sketched in Fig. 2. This twin model consists of a $`Pnma`$ unit cell, transformed by rotation along the ’cubic’  axis. We still have the freedom to choose the $`a`$ and $`c`$ axes, perpendicular to the $`b`$ axis, and rotated 45 with respect to the cubic axes. Therefore, this model yields six different orientations of the orthorhombic unit cell. As the differences between $`a/\sqrt{2}`$, $`b/2`$ and $`c/\sqrt{2}`$ are small, we observe the reciprocal superposition of the six orientations as metrically cubic.

The standard refining programs can work with twin models, but only if they consists of merohedral twins. In our case, the reflections do not correspond with an orthorhombic unit cell, but with a, twice as large, cubic unit cell. In appendices A and B, we will elaborately explain the transformation from orthorhombic to double cubic and the way we worked around the refining program.

We first measured a hemisphere in $`hkl`$ space with -$`20<h<20`$, -$`20<k<20`$, and $`0<l<20`$. On this dataset, we could refine our twin model with the six fractions, including the Ca concentration on the A site. Measuring a large range in $`hkl`$ space requires a large amount of time, roughly 10 days in the present set-up. Conventionaly, for crystals with orthorhombic unit cells measuring one octant is sufficient. From refinements on selected parts of the dataset, we concluded that we could investigate the structure at different temperatures by only measuring the positive octant of the main fraction, thereby limiting the measuring time considerably.

## V Discussion

The observed reflections could be indexed on a cubic lattice with $`a7.8`$ Å. Although the three axes had equal lengths, we could not observe the required threefold symmetry axis along $`111`$ for cubic symmetry. This is shown for eeo reflections in Fig. 3. Furthermore, studying the intensity distribution of $`hkl`$ planes with constant $`h`$, $`k`$ or $`l`$ showed much regularity, as is shown in appendix A. Rodr$`\stackrel{´}{\mathrm{i}}`$gruez-Carvajal et al. noted that their neutron powder spectra of pure LaMnO3 could be indexed in a cubic $`2a_p`$x$`2a_p`$x$`2a_p`$ unit cell, above the Jahn-Teller transition temperature. They did not observe a splitting of the peaks. Nevertheless, they could not refine the structure in a cubic unit cell. Note that in powder diffraction spectra there is no direct method to observe the threefold symmetry along the body diagonal. The refinement of the spectra has been done in the conventional $`Pnma`$ setting. They suggested that the powder is twinned but could not give evidence for that, as they only studied powders. The observed reflections, the systematic extinctions and the non-cubic intensity distribution suggest that our single crystals are also twinned.

Typically, Ca-doped LaMnO3 has space group $`Pnma`$. We propose a twinned structure consisting of coherent $`Pnma`$ domains, which results in a metric cubic system with $`a7.8`$ Å. The model is identical to the twinning observed by electron microscopy in CaSnO3. They could however not refine their $`2a_p`$x$`2a_p`$x$`2a_p`$ data on different $`Pnma`$ domains, and deduced the structure from phenomenological relations and the similarities with SrSnO3. In the present work, the different contributions of the twin fractions to the total integrated peak intensity are taken into account. The refined model consisted of the regular parameters in single crystal diffractometer. We considered for each observed reflection the related orthorhombic reflections for the applicable twin fractions and used an identical unit cell. This way, we can refine both the atomic positions of the asymmetric unit and the volume ratio of the twin domains simultaneously.

Several crystals were measured and refined with this detwin model. The refinements showed that for every crystal the distribution of the volume over the twin fractions is different. This is an indication that our crystals are really twinned. Other structural deformations would lead to a constant, not sample dependent, effect on the structure factors and therefore to the same fractions of the twin domains. We can also conclude that the size of the twin domains must be slightly smaller than the magnitude of the measured crystals, i.e. several tenths of a millimeter. Larger twin domains would give rise to crystals of one single domain and smaller domains are more likely to produce a constant spreading of the domains. Non-regular, though constant, distributions could signal a preferential growth direction of the crystal.

### A Intensity distribution

Careful analysis of the intensity distribution showed extra evidence for the proposed models. The full dataset shows that reflections $`hkl`$ of type eee)e denotes an even value for the Miller indices, o stands for odd) are by far the strongest. Reflections having one odd Miller index, oee, eoe, eeo, are second in intensity. eoe has the highest intensity of these reflections as seen in Fig. 3. To understand the intensity distribution in the $`2a_p`$x$`2a_p`$x$`2a_p`$ unit cell, we examine the structural deformations. Despite the effect of Jahn-Teller distortions, GdFeO3 rotations and twinning, the structure remains in origin cubic. This means that the intensity distribution of the peaks with high intensity will always mimic the intensity distribution of the peaks of the undistorted cubic perovskite. An arbitrary reflection, in the cubic unit cell, $`hkl`$, is transformed to the orthorhombic unit cell as $`h+l`$, $`2k`$, $`hl`$. Thus orthorhombic reflections with $`h^{}`$ and $`l^{}`$ both even or both odd and $`k^{}`$ even stem from the original cubic unit cell. Although the fact that the $`h^{}k^{}l^{}=`$eoe reflections do not originate from the original cubic planes, they are allowed in the $`Pnma`$ symmetry. Reflections $`h^{}k^{}l^{}`$ that satisfy these conditions have higher intensities than those that do not.

To double the $`Pnma`$ unit cell to a $`2a_p`$x$`2a_p`$x$`2a_p`$ unit cell, a reflection $`h^{}k^{}l^{}`$ is transformed to $`h^{\prime \prime }k^{\prime \prime }l^{\prime \prime }`$ as $`h^{}+l^{}`$, $`k^{}`$, $`h^{}l^{}`$. Please note that we do not explicitly mention the different orientations of the orthorhombic unit cell. If the regarded $`b`$ axis is not parallel with $`k^{\prime \prime }`$, than the indices should be cycled appropriately. The transformation yields that double cubic reflections with $`k^{\prime \prime }=k^{}=2n`$ originate from single cubic reflections, with $`k=n`$. However, $`k^{\prime \prime }=2n+1`$ reflections originate from single cubic non-integer indices $`k`$, and will therefore have less intensity. This explains why the measured $`h^{\prime \prime }k^{\prime \prime }l^{\prime \prime }=`$eee reflections are strongest and eoe are one order of magnitude less strong. We have seen that $`h^{\prime \prime }=h^{}+l^{}`$ and $`l^{\prime \prime }=h^{}l^{}`$ thus $`h^{\prime \prime }`$ and $`l^{\prime \prime }`$ are always both even or both odd. Therefore eoe has only contribution of the fractions with $`k^{\prime \prime }`$ corresponding to the orthorhombic $`b`$ axis. oee and eeo originate from fractions with their $`b`$ axis parallel to $`h^{\prime \prime }`$ and $`l^{\prime \prime }`$, respectively. The fact that our observed eoe is stronger than oee and eeo, indicates that, by chance, the observed $`b`$ axis, or $`k^{\prime \prime }`$ direction, is parallel with the $`b`$ axis or $`k^{}`$ direction of the largest twin fraction.

### B Refinement

We determined both the Ca content and the volume fractions of the twins by measuring a full hemisphere of $`hkl`$ space. We compared the results with a refinement on only one octant of these measured reflections. The outcome was equivalent within the error bars. We refined again fixing the twin fractions and Ca concentration to the values found for the first refinement on the large dataset, and allowing only the atomic co-ordinates and anisotropic displacement parameters (adp’s) to change. This resulted in the same structural model as found with the complete structure determination on the largest dataset. We concluded that apparently the refinement with fixed fractions and Ca concentration is insensitive for the size of the dataset. This can be understood if we view the symmetry relations and Friedel pairs. Due to the Patterson symmetry, $`mmm`$ the relation $`F(\overline{h}kl)=F(hkl)`$ is valid for all structure factors. Friedel pairs are reflections that have intrinsically the same structure factor, related thus by symmetry and ignoring absorption and anomalous scattering effects. Although we have six different orientations, one octant chosen for a particular orientation will still contain data of all fractions that can be transformed in such a way that every $`h^{}k^{}l^{}`$ is included for all three $`b`$ directions.

## VI Conclusion

We have shown that crystals of La1-xCaxMnO3, with $`x=0.19`$, have the orthorhombic $`Pnma`$ space group, although they appear in single crystal diffractometry metrically cubic. The metric cubic appearance results from twinning in which the doubled $`b`$ axis of the $`Pnma`$ cell is oriented along all three cubic axes of the single cubic unit cell during crystal growth or cooling. This twin model allows successful structure refinements of single crystal x-ray diffraction data.

## A Data analysis

Here we describe the effect of the twin model on the visibility of the standard $`Pnma`$ reflection conditions as well as a detailed analysis of the intensity distribution of a model crystal with $`Pnma`$ symmetry and of the measured reflections in double cubic setting. The first step in space group determination is to look for systematic absences in the list of reflections. We observe no intensity for the reflections shown in table I. Note that the indices are in double cubic setting and therefore allow cyclic permutation, i.e. $`h00`$, $`0k0`$ and $`00l`$ are all represented by $`h00`$. These reflection conditions allow only the space groups $`P2_13`$ and $`P4_232`$. However, these space groups have as reflection condition only $`h00`$ : $`h=2n`$. This suggests that we have more information than can be attributed to these space groups. Attempts to determine the structure using these space groups were unsuccessful. Furthermore, the observed and, for cubic systems, anomalous, reflection conditions suggested that we might be looking at a twinned crystal. In the next paragraph, we present the transformation of the reflection conditions for orthorhombic $`Pnma`$ according to the presented twin model.

### 1 Extinctions

The $`Pnma`$ reflection condition $`0k^{}0`$ : $`k^{}=2n`$ is transformed to double cubic as $`0k^{\prime \prime }0`$ : $`k^{\prime \prime }=2n`$ and cyclic permutation. This corresponds to the reflection condition $`h^{\prime \prime }00`$ : $`h^{\prime \prime }=2n`$ in table I. $`Pnma`$ reflections of the type $`h^{}0h^{}`$ are transformed to $`h^{\prime \prime }00`$, with $`h^{\prime \prime }=2h^{}`$. This implies that $`h^{}0h^{}`$ never contributes intensity on reflections $`h^{\prime \prime }00`$, with $`h^{\prime \prime }`$ odd. Conversely, if we do not observe intensity on all reflections $`h^{\prime \prime }00`$ : $`h^{\prime \prime }2n`$, then all of the constituting $`Pnma`$ reflections should be absent.

To transform and understand the other reflection conditions of $`Pnma`$ we have to take into account that an orthorhombic, and also a cubic, unit cell implies that $`hkl`$ is equivalent upon sign reversal of each of the Miller indices. The $`Pnma`$ reflection condition $`h^{}00`$ : $`h^{}=2n`$ becomes $`h^{\prime \prime }0h^{\prime \prime }`$ : $`h^{\prime \prime }=2n`$ and similarly $`00l^{}`$ transforms to $`l^{\prime \prime }0\overline{l^{\prime \prime }}`$. But as we just stated, $`h^{\prime \prime }0h^{\prime \prime }`$ is identical to $`h^{\prime \prime }0\overline{h^{\prime \prime }}`$. If we combine reflections then a reflection condition is fulfilled if any of the contributing parts has intensity. To obtain extinction, all contributing reflections should be extinct. As $`00l^{}`$ has no extinction condition, i.e. intensity for all $`l^{}`$, they will contribute to all $`h^{\prime \prime }0h^{\prime \prime }`$/$`h^{\prime \prime }0\overline{h^{\prime \prime }}`$ reflections, and the $`Pnma`$ reflection condition $`h^{}00`$ : $`h^{}=2n`$ is masked.

The overlap of different reflections occurs also for $`0k^{}h^{}`$ and $`h^{}k^{}0`$. These transform to $`h^{\prime \prime }k^{\prime \prime }\overline{h^{\prime \prime }}`$, and $`h^{\prime \prime }k^{\prime \prime }h^{\prime \prime }`$. Therefore, we cannot disentangle $`0k^{}h^{}`$ and $`h^{}k^{}0`$ in the measured double cubic cell. The $`Pnma`$ reflection conditions for $`0k^{}l^{}`$ and $`h^{}k^{}0`$ are $`k^{}`$+$`l^{}=2n`$ and $`h^{}=2n`$, respectively. $`h^{\prime \prime }k^{\prime \prime }h^{\prime \prime }`$ and $`h^{\prime \prime }k^{\prime \prime }\overline{h^{\prime \prime }}`$ are found to be extinct for $`h^{\prime \prime }2n`$ and $`k^{\prime \prime }=2n`$ in the double cubic setting. To have these reflection extinct, all contributing $`Pnma`$ reflection conditions must be extinct. This implies that *both* $`h^{}=2n`$ *and* $`k^{}`$+$`h^{}=2n`$, the reflection conditions for the contributing $`0k^{}h^{}`$ and $`h^{}k^{}0`$, may not be fulfilled. Thus $`h^{\prime \prime }2n`$ and this yields $`k^{\prime \prime }=2n`$.

Now we consider reflections of the type $`h^{\prime \prime }k^{\prime \prime }h^{\prime \prime }`$ that are not extinct. Our twin model yields six fractions that can contribute to the considered reflections.

1. $`h^{\prime \prime }=2n`$ and $`k^{\prime \prime }=2n`$. Reflections eee have contributions from the four fractions with the doubled $`b^{}`$-axis parallel to either $`a^{\prime \prime }`$ or $`c^{\prime \prime }`$, this is also true for reflections ooo. In addition, we also have contributions from the two fractions with $`b^{}`$ parallel to $`b^{\prime \prime }`$, *i.e.* from $`h^{}k^{}0`$ and $`0k^{}h^{}`$.
2. $`h^{\prime \prime }2n`$ and $`k^{\prime \prime }=2n`$. Reflections ooo, in addition to the four fractions with the doubled $`b^{}`$-axis parallel to either $`a^{\prime \prime }`$ or $`c^{\prime \prime }`$, have one additional contribution from $`0k^{}h^{}`$. In the special case that $`h^{\prime \prime }=k^{\prime \prime }`$, for both eee and ooo, then two of the four fractions, parallel to $`a`$ and $`c`$ contribute.
3. $`h^{\prime \prime }=2n`$ and $`k^{\prime \prime }2n`$. Reflection eoe, have only a contribution from the $`h^{}k^{}0`$ reflection, with $`b^{}`$ parallel to the observed $`b^{\prime \prime }`$ axis.

We have shown that the intensity of observed $`h^{\prime \prime }k^{\prime \prime }h^{\prime \prime }`$ reflections depend on odd or even. We conclude that $`h^{\prime \prime }k^{\prime \prime }h^{\prime \prime }`$ reflections are extinct for $`OEO`$, while the other reflections have contributions from up to six twin fractions.

### 2 Intensity in double cubic setting compared with $`Pnma`$

We can also learn something about the twinned origin by studying the intensity distribution. Our twin model provides a natural explanation for the hierarchy in intensities of the reflections. First, we consider the known transition from the single cubic unit cell to the orthorhombic one to elucidate the patterns that are inherent to perovskites. Then we proceed to show how these patterns are convoluted in the double cubic, twinned model.

The intensity distribution of a normal $`Pnma`$ perovskite is shown in Fig. 4. In conventional orthorhombic $`Pnma`$ the reflections $`h^{}k^{}l^{}`$ with $`h^{}+l^{}=2n`$ are stronger than $`h^{}+l^{}2n`$. From these, the $`k=2n`$ reflections are stronger than the $`k2n`$ ones. This can be attributed to their origin in the 3.9 Å cubic structure. The $`h^{}k^{}l^{}`$ reflections in orthorhombic setting originate from $`hkl`$ in the single cubic cell for $`h^{}=h`$+$`l`$ and $`l^{}=l`$-$`h`$. Therefore, $`h^{}`$ and $`l^{}`$ will always be both even or both odd. In other words if $`h^{}`$-$`l^{}2n`$ than the reflections originate from a crystal plane in the single cubic setting with a non-integer Miller index, *i.e.* a superlattice reflection and their intensities are usually weak. The same argument can be applied to $`k`$. As $`k^{}=2k`$, thus $`k^{}=2n`$ originate from a regular crystal plane, $`k^{}2n`$ is the superlattice reflection. These $`k^{}2n`$ however appear to have somewhat more intensity than the $`h^{}`$-$`l^{}2n`$ reflections.

We observe the following patterns in the intensity distribution in the double cubic cell. First, reflections eee are generally the strongest, eoo usually absent, or very weak. If we ignore the allowed permutations in cubic setting, the following hierarchy can be made. eoe reflections were much stronger than eeo, which are again stronger than oee. Likewise, oeo reflections were, although weak, stronger than ooe and again stronger than eoo.

In conventional $`Pnma`$ structure we found that if $`h^{}`$+$`l^{}2n`$ than reflections with $`k^{}2n`$ are stronger than those with $`k^{}=2n`$. In our notation, for the $`Pnma`$ structure this yields (eee=oeo) $`>>`$ (ooo=eoe) $`>>`$ (ooe=eoo) $`>`$ (oee=eeo). So there are roughly four intensity groups. We can correlate these four groups with the four types of reflections we measured in the supercubic setting, with the following intensity hierarchy: eee$`>>`$eeo$`>>`$ooo$`>`$eoo.

We consider here the transformation from the orthorhombic to the double cubic setting for the in-plane indices, $`h^{}`$ and $`l^{}`$. Double cubic indices are calculated with $`h^{}`$+$`l^{}`$ and $`h^{}`$-$`l^{}`$, therefore $`h^{\prime \prime }`$ and $`l^{\prime \prime }`$ will be both even or both odd. Reflections of the type e?o or o?e will have no intensity of the twin with $`b^{}`$ parallel to $`b^{\prime \prime }`$, as they are not originating from an ’integer’ crystal plane. But we have a mixing of different orientations of the $`b^{}`$ axis along the three cubic axes $`a^{\prime \prime }`$, $`b^{\prime \prime }`$ and $`c^{\prime \prime }`$. We can cycle the indices of these reflections such that we get either e?e or o?o, and from the corresponding twin, $`b^{}`$ parallel to $`a^{\prime \prime }`$ or $`c^{\prime \prime }`$ we do have intensity on these reflections. eoe and oeo will only have contributions from this particular setting. eee and ooo can be cycled and in general have contributions from all six orientations.

### 3 distribution of the twins

If we consider that the main part of the intensity in the double cubic setting originates from the largest twin fraction, we can search for the orientation of the largest fraction. That $`b^{}`$ axis is expected to be oriented along one of the measured double cubic axis. We can not differentiate the two fractions with parallel $`b^{}`$ and perpendicular $`a^{}`$ and $`c^{}`$. We considered groups of reflections with cycled indices, having either one odd, oee, or one even index, eoo. Roughly, the intensities of the reflections within these groups occurred with ratio 70:20:10. This suggested that one twin had 70% of the volume, the others 20% and 10%. Now we sorted these reflections to find the corresponding orientations. The sort parameter is the observed double cubic axis that corresponded with the odd indices in oee, and vice versa. Only cycled reflections that occurred three times with measurable intensity were taking into account. We found that $`I_{a^{\prime \prime }}:I_{b^{\prime \prime }}:I_{c^{\prime \prime }}=5:80:15`$. The constraint that all three intensities have positive values, ignores the weakest reflections. It also ignores Bijvoet pairs, as we only measured $`20<l^{\prime \prime }<2`$. If we indeed had a cubic system then these variations should be zero within standard deviation. The distribution of the intensities with respect to the different axis, strongly suggests that the crystal is twinned.

## B Implementation in refinement

Here we describe the model and the application by SHELXL. The standard refining programs can work with merohedral twin models. However, in our case, the reflections do not coincide with one orthorhombic unit cell. Every observed reflection has contributions of up to six twin fractions. We measured both in the double cubic setting as in the orthorhombic setting of the main twin fraction. We transformed the $`h^{\prime \prime }k^{\prime \prime }l^{\prime \prime }`$ indices to the six possible twin orientations. Three of the possibilities are given by

$`h^{}k^{}l^{}={\displaystyle \frac{1}{2}}(h^{\prime \prime }+l^{\prime \prime }),k^{\prime \prime },{\displaystyle \frac{1}{2}}(h^{\prime \prime }l^{\prime \prime })`$ (B1)
$`h^{}k^{}l^{}={\displaystyle \frac{1}{2}}(l^{\prime \prime }+k^{\prime \prime }),h^{\prime \prime },{\displaystyle \frac{1}{2}}(l^{\prime \prime }k^{\prime \prime })`$ (B2)
$`h^{}k^{}l^{}={\displaystyle \frac{1}{2}}(k^{\prime \prime }+h^{\prime \prime }),l^{\prime \prime },{\displaystyle \frac{1}{2}}(k^{\prime \prime }h^{\prime \prime })`$ (B3)
The other three can be acquired by changing $`h^{}`$ for $`l^{}`$ and $`l^{}`$ for -$`h^{}`$.

A new software program TWINSXL was developed to transform the standard data file, HKLS, by using the appropriate transformation matrices from a second input file. The new data file constitutes of lines of $`h,k,l`$, intensity plus standard deviation and twin fraction number. We used the ”HKLF 5” option of SHELXL to refine data. The refinement uses a crystal model for the orthorhombic structure with the normal, adjustable variables. Five variables for the twin fractions were added, the sixth fraction is calculated as the complement of the other five fractions. The sum of the appropriate calculated intensities for all fractions was compared with the observed integrated intensities.


format: text

text:
# 1 Introduction

## 1 Introduction

A symplectic instanton bundle on $`\text{}_{\text{}}^{2n+1}`$ is a bundle of rank $`2n`$ defined as the cohomology bundle of a well known monad (see Definition 2.2).

In \[ADHM78\] it was shown that instanton bundles on $`\text{}^3`$ satisfying a reality condition correspond to self dual Yang Mills $`Sp(1)`$-connections over the $`4`$-dimensional sphere $`S^4=\text{}_{\text{}}^1`$. This correspondence was generalized by Salamon (\[Sal84\]) who showed that instanton bundles on $`\text{}^{2n+1}`$ which are trivial on the fiber of the twistor map $`\text{}^{2n+1}\text{}_{\text{}}^n`$ correspond to $`Sp(n)`$-connections which minimize a certain Yang Mills functional over $`\text{}_{\text{}}^n`$. We denote by $`MI_{\text{}^{2n+1}}(k)`$ the moduli space of symplectic instanton bundles on $`\text{}^{2n+1}`$ with $`c_2=k`$ (see Definition 2.4). Few things are known on $`MI_{\text{}^{2n+1}}(k)`$ in general. It is known (\[ADHM78\]) that it has a component of dimension $`8k3`$ for $`n=1`$, it is smooth for $`n=1`$ and $`k5`$ (\[KO99\]) and it is expected to be singular and reducible for $`n2`$ and $`k4`$ (see \[AO00\]).

Our main result, proved in section 2, is the following

###### Theorem 1.1.

$`MI_{\text{}^{2n+1}}(k)`$ is affine.

In particular $`MI_{\text{}^{2n+1}}(k)`$ does not contain any complete subvariety of positive dimension. For other results of this type see \[HH86\]. The technique we use to prove the main theorem is to exhibit $`MI_{\text{}^{2n+1}}(k)`$ as the GIT-quotient of an affine variety $`𝒬^0`$ and then use standard results about invariant theory. The fact that $`𝒬^0`$ is affine is a consequence of the existence of a invariant of multidimensional matrices representing the instanton bundles, which generalizes the hyperdeterminant (see \[GKZ94\] and \[AO99\]).

The first named author would like to thank the Dipartimento di Matematica U. Dini for their hospitality and support at the time of the preparation of this paper.

## 2 Notations and preliminaries

We will start fixing some notation and recalling some facts about $`k`$-instanton bundles on $`\text{}^{2n+1}=\text{}(V)`$, where $`V`$ is a complex vector space of dimension $`2n+2`$ (See for instance \[OS86\] and \[AO94\]).

###### Notation 2.1.

$`𝒪(d)=𝒪_{\text{}^{2n+1}}(d)`$ denotes the invertible sheaf of degree $`d`$ on $`\text{}^{2n+1}`$ and for any coherent sheaf $`E`$ on $`\text{}^{2n+1}`$ we denote $`E(d)=E𝒪_{\text{}^{2n+1}}(d)`$ .

###### Definition 2.2.

A symplectic instanton bundle $`E`$ over $`\text{}^{2n+1}=\text{}(V)`$ is a bundle of rank $`2n`$ which appears as a cohomology bundle of a monad

$$(1)I^{}𝒪(1)\stackrel{A}{}W𝒪\stackrel{A^t}{}I𝒪(1)$$

where $`(W,J)`$ is a symplectic complex vector space of dimension $`2n+2k`$ and $`I`$ is a complex vector space of dimension $`k`$.

We do not assume in the definition that $`E`$ is stable, so we have to recall some results.

The monad condition means that $`A`$ is injective (as a bundle morphism), $`A^t`$ is surjective and $`\mathrm{im}A\mathrm{ker}A^t`$ so that $`E\mathrm{ker}A^t/\mathrm{im}A`$. The fact that the map $`W𝒪\stackrel{A^t}{}I𝒪(1)`$ is surjective, is equivalent to the fact that the matrix $`AHom(V^{}I^{},W)`$ representing $`E`$ is nondegenerate according to \[GKZ94\] (see Definition 2.3 for the precise definition).

$`Hom(V^{}I^{},W)`$ contains the subvariety $`𝒬`$ given by matrices $`A`$ for which the sequence (2.2) is a complex, that is, such that $`A^tJA=0`$. $`GL(I)\times Sp(W)`$ acts on $`𝒬`$ by $`(g,s)A=sAg`$.

###### Definition 2.3.

A matrix $`AHom(V^{}I^{},W)`$ is called degenerate if the multilinear system $`A(vi)=0`$ has a solution such that $`0vV^{}`$ and $`0iI^{}.`$

By \[GKZ94\], Theorem 14.3.1, this is equivalent to the standard definition of degeneracy given in chapter 14.1 of \[GKZ94\]. It is easy to check that nondegenerate matrices fill an irreducible subvariety $`N`$ of $`Hom(V^{}I^{},W)`$ of codimension $`k`$ (see \[WZ96\]). Hence only in the case $`k=1`$ it is well defined a hyperdeterminant according to \[GKZ94\]. In the next section we will define a $`SL(I)\times Sp(W)`$-invariant on $`Hom(V^{}I^{},W)`$ called $`D`$ which generalizes the hyperdeterminant and it is suitable for our purposes.

It was shown in \[AO94\] that all instanton bundles are simple, so that they carry a unique symplectic form. Moreover for $`n=1,2`$ it was proved in \[AO94\] that all instanton bundles are stable, and it is expected that the same result is true for $`n3`$. In \[BH78\] it was essentially proved that there is a natural 1:1 correspondence between

i) isomorphism classes of symplectic instanton bundles, and

ii) orbits of $`GL(I)\times Sp(W)`$ on the open subvariety $`𝒬^0`$ of $`𝒬`$ given by nondegenerate matrices.

In fact, using the quoted results of \[AO94\], one can see that \[BH78\]; Section $`4`$ and Theorem at page 19, adapt literally to our situation. Moreover, \[BH78\]; Remark 2 at page 19 shows that, if $`G`$ denotes the quotient of $`GL(I)\times Sp(W)`$ by $`\pm (id,id)`$, then $`G`$ acts freely on $`𝒬^0`$. In particular all points of $`𝒬^0`$ are stable for the action of $`GL(I)\times Sp(W)`$.

###### Definition 2.4.

The GIT-quotient $`𝒬^0/GL(I)\times Sp(W)`$ is denoted by $`MI_{\text{}^{2n+1}}(k)`$ and it is called the moduli space of symplectic $`k`$-instanton bundles on $`\text{}^{2n+1}`$. It is a geometric quotient.

The above discussion shows that $`MI_{\text{}^{2n+1}}(k)`$ coincides for $`n=1,2`$ with the open subset $`_{\text{}^{2n+1}}(k)`$ of the Maruyama scheme of symplectic stable bundles on $`\text{}^{2n+1}`$ of rank $`2n`$ and Chern polynomial $`\frac{1}{(1t^2)^k}`$ which are instanton bundles (this is an open condition because by Beilinson theorem is equivalent to certain vanishing in cohomology, see \[OS86\] ). In particular our notation for $`MI_\text{}^3(k)`$ is consistent with the usual one. For $`n3`$ it is expected that the same result is true, but at present we can only say that $`_{\text{}^{2n+1}}(k)`$ is an open subset of $`MI_{\text{}^{2n+1}}(k)`$.

## 3 The invariant $`D`$ and the proof of the main result

First of all we remark that the vector spaces $`WS^nI`$ and $`VS^{n+1}I`$ have the same dimension $`(2n+2k)\left(\genfrac{}{}{0pt}{}{k+n1}{n}\right)=(2n+2)\left(\genfrac{}{}{0pt}{}{k+n}{n+1}\right)`$. We can construct from

$$W\stackrel{A^t}{}VI$$

the following morphisms

$$\begin{array}{c}A^tid_{S^nI}:WS^nIVIS^nI,\\ id_V\pi :VIS^nIVS^{n+1}I,\end{array}$$

where $`\pi `$ is the natural projection, and we consider the composition

$$(2)\mathrm{\Delta }_A=(id_V\pi )(A^tid_{S^nI}):WS^nIVS^{n+1}I.$$

###### Definition 3.1.

Let $`AHom(V^{}I^{},W)`$. We define $`D(A)`$ to be the usual determinant of the morphism $`\mathrm{\Delta }_A`$ in (3) induced by $`A`$.

Notice that

$$D:Hom(V^{}I^{},W)(detI)^\alpha (detV)^\beta $$

where $`\alpha =2\left(\genfrac{}{}{0pt}{}{k+n}{n}\right)`$ and $`\beta =\left(\genfrac{}{}{0pt}{}{k+n}{n+1}\right)`$ is a $`GL(V)\times GL(I)\times Sp(W)`$-equivariant map and $`D(A)=0`$ defines a homogeneous hypersurface of degree $`(2n+2k)\left(\genfrac{}{}{0pt}{}{k+n1}{n}\right)=(2n+2)\left(\genfrac{}{}{0pt}{}{k+n}{n+1}\right)`$. After a basis has been fixed in each of the vector spaces $`V`$, $`I`$ and $`W`$, the map $`D`$ can be seen as a $`SL(V)\times SL(I)\times Sp(W)`$-invariant.

In fact this definition generalizes the hyperdeterminant of boundary format as introduced in Theorem 14.3.3 of \[GKZ94\].

###### Lemma 3.2.

If $`A`$ is degenerate then $`D(A)=0`$.

Proof. There are $`0vV^{}`$ and $`0iI^{}`$ such that $`A(vi)=0`$. Hence $`vS^{n+1}iV^{}S^{n+1}I^{}`$ goes to zero under the dual of (3). ∎

If $`A`$ is nondegenerate we get $`D(A)0`$ only in the case $`k=1`$ and in general it can happen that $`D(A)=0`$, because the codimension of $`N`$ is $`k`$. Our main technical result is the following

###### Theorem 3.3.

If $`A`$ defines an instanton (that is $`A`$ belongs to $`𝒬^0`$) then $`D(A)0`$.

Proof. From (2.2) we get the exact sequence

$$(3)0KW𝒪I𝒪(1)0.$$

The $`(n+1)`$-th wedge power twisted by $`𝒪(n)`$ gives the exact sequence

$$0^{n+1}K(n)^{n+1}W(n)\mathrm{}$$

$$\mathrm{}^2WS^{n1}I(1)WS^nIS^{n+1}I(1)0$$

where the $`H^0`$ of the last morphism corresponds to $`\mathrm{\Delta }_A`$ in (3). Taking cohomology, it is enough to prove

$$(4)H^n(^{n+1}K(n))=0.$$

The $`(n+1)`$-th wedge power twisted by $`𝒪(n)`$ of the sequence

$$0I^{}𝒪(1)KE0$$

gives the sequence

$$\begin{array}{c}0S^{n+1}I^{}K(2n1)\mathrm{}^{n1}KS^2I^{}(n2)\\ ^nKI^{}(n1)^{n+1}K(n)^{n+1}E(n)0.\end{array}$$

In order to prove (3), taking cohomology, we need $`H^{n+i}(^{ni}K(ni1))=0`$ for $`i=0,\mathrm{},n`$ and $`H^n(^{n+1}E(n))=0`$. The first group of vanishing is easily obtained by taking suitable wedge powers of (3). The crucial point to get the last vanishing is the isomorphism $`^{n+1}E^{n1}E`$, that it is true because $`E`$ is a rank $`2n`$ vector bundle with $`c_1=0`$. From the sequence

$$0S^{n1}I^{}(2n1)S^{n2}I^{}K(2n)\mathrm{}$$

$$\mathrm{}^{n1}K(n)^{n1}E(n)0,$$

in order to prove $`H^n(^{n1}E(n))=0`$, we only need to see that

$$H^{n+i}(^{n1i}K(ni))=0\text{ for }i=0,\mathrm{},n$$

which follows by using the exact sequence (3) exactly as above. ∎

###### Corollary 3.4.

$`MI_{\text{}^{2n+1}}(k)`$ is affine.

Proof. By the Theorem 3.3 we get that $`𝒬N=𝒬^0=𝒬\{D=0\}`$ is affine. It follows that $`MI_{\text{}^{2n+1}}(k)`$ is affine too, because it is the quotient of an affine variety by a reductive group, see e.g. \[PV89\], section 4.4. ∎

###### Remark 3.5.

The invariant $`D`$ is meaningful even in the case $`n=0`$. In this case it corresponds to the usual determinant of the map $`\text{}^{2k}\text{}^2\text{}^k`$. For example for $`n=0`$ and $`k=2`$ the degenerate $`2\times 2\times 4`$ matrices fill a variety of codimension $`2`$ and degree $`12`$ (\[BS78\]) in $`\text{}^{15}`$ whose ideal is generated by one quartic (which is our invariant $`D`$), $`10`$ sextics and one octic. We remark that the case $`2\times 2\times 3`$ is of boundary format. The case $`2\times 2\times 5`$ is interesting. Here degenerate matrices fill a variety of codimension $`3`$ and degree 20, and its ideal is generated (at least) by $`5`$ quartics, $`50`$ sextics and $`12`$ octics. The $`5`$ quartics define a variety of codimension $`2`$ and degree $`10`$. Hence in this case no analog of the invariant $`D`$ can exist.

## 4 Instanton bundles with structure group $`GL(2n)`$

###### Definition 4.1.

A $`k`$-instanton bundle $`E`$ on $`\text{}^{2n+1}`$ is the cohomology bundle of a monad

$$(5)K𝒪(1)\stackrel{A}{}W𝒪\stackrel{B}{}I𝒪(1)$$

where $`W`$ is a complex vector space of dimension $`2n+2k`$ and $`I,K`$ are complex vector spaces of dimension $`k`$.

Notice that $`E`$ is not necessarily symplectic and that this notion is a true generalization of the above one only for $`n2`$ because all rank $`2`$ bundles on $`\text{}^3`$ with $`c_1=0`$ are symplectic.

Let $`(A,B)Hom(KV^{},W)\times Hom(W,IV)`$ defining $`E`$. The monad condition is now equivalent to the fact that the matrices $`A`$ and $`B`$ are both nondegenerate and $`BA=0`$.

$`Hom(KV^{},W)\times Hom(W,IV)`$ contains the subvariety $`𝒫`$ given by pairs of matrices $`(A,B)`$ for which the sequence (4.1) is a complex, that is, such that $`BA=0`$. $`GL(I)\times GL(K)\times GL(W)`$ acts on $`𝒫`$ by $`(a,b,c)(A,B)=(cAb,aBc^1)`$.

As in the previous section we can see that there is a natural 1:1 correspondence between

i) isomorphism classes of instanton bundles, and

ii) orbits of $`GL(I)\times GL(K)\times GL(W)`$ on the open subvariety $`𝒫^0`$ of $`𝒫`$ given by pairs of nondegenerate matrices.

Moreover, if we denote by $`H`$ the quotient of $`GL(I)\times GL(K)\times GL(W)`$ by $`(\lambda id,\lambda ^1id,\lambda id)`$, then $`H`$ acts freely on $`𝒫^0`$. In particular, all points of $`𝒫^0`$ are stable for the action of $`GL(I)\times GL(K)\times GL(W)`$.

###### Definition 4.2.

The GIT-quotient $`𝒫^0/GL(I)\times GL(K)\times GL(W)`$ is denoted by $`I_{\text{}^{2n+1}}(k)`$ and it is called the moduli space of $`k`$-instanton bundles on $`\text{}^{2n+1}`$. It is a geometric quotient.

$`I_{\text{}^{2n+1}}(k)`$ coincides for $`n=1,2`$ with the open subset $`_{\text{}^{2n+1}}(k)`$ of the Maruyama scheme of stable bundles on $`\text{}^{2n+1}`$ of rank $`2n`$ and Chern polynomial $`\frac{1}{(1t^2)^k}`$ which are instanton bundles. For $`n3`$ we can say that $`_{\text{}^{2n+1}}(k)`$ is an open subset of $`I_{\text{}^{2n+1}}(k)`$. We remark that $`MI_\text{}^3(k)=I_\text{}^3(k)`$. $`_{\text{}^{2n+1}}(k)`$ is known to be singular for $`n2`$ and $`k3`$ (see \[MO97\]) and reducible for $`n4`$ (see \[AO00\]).

###### Definition 4.3.

Let $`(A,B)Hom(KV^{},W)\times Hom(W,IV)`$. We define

$$\stackrel{~}{D}(A,B):=detS(A)detR(B)$$

where $`det`$ denotes the usual determinant and $`S(A)`$, $`R(B)`$ are the following morphisms

$$\begin{array}{c}S(A):S^{n+1}KV^{}S^nKW,\hfill \\ R(B):S^nIWS^{n+1}IV,\hfill \end{array}$$

induced by $`A`$ and $`B`$ respectively, as in Definition 3.1.

###### Theorem 4.4.

If $`(A,B)`$ defines an instanton (that is $`(A,B)`$ belongs to $`𝒫^0`$) then $`\stackrel{~}{D}(A,B)0`$.

Proof. First of all we will see that $`detS(A)0`$. From (4.1) we get the exact sequence

$$(6)0K𝒪(1)W𝒪Q0.$$

The $`(n+1)`$-th wedge power twisted by $`𝒪(n2)`$ gives the exact sequence

$$\begin{array}{c}0S^{n+1}K𝒪(2n3)S^nKW𝒪(2n2)\mathrm{}\\ ^{n+1}W𝒪(n2)^{n+1}Q(n2)0\end{array}$$

where the $`H^{2n+1}`$ of the first morphism corresponds to $`S(A)`$. Hence, taking cohomology, it is enough to prove

$$H^n(^{n+1}Q(n2))=0.$$

This is shown considering the $`(n+1)`$-wedge sequence of the exact sequence

$$0EQI𝒪(1)0$$

and arguing as in the proof of Theorem 3.3.

In order to prove $`detR(B)0`$ we proceed exactly as in Theorem 3.3 and we leave the details to the reader. ∎

###### Corollary 4.5.

$`I_{\text{}^{2n+1}}(k)`$ is affine.

Proof. First of all notice that given $`(A,B)𝒫`$, if $`A`$ or $`B`$ is degenerate then $`detS(A)detR(B)=0`$. Hence, by Theorem 4.4 we get that $`𝒫^0=𝒫\{\stackrel{~}{D}=0\}`$ is affine. Therefore, by \[PV89\] section 4.4, $`I_{\text{}^{2n+1}}(k)`$ is affine too. ∎

Authors’ addresses

Laura Costa Dept. Algebra y Geometria Universitat de Barcelona Gran Via, 585 08007 Barcelona SPAIN costa@mat.ub.es    Giorgio Ottaviani Dipartimento di Matematica ”U. Dini” Università di Firenze viale Morgagni 67/A I 50134 FIRENZE ITALY ottavian@math.unifi.it


text:
# References

Lengths of the First Days of the Universe

Moshe Carmeli

Department of Physics, Ben Gurion University, Beer Sheva 84105, Israel

Email: carmelim@bgumail.bgu.ac.il

## Abstract

The early stage of the Universe is discussed and the time lengths of its first days are given. If we denote the Hubble time in the zero-gravity limit by $`\tau `$ (approximately 12.16 billion years), and $`T_n`$ denotes the length of the $`n`$-th day, then we have the very simple relation $`T_n=\tau /(2n1)`$. Hence we obtain for the first days the following lengths of time: $`T_1=\tau `$, $`T_2=\tau /3`$, $`T_3=\tau /5`$, etc.

In this Note we calculate the lengths of days of the early Universe, day by day, from the first day after the Big Bang on up to our present time. We find that the first day actually lasted the Hubble time in the limit of zero gravity. If we denote the Hubble time in the zero-gravity limit by $`\tau `$ which equals about 12.16 billion years and $`T_n`$ denotes the length of the $`n`$-th day in units of times of the early Universe, then we have a very simple relation

$$T_n=\frac{\tau }{2n1}.$$
$`(1)`$
Hence we obtain for the first few days the following lengths of time:

$$T_1=\tau ,T_2=\frac{\tau }{3},T_3=\frac{\tau }{5},T_4=\frac{\tau }{7},T_5=\frac{\tau }{9},T_6=\frac{\tau }{11}.$$
$`(2)`$
It also follows that the accumulation of time from the first day to the second, third, fourth, …, up to now is just exactly the Hubble time. The Hubble time in the limit of zero gravity is the maximum time allowed in nature.

Using Cosmological Special Relativity \[1-4\], the calculation is very simple. We assume that the Big Bang time with respect to us now was $`t_0=\tau `$, the time of the first day after that was $`t_1`$, the time of the second day was $`t_2`$, and so on. In this way the time scale is progressing in units of one day (24 hours) in our units of present time. The time difference between $`t_0`$ and $`t_1`$, denoted by $`T_1`$, is the time as measured at the early Universe and is by no means equal to one day of our time. In this way we denote the times elapsed from the Big Bang to the end of the first day $`t_1`$ by $`T_1`$, between the first day $`t_1`$ and the second day $`t_2`$ by $`T_2`$ and so on. According to the rule of the addition of cosmic times one has, for example,

$$t_6+1(day)=\frac{t_6+T_6}{1+t_6T_6/\tau ^2}.$$
$`(3)`$
A straightforward calculation then shows that

$$T_6=\frac{\tau ^2}{\tau ^2\left(\tau 6\right)\left(\tau 5\right)}=\frac{\tau ^2}{11\tau 30}.$$
$`(4)`$
In general one finds that

$$T_n=\frac{\tau ^2}{\tau ^2\left(\tau n\right)\left(\tau n+1\right)},$$
$`(5)`$
or

$$T_n=\frac{\tau }{n+\left(n1\right)n\left(n1\right)/\tau }.$$
$`(6)`$
As is seen from the last formula one can neglect the last term in the denominator in the first approximation and we get the simple Eq. (1).

From the above one reaches the conclusion that the age of the Universe exactly equals the Hubble time in vacuum $`\tau `$, i.e. 12.16 billion years, and it is a universal constant . This means that the age of the Universe tomorrow will be the same as it was yesterday or today.

But this might not go along with our intuition since we usually deal with short periods of times in our daily life, and the unexperienced person will reject such a conclusion. Physics, however, deals with measurements.

In fact we have exactly a similar situation with respect to the speed of light $`c`$. When measured in vacuum, it is 300 thousands kilometers per second. If the person doing the measurement tries to decrease or increase this number by moving with a very high speed in the direction or against the direction of the propagation of light, he will find that this is impossible and he will measure the same number as before. The measurement instruments adjust themselves in such a way that the final result remains the same. In this sense the speed of light in vacuum $`c`$ and the Hubble time in vacuum $`\tau `$ behave the same way and are both universal constants.

The similarity of the behavior of velocities of objects and those of cosmic times can also be demonstrated as follows. Suppose a rocket moves with the speed $`V_1`$ with respect to an observer on the Earth. We would like to increase that speed to $`V_2`$ as measured by the observer on the Earth. In order to achieve this, the rocket has to increase its speed not by the difference $`V_2V_1`$, but by

$$\mathrm{\Delta }V=\frac{V_2V_1}{1V_1V_2/c^2}.$$
$`(7)`$
As can easily be seen $`\mathrm{\Delta }V`$ is much larger than $`V_2V_1`$ for velocities $`V_1`$ and $`V_2`$ close to that of light $`c`$. This result follows from the rule for the addition of velocities,

$$V_{1+2}=\frac{V_1+V_2}{1+V_1V_2/c^2},$$
$`(8)`$
a consequence of Einstein’s famous Special Relativity Theory . In cosmology, we have the analogous formula

$$T_{1+2}=\frac{T_1+T_2}{1+T_1T_2/\tau ^2}$$
$`(9)`$
for the cosmic times.


text:
# Submillimeter CO emission from shock-heated gas in the L1157 outflow

## 1 INTRODUCTION

There is now increasing evidence that the extremely young outflows driven by class 0 protostars interact with the surrounding ambient gas and produce strong shock waves (e.g., Bachiller & Gómez-González 1992; Bachiller 1996). Since the shocked gas is significantly heated and compressed, various chemical reactions that cannot proceed in the cold and quiescent dark clouds are triggered off (e.g., Bachiller & Pérez Gutiérrez 1997; Garay et al. 1998). One prototypical and well studied example of an outflow with strong shocks is the well-collimated bipolar outflow in the L1157 dark cloud (e.g., Umemoto et al. 1992), which is located at 440 pc from the Sun (Viotti 1969). The blue lobe of this outflow is known to consist of two limb-brightened cavities (Gueth, Guilloteau, & Bachiller 1996); toward the tips of these two cavities, strong SiO and NH3 (3, 3) emission lines which are believed to be good tracers of shocked molecular gas are observed (Gueth, Guilloteau, & Bachiller 1998; Zhang et al. 1995; Tafalla & Bachiller 1995). Recently, Umemoto et al. (1999) have detected the highly excited NH3 emission lines of ($`J,K`$) = (5, 5) and (6, 6) in the blue lobe. They have shown that the blueshifted gas in the L1157 outflow is heated to 130–140 K, which is at least by a factor of $``$10 over the temperature of the quiescent gas even though the outflow is excited by a low-luminosity source of only $``$11 $`L_{}`$. However, their NH3 observations have been performed with rather low angular resolution, $``$72′′, and therefore their results do not tell us the spatial distribution of the highly excited gas and the variation of the gas kinetic temperature in the lobe.

In this Letter, we present the results of our new CO $`J`$=6–5, 4–3, and 3–2 observations of the high-velocity gas in the blue lobe of the L1157 outflow. Since the upper level of the CO $`J`$=6–5 transition lies 117 K above the ground, this emission line is expected to be a good tracer of the warm molecular gas. Furthermore, ratios of this line to the lower energy transitions can be used to probe the kinetic temperature and density of the high-velocity gas. Our spatial resolution of $``$14′′ is high enough to explore the spatial variation of the kinetic temperature in the outflow lobe whose extent is $``$2.

## 2 OBSERVATIONS

The observations of the CO $`J`$ = 6–5, 4–3, and 3–2 transitions were carried out with the James Clerk Maxwell Telescope (JCMT)11affiliation: The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre, on behalf of the Particle Physics and Astronomy Research Council of the United Kingdom, the Netherlands Organization for Science Research and the National Research Council of Canada in 2000 June. The half-power beam widths (HPBW) of the telescope were 8′′ at the frequency of CO 6–5 (691.473 GHz), 12′′ at that of CO 4–3 (461.041 GHz), and 14′′ at the CO 3–2 frequency (345.376 GHz). The main beam efficiencies at 691, 461, and 345 GHz were 0.3, 0.5, and 0.62, respectively. The CO 6–5 and 4–3 observations were carried out using the dual-channel SIS receiver RxW in the D-band and C-band, respectively, with an atmospheric opacity at 225 GHz of $``$0.06. The single-sideband system temperature was 7000–10000 K in the D-band and was 2500–3000 K in the C-band. We used the SIS receiver RxB3 in the single-sideband mode for the CO 3–2 observations. The system noise temperature was 520–570 K. The backend was a digital autocorrelation spectrometer (DAS) set up with a bandwidth of 500 MHz for the CO 6–5 observations and 250 MHz for the CO 4–3 and 3–2 observations. The frequency resolution of the DAS with 500 MHz bandwidth mode was 625 kHz, corresponding to a velocity resolution of 0.27 km s-1 at 691 GHz. The frequency resolution of the DAS with 250 MHz mode was 312.5 kHz, which corresponds to the velocity resolution of 0.27 km s-1 at 345 GHz and of 0.20 km s-1 at 461 GHz.

We obtained spectra toward three positions along the blueshifted lobe of the CO outflow. The position offsets from the driving source, IRAS 20386+6751 (hereafter referred to as L1157 IRS), $`\alpha `$(1950) = 20h38m39.s6, $`\delta `$(1950) = 675133′′, are (0′′, 0′′), (20′′, -60′′), and (45′′, -105′′), which were referred to as positions A, B, and C, respectively, by Mikami et al. (1992). Positions B and C correspond to the tips of two CO cavities (Gueth et al. 1996). The observed positions are shown in Figure 1 superposed on a CO ($`J`$=1–0) outflow map observed with the 45 m telescope of the Nobeyama Radio Observatory (NRO) by Umemoto et al. (private communication). We made five-point maps with 5′′ spacing centered at the three positions, A, B, and C, in the CO $`J`$=6–5. In the CO $`J`$=3–2 transition, we obtained the spectra along the line through the three positions. The sampling intervals were 6.3′′ between positions A and B and 6.8′′ between positions B and C. We used the position switching mode to obtain the spectra, and integrated 10 minutes per position in the CO $`J`$=6–5, 5 minutes in the $`J`$=4–3, and 3 minutes in the $`J`$=3–2. In the CO $`J`$=6–5 and $`J`$=4–3, the signals collected in two channels were co-added in order to improve the signal-to-noise ratio. The telescope pointing accuracy was checked by observing W75N and NGC7538 IRS1, and was better than 3′′. We have obtained the spectra of W75N and NGC7538 IRS1, and used them as intensity calibrators. The intensities and line shapes of the calibrator spectra observed by us were consistent with those of the previously observed ones in the CO 3–2 and 6–5 lines. However, the intensity and line profile of the CO 4–3 line did not agree well with the previous ones. We also compared the CO 4–3 spectrum observed at position B with that observed in 1995 which is available in the JCMT data archive. We found that our peak temperature ($``$33 K) is lower than the previous measurements ($``$40 K). Therefore, the CO 4–3 intensity values of our measurements may be underestimated by $``$20%. Typical rms noises per channel were $``$0.3 K in $`T_{\mathrm{mb}}`$ for CO 3–2, 0.9–1.2 K for CO 4–3, and 2.0–3.0 K for CO 6–5 spectra.

## 3 RESULTS AND DISCUSSION

### 3.1 Strong CO Emission Lines in the Blue Lobe

In Figure 1, we show the CO $`J`$=6–5, 4–3, and 3–2 spectra observed at positions A, B, and C along with those of the CO and 13CO $`J`$=1–0. The CO $`J`$=6–5 spectra presented in Figure 1 were obtained by convolving the data of five points centered at each position with the 14′′ Gaussian beam. The spectra of the CO and 13CO $`J`$=1–0 were observed with the NRO 45m telescope (the beam size was 15′′) by Umemoto et al. (private communication). It is shown that the CO $`J`$=6–5, 4–3, and 3–2 emission lines are extremely strong ($`T_{\mathrm{mb}}>`$ 30 K) at position B which is 63′′ ($``$0.13 pc) away from the driving source. The peak temperatures of the CO $`J`$=6–5, 4–3, and 3–2 emission lines observed at position B are 32, 33, and 37 K, respectively, which are much higher than $`T_{\mathrm{mb}}`$ 20 K in the CO $`J`$=1–0.

Figure 1 shows that the upper three transitions of CO (i.e., $`J`$=6–5, 4–3, and 3–2) have similar line profiles, although the absorption dips in the $`J`$=6–5 lines are less prominent than those of the $`J`$=4–3 and 3–2 ones. The absorption dips are likely to be caused by the foreground material with lower excitation because their velocities coincide with the cloud systemic velocity of $``$2.7 km s-1 determined from the intensity peak of the 13CO emission. A comparison of the profiles of the CO $`J`$=6–5, 4–3, and 3–2 lines with the 13CO $`J`$=1–0 ones implies that the upper three transitions of CO lines are emitted by the outflowing gas even in the velocity ranges close to the cloud systemic velocity. The spectra of CO $`J`$=1–0 show different characteristics from those of the upper transitions; no absorption dips, smaller velocity extent at positions B and C, and less prominent blueshifted wing at position A. The shapes of the CO $`J`$=1–0 spectra are also affected by the self-absorption because their peak velocities are blueshifted by $``$0.8 km s-1 with respect to the cloud systemic velocity.

### 3.2 Physical condition of the outflowing gas

In Figure 2, we show the line intensity ratios of CO 6–5/3–2 and 3–2/1–0 measured at position B. We used the CO $`J`$=6–5 data convolved with the 14′′ beam in order to match their angular resolution with the CO $`J`$=3–2 and 1–0 data. Because of the uncertainty of the CO $`J`$=4–3 intensity values, the ratios against the CO $`J`$=4–3 are not used for the following analyses. It is shown that the 6–5/3–2 line ratio is almost constant to $``$0.7 throughout the blueshifted wings except the innermost velocities of a few km s-1 from the cloud systemic velocity. On the other hand, the 3–2/1–0 ratio increases from $``$3 at $`V_{\mathrm{LSR}}`$ = 0 km s-1 to $``$6 at $`V_{\mathrm{LSR}}`$ = $``$7 km s-1. In the following analyses, we use the intensity weighted value of $``$4 as a 3–2/1–0 ratio at position B. At positions A and C, variations of the 6–5/3–2 and 3–2/1–0 ratios as a function of the velocity show the same trends as those of position B. At these two positions, typical values of the 6–5/3–2 and 3–2/1–0 ratios are $``$0.5 and $``$3, respectively.

In order to derive the physical parameters of the outflowing gas, we carried out Large Velocity Gradient (LVG) statistical equilibrium calculations (Goldreich & Kwan 1974; Scoville & Solomon 1974). Detailed descriptions of the calculations are presented in Appendix A of Oka et al. (1998). We adopted a CO to H2 abundance ratio $`X`$ of 10-4. The velocity gradient $`dv/dr`$ is assumed to be 180 km s-1 pc-1, which is derived from the velocity dispersion of the CO $`J`$=3–2 line, 15 km s-1, and the lobe width of 0.085 pc determined from the half-power contour of the CO $`J`$=1–0 map presented in Figure 1.

Figure 3 displays the CO 6–5/3–2 and 3–2/1–0 intensity ratios as a function of density and kinetic temperature. It is shown that the line ratios observed at positions A and C correspond to $`n(\mathrm{H}_2)=2\times 10^4`$ cm-3 and $`T_{\mathrm{kin}}`$ = 75 K, and those at position B to $`n(\mathrm{H}_2)=3\times 10^4`$ cm-3 and $`T_{\mathrm{kin}}`$ = 120 K. If we assume a larger (smaller) value as the velocity gradient, the solutions of the LVG calculations give us slightly higher (lower) densities and lower (higher) kinetic temperatures. When we adopted the lobe widths determined from the interferometric CO map of Gueth et al. (1996; 0.064 pc at position B and 0.032 pc at position C), the CO $`J`$=3–2 line widths at positions B and C (15 km s-1 and 7.5 km s-1, respectively) give a velocity gradient of $`dv/dr`$ $``$ 230 km s-1 pc-1. In this case, we obtained $`n(\mathrm{H}_2)=3\times 10^4`$ cm-3 and $`T_{\mathrm{kin}}`$ = 65 K for positions A and C, and $`n(\mathrm{H}_2)=3\times 10^4`$ cm-3 and $`T_{\mathrm{kin}}`$ = 100 K for position B. The kinetic temperature at position B is consistent with the rotational temperature derived from the NH3 ($`J`$, $`K`$) = (3, 3) to (6, 6) data by Umemoto et al. (1999). On the other hand, the kinetic temperature at position C derived from the CO data is lower than the NH3 rotational temperature at the same position. Because of the large beam size of 72′′, the NH3 spectra observed at position C by Umemoto et al. may have been contaminated by the higher temperature component located at the position B.

The LVG calculations indicate that the density of the outflowing gas traced by the CO emission is typically in the range of (1–3) $`\times `$ 104 cm-3 and does not vary among the three observed positions although the density at position B is slightly higher than those of the other two positions. Because of this rather low density, the CO transitions are not thermalized. Especially, the CO $`J`$=1–0 transition falls in the parameter range of superthermally excited or inversely populated. Using the model parameters derived from the LVG analysis and the observed total integrated intensities of the $`J`$=6–5, 3–2, and 1–0 lines, we estimated the total energy released by the rotational transitions of CO to be 1.2$`\times `$1031 erg s-1 at position B and 1.4$`\times `$1030 erg s-1 at position C. The LVG model predicts that the $`J`$=6–5 line is the most luminous CO transition at position B and accounts for $``$25% of the total CO luminosity. On the other hand, the $`J`$=1–0 transition contributes only 0.5 % of the CO luminosity at the same position. The $`J`$=4–3 and 3–2 lines account for 13% and 5% of the total CO luminosity at this position, respectively. At position C with lower gas temperature, the $`J`$=6–5, 4–3, 3–2, and 1–0 transitions contribute 26 %, 23 %, 11 %, and 1 % of the total CO luminosity, respectively. It should be noted that the values estimated here (i.e., the energy released by the CO lines and the relative luminosity of each transition) depend little on the assumed $`dv/dr`$ value.

### 3.3 Temperature variation along the blue lobe

In Figure 4a, we show the position-velocity diagram of the CO $`J`$= 3–2 emission along the line through A, B, and C. It is shown that the velocity dispersion of the CO $`J`$=3–2 line increases with increasing distance from the protostar toward position B, which corresponds to the tip of the CO cavity (Gueth et al. 1996). This velocity pattern can be explained if the cavity was created by the passage of bow shock (e.g., Bence, Richer, & Padman 1996; Smith, Suttner, & Yorke 1997). The velocity dispersion suddenly decreases downstream from position B, and slightly increases again toward position C. Such a double-humped velocity dispersion shown in Figure 4a can be explained in terms of multiple eruption events, each of which having produced a bow shock.

In Figure 4b, we plot the variation of the CO 3–2/1–0 intensity ratio averaged over the blueshifted wing ($`V_{\mathrm{LSR}}`$ = $``$15 – 1 km s-1) along the same line as Figure 4a. The CO $`J`$=1–0 spectra used to obtain the 3–2/1–0 ratio were reproduced from the mapped data obtained every 10′′ spacing by resampling them with the 14′′ Gaussian beam. If we assume that the density of the outflowing gas is almost constant along the lobe, $`n(\mathrm{H}_2)=2\times 10^4`$ cm-3, which is suggested from the ratios among the $`J`$=6–5, 3–2, and 1–0 lines at positions A, B, and C, the variation of the 3–2/1–0 ratio corresponds to that of the kinetic temperature (see Figure 3). At position A and the southernmost position of Figure 4b, the values of the 3–2/1–0 ratio contain large uncertainties because the CO $`J`$=1–0 blueshifted wings are too weak. If we exclude these two points at both ends, we can find that the variation of the 3–2/1–0 ratio, that is, the variation of the kinetic temperature correlates well with that of the velocity dispersion of the CO $`J`$=3–2 emission. Figure 4b shows that the temperature of the outflowing gas, which is $``$80 K near the central source, increases toward position B and reaches up to $`T_{\mathrm{kin}}`$ = 170 K (under the assumption of $`dv/dr`$ = 180 km s-1) at the position where the largest velocity dispersion is observed. Such a remarkable agreement between the velocity dispersion and kinetic temperature can be explained by the model of jet-driven bow shock (e.g., Masson & Chernin 1993), in which both momentum and energy input from the jet to the swept-up medium are maximal behind the bow shock. A model calculation including the balance of heating and cooling (Hatchell, Fuller, & Ladd 1999) has shown that the kinetic energy loss at the bow shock is sufficient to heat the swept-up molecular gas and produces the temperature gradient that increases from the protostar towards the bow shock. At the downstream side of position B, the temperature once drops to $``$50 K and increases again to $``$80 K at 10′′ upstream side of position C. The temperature increase toward the shock near position C is less significant as compared to that toward the shock at position B.

Toward the shock at position B, the emission from the ro-vibrational H2 lines in the near-infrared and that from the pure rotational lines of H2 in the mid-infrared have been detected (Davis & Eislöffel 1995; Cabrit et al. 1998). Since the required temperatures to excite the ro-vibrational and pure rotational lines of H2 are $``$2000 K and 700-800 K, respectively, the detection of these lines implies that the molecular gas in the shocked region is heated to above 2000 K. Because of the very short cooling timescale of such hot molecular gas, $``$10 yr (e.g., Neufeld, Lepp, & Melnick 1995), the presence of hot gas means that the energy input was made very recently in the shocked gas at position B. On the other hand, no H2 emission appears to be associated with the other shock near position C, suggesting that the major energy input process has been already suspended there. By using the total CO luminosity expected at position B, 1.2$`\times `$ 1031 erg s-1, and the H2 column density derived at the same position, $``$6$`\times `$1021 cm-2, we can estimate the cooling rate per H2 molecule by the rotational CO line emission to be $``$3$`\times `$10-25 erg s-1. If we assume that the rotational lines of CO are the major coolant of the gas with temperature less than several hundred K as in the case of the gas with $`T_{\mathrm{kin}}<`$ 40 K (Goldsmith & Langer 1978), it takes $``$1000 yr to cool the gas from the temperature at position B, 120 K, to 75 K derived at position C. Since this cooling time scale is comparable to the difference of time scale between the two shocks (Gueth et al. 1996), the lower gas kinetic temperature at the shock near position C can be explained if this shock is in a later evolutionary stage, in which the gas has been cooled through radiation of the CO rotational lines. Theoretical studies of the thermal balance in dense \[$`n`$(H2)$``$103 cm-3\] molecular clouds predicted that H2O and O i would be important gas coolants (e.g., Neufeld et al. 1995). However, because of the large dipole moment of the water molecule, the H2O lines cannot be the major coolant of the gas unless the density is higher than 106 cm-3. On the other hand, contribution of the O i fine structure lines to the cooling of the gas with $`T_{\mathrm{kin}}>`$ 100 K could be comparable to that of the CO rotational lines. However, the gas cooling rate through O i emission is expected to decrease rapidly with decreasing gas temperature, because the upper energy state of the lowest transition, $`{}_{}{}^{3}P_{1}^{}`$, is 228 K above the ground.

The total energy released by the rotational CO lines estimated toward position C is one order of magnitude lower than that of position B. Therefore, the cooling time scale at position C where the gas kinetic temperature is $``$80 K is estimated to be of the order to 104 yr, which is comparable to the dynamical time scale of the outflows driven by the Class 0 sources. This implies that the kinetic temperature of the outflowing gas can be maintained at several tens of K, which is a factor of 3–5 higher than the ambient temperature, during the life time of the Class 0 stage

We would like to thank the staff of the JCMT for the operation of our observations, and Dr. R. Phillips for assisting us with our observations and data reduction. We acknowledge T. Oka for providing the LVG calculation code of good performance, and T. Umemoto for providing us with the CO and 13CO $`J`$=1–0 data observed with the NRO 45m telescope. We also acknowledge the helpful comments of an anonymous referee whose suggestions improved the paper.


text:
# Novel continuum modeling of crystal surface evolution

\[

## Abstract

We propose a novel approach to continuum modeling of the dynamics of crystal surfaces. Our model follows the evolution of an ensemble of step configurations, which are consistent with the macroscopic surface profile. Contrary to the usual approach where the continuum limit is achieved when typical surface features consist of many steps, our continuum limit is approached when the number of step configurations of the ensemble is very large. The model can handle singular surface structures such as corners and facets. It has a clear computational advantage over discrete models.

\]

The behavior of classical physical systems is typically described in terms of equations of motion for discrete microscopic objects (e.g. atoms). In many cases, the behavior of such systems is smooth when observed on macroscopic length and time scales. It is useful to describe these systems in terms of continuum, coarse-grained models, which treat the dynamics of the macroscopic, smoothly varying, degrees of freedom rather than the microscopic ones. Such models are more amenable to analytical treatments and have enormous computational advantages over their discrete counterparts.

Many physical systems exhibit a macroscopically smooth behavior everywhere, except in small regions of space where their behavior is singular. Examples are dislocations in a crystal, cracks in crystalline material, facet edges on crystal surfaces, etc. These singular regions are of interest because they frequently drive the dynamics of the whole system. It is an interesting and important challenge to develop continuum descriptions of such systems, since standard phenomenological continuum models completely fail in the singular regions.

In this Letter we address the above problem in the context of the dynamics of crystal surfaces. Below the roughening temperature the evolution of these surfaces proceeds by the nucleation, flow and annihilation of atomic steps. These steps originate either from a miscut of the initial surface with respect to a high symmetry plane of the crystal or as the boundaries of islands and vacancies which nucleate on the surface during morphological evolution (we ignore screw dislocations as sources of steps). The steps are separated by terraces, which are parallel to a high symmetry crystal orientation. Each step consists of atomic kinks separated by straight portions along closed-packed orientations. However, on length scales larger than the typical distance between kinks the step is a smooth line. Thus the evolution of the surface corresponds to the motion of discrete smooth lines.

One can model step motion by solving the diffusion problem of adatoms on the terraces with boundary conditions at step edges. These conditions account for attachment and detachment of atoms to and from step edges. The local flux of adatoms at step edges and the resulting step velocities can then be calculated. This approach was introduced long ago by Burton, Cabrera and Frank , and was further developed by other authors . Given the arrangement of few neighboring steps, these models allow a calculation of the normal velocity of a step. The number $`M`$ of neighboring steps on each side required for this calculation depends on the step kinetics and on the interaction range between steps (in the examples below $`M=2`$). Step models are capable of describing surface evolution on the mesoscopic scale with significant success . Such models pose a serious challenge for numerical computations, and can be solved only for small systems.

Several attempts were made to construct continuum models of stepped surfaces , to study their large scale properties. The general idea behind these attempts is that step flow can be treated continuously in regions where every morphological surface feature is composed of many steps. If we label surface steps by the index $`n`$, the continuum limit in these models is obtained by taking $`n`$ to be continuous. The outcome of these attempts are partial differential equations for surface evolution. Such continuum models are fairly successful in describing the evolution of smooth surfaces of very simple morphology. However, below the roughening temperature, crystal surfaces have singularities in the form of corners and macroscopic facets. The latter are a manifestation of the cusp singularity of the surface free energy at high symmetry crystal orientations. The assumption that every surface feature is composed of many steps clearly breaks down on macroscopic facets where there are no steps at all. Thus, existing continuum models fail near singular regions.

A possible resolution of this problem is to solve the continuum model in non-singular regions and use matching conditions at the singular facet edges . This is a good solution in extremely simple cases (see, e.g. ), but generally it is not possible to derive the appropriate matching conditions from the microscopic models . As a result this approach is not very useful.

Another suggested solution is to round the surface free energy cusp . The resulting continuum models replace facets by relatively flat but analytic regions. As explained below, this approach does not correctly reproduce the evolution of facet edges. Thus, all existing continuum models fail in the description of the evolution of surfaces with singularities. This failure is both conceptual and practical, in the sense that it leads to qualitatively as well as quantitatively incorrect results.

We now propose a conceptually new definition of the continuum limit, which we term configurational continuum. This definition allows construction of continuum models, free of all the limitations of standard continuum models. Our key observation is that a continuous surface height profile can be represented by many similar, but not identical, step configurations. Figure 1 is a demonstration of this point for a one dimensional surface. It shows a continuous height function, $`h(x)`$, as a function of position $`x`$ (thick solid line), and three valid discrete representations of this profile as step configurations. In configurational continuum we define the continuous height profile as the upper envelope of the discrete height functions of an ensemble of many such step configurations. To construct these configurations, let $`a`$ be the height of a single step and $`N`$ the number of configurations in the ensemble. We construct the ensemble so that the height difference $`\delta h`$, between two adjacent configurations is $`a/N`$, as depicted in Fig. 1. The configurational continuum limit is obtained when $`N\mathrm{}`$. The generalization to higher dimensions is straightforward.

The dynamics of the continuum model is as follows. Each step configuration of the ensemble evolves according to the microscopic dynamics. As a result, the envelope of discrete height functions changes with time, thus defining the evolution of the continuous height function.

The crucial assumption hidden in this definition of the continuum model is that our construction leads to a mathematically well-defined height function at all times. When does this assumption hold? Consider two initially similar configurations of the ensemble. Our continuum limit is well-defined provided these two configurations have similar microscopic dynamics and hence remain similar at later times . This assumption has to hold in the standard continuum definition as well. If two initially similar configurations evolve very differently, one must follow the specific microscopic configuration of interest.

We now derive the evolution equation for the continuous height of a two dimensional surface, $`h(\stackrel{}{r},t)`$, at position $`\stackrel{}{r}`$ and time $`t`$. Here $`\stackrel{}{r}`$ is a vector in the $`xy`$ plane, defined as the high symmetry crystal plane parallel to the terraces. $`h`$ changes with time due to the flow of steps through $`\stackrel{}{r}`$, and due to nucleation and annihilation of steps. At this stage we disregard nucleation processes and include them later. First, we consider positions, which are not local extrema of the height profile. It is obvious from the construction of the configurational continuum, that for each point $`\stackrel{}{r}`$ there is exactly one configuration in the ensemble, which has a step that passes through $`\stackrel{}{r}`$ at the time $`t`$. That step lies along the unique equal-height contour line, which passes through $`\stackrel{}{r}`$. As is demonstrated in Fig. 1, the exact positions of neighboring steps in the same configuration can be calculated from the knowledge of $`h(\stackrel{}{r},t)`$, and the fact that in this configuration there is a step at $`\stackrel{}{r}`$. Hence, we can use the discrete step model to calculate the normal velocity of the step at position $`\stackrel{}{r}`$ at time $`t`$. Let us denote this velocity by $`\stackrel{}{v}_s(\stackrel{}{r},t)`$. Note that at different positions, $`\stackrel{}{v}_s`$ is the normal velocity of steps which belong to different configurations in the ensemble.

Next we define a two dimensional unit-vector $`\widehat{n}(\stackrel{}{r},t)`$ in the $`xy`$ plane, which is normal to the contour line at $`\stackrel{}{r}`$, and points in the direction opposite to $`\stackrel{}{v}_s(\stackrel{}{r},t)`$. This allows us to define the directional gradient

$$h_{\widehat{n}}(\stackrel{}{r},t)\widehat{n}(\stackrel{}{r},t)\underset{ϵ0_+}{lim}\frac{h(\stackrel{}{r}+ϵ\widehat{n},t)h(\stackrel{}{r},t)}{ϵ},$$
(1)
where the directionality is needed in order to deal with corners and facet edges. The gradient is in the direction from which steps flow towards the point $`\stackrel{}{r}`$. This is useful for the calculation of the current of steps arriving at $`\stackrel{}{r}`$:

$$J(\stackrel{}{r},t)=\frac{N}{a}\left|h_{\widehat{n}}(\stackrel{}{r},t)\stackrel{}{v}_s(\stackrel{}{r},t)\right|,$$
(2)
where we have used the fact that the local step density is $`\left|h_{\widehat{n}}(\stackrel{}{r},t)\right|N/a`$. Note that $`J`$ is the current of steps belonging to all configurations in the ensemble, and not to one particular configuration. Since each step (from any configuration), which passes through $`\stackrel{}{r}`$ changes the height of the ensemble envelope by $`a/N`$ (see Fig. 1), the continuous height profile obeys the evolution equation

$$\frac{h(\stackrel{}{r},t)}{t}=h_{\widehat{n}}(\stackrel{}{r},t)\stackrel{}{v}_s(\stackrel{}{r},t).$$
(3)
The above derivation of the evolution in the continuum is not valid at local extrema of the surface, because generally one cannot define a unique equal-height contour line which passes through such a point. To avoid the problem, we define $`h/t`$ at local extrema as the limit of the height time derivative as one approaches these points. This is justified, since there are no microscopic realizations of the surface with steps exactly at the local extrema.

At this point we emphasize that the configurational continuum evolution is formally identical to the evolution of the discrete step model. This statement is almost trivial, since by definition Eq. (3) follows the envelope defined by the ensemble of configurations, and each configuration evolves with step velocities calculated from the discrete step model. Thus Eq. (3) is exact. Similarly to other continuum models it is solved numerically by discretization of space, which is the only approximation involved in such solutions.

What is the relation between configurational continuum and other existing continuum models? We will show elsewhere that in regions where $`h(\stackrel{}{r},t)`$ is a smooth function of $`\stackrel{}{r}`$, the evolution equation (3) reduces to the differential equation of standard continuum models. However, near corners or facets (3) accounts for the discrete nature of steps and cannot be written as a differential equation, due to the singular nature of the profile. Contrary to other continuum models, configurational continuum is exact even at singular points and is capable of describing both smooth and singular behaviors without artificial smoothing or matching conditions. To the best of our knowledge, configurational continuum is the only continuum model capable of correctly treating the dynamics of singular lines .

The fact that no matching conditions are needed in configurational continuum is truly important for possible generalizations. A matching condition approach has been used in other areas, such as plasticity of crystals. The idea there is to solve continuum elasticity far from dislocations and a discrete model near dislocations. The inaccurate matching conditions between the different regions lead in this case to artificial reflection of phonons from the boundaries between the regions . A proper generalization of configurational continuum to atomistic systems may eliminate these artificial phonon reflections.

In terms of computational efficiency, our model is comparable to standard continuum models. The computational advantage over discrete models arises from the fact that in macroscopically smooth regions the density of steps is much larger than the needed density of discretization points. Some care should be taken near singular lines, and discretization there should be dense enough to allow an accurate interpolation of the singular profile. This, however, should not affect the overall efficiency, and when most of the surface is smooth the number of steps needed to be simulated is usually orders of magnitude smaller than the number of steps in the system.

We now apply the configurational continuum approach to a few simple cases in order to demonstrate its validity. To this end, we consider the evolution of a conic structure which consists of circular concentric steps. The adatom concentration on the terraces can be integrated out due to the radial symmetry of the system, yielding a set of coupled non-linear differential equations for the evolution of step radii . We have solved this discrete model numerically with and without growth flux. In the absence of flux, the cone decays by step flow. During the decay, a facet develops at the top of the cone, and its radius grows as $`t^{1/4}`$. With flux, the cone grows except at the peak which initially decays and then saturates. A facet forms at the peak after saturation.

The continuum model we solved for this example is a two dimensional model, which can, in principle, develop non radially symmetric morphology. The microscopic dynamics we used was a two dimensional generalization of the microscopic equations for step radii of the discrete one dimensional model. Any microscopic dynamics, such as a full solution of the diffusion equation on each terrace, can be used in the framework of configurational continuum. For the sake of demonstrating the validity of our approach the simple dynamics we chose is sufficient.

Figures 2(a) and (b) show a comparison between a numerical solution of the discrete step model with growth flux and a cross section of the two dimensional surface obtained from the solution of our continuum model at various times. The agreement is quite impressive. In particular, the width and height of the facet obtained from the continuum model are equal (within numerical accuracy) to those of the discrete model at all times. A similar agreement between the two models is achieved in the absence of growth flux. For comparison, we also investigated the dynamics of the decaying cone using the continuum model of . Within this model, the cusp in the free energy is rounded, and the evolution is studied with various degrees of rounding. We found that even in the limit of small rounding the height profile is very different from the profile of the discrete model. Although a facet does develop, its size is about 40% larger than the facet of the discrete model at all times. In addition, the shapes of the height profile differ qualitatively.

Island nucleation can also be included in our continuum model. To this end we add a microscopic nucleation scheme, which defines the nucleation probability at any given position on the surface, given a specific microscopic step configuration. Within our continuum approach, the nucleation probability at a point on the continuous surface is the ensemble average of the microscopic nucleation probabilities at this point. For demonstration purposes we included the simple scheme where the nucleation probability in any given step configuration is proportional to the square of the local adatom concentration. In fact, island nucleation is already included in the simulations presented in Fig. 2. Nucleation events occur on the top facet once it becomes large enough, but there is hardly any nucleation on the finite slope regions, because the terraces there are very narrow. Figure 2(c) shows an island on the top facet, which nucleated and started growing. Vacancy nucleation during evaporation can be included in a similar manner.

We now turn to the more demanding example of bidirectional sinusoidal grating relaxation. Here the initial surface height profile of wave length $`L`$ is given by $`h_L(\stackrel{}{r},t=0)=h_0\mathrm{sin}\left(2\pi x/L\right)\mathrm{sin}\left(2\pi y/L\right)`$. The relaxation of this profile towards a flat surface was studied by Rettori and Villain , who gave an approximate solution to a step flow model, in the limit where the interaction between steps can be neglected in comparison with the step line tension. The decay of bidirectional sinusoidal profiles was also studied numerically . We now apply our model to this problem assuming diffusion limited kinetics without deposition flux.

For weak interactions between steps, the surface height evolves according to $`h_L/t^2\kappa L^3`$, where $`\kappa `$ is the step curvature. We therefore expect the following scaling law: $`h_L(\stackrel{}{r},t)=h_{L=1}(\stackrel{}{r}/L,t/L^3)`$. Figure 3 shows the data collapse of cross sections of profiles from the configurational continuum model. The different symbols correspond to different wave lengths at time $`t=t_0L^3`$ for a fixed $`t_0`$. The quality of the data collapse shows that the scaling scenario holds very accurately. Note that large facets have developed at the surface extrema, and they are connected by very steep slopes. This shape does not agree with Rettori and Villain’s heuristic argument , which predicts that facets appear also near $`h=0`$ lines. Nevertheless, their prediction that the decay time of the height profile goes as $`L^3`$ is in agreement with both the scaling law and our numerical solutions. Our simulations show that for much stronger repulsive interactions between steps, the scaling law does not hold.