text:
# L-Band Photometry of L and T Dwarfs

## 1 Introduction

The discovery of hundreds of ultra low-mass dwarfs and substellar objects by the 2 Micron All Sky Survey (2MASS), the Deep Near Infrared Southern Sky Survey (DENIS) and the Sloan Digital Sky Survey (SDSS) has provided a large data base from which to select L and T dwarfs for further study. The L dwarfs identified by these surveys have been classified using far optical ($``$ 0.6-1.0 $`\mu \mathrm{m}`$) spectroscopy and most have published J (1.1-1.4 $`\mu \mathrm{m}`$), H (1.5-1.8 $`\mu \mathrm{m}`$) and K band (2.0-2.4 $`\mu \mathrm{m}`$) magnitudes (Delfosse et al., 1997; Kirkpatrick et al., 1999; Martín et al., 1999; Fan et al., 2000; Gizis et al., 2000; Kirkpatrick et al., 2000). In addition several groups have used near-infrared spectroscopy from 1 to 2.5 $`\mu \mathrm{m}`$ to classify and study L dwarfs (Delfosse et al., 1997, 1999; Tokunaga & Kobayashi, 1999; McLean et al., 2000; Leggett et al., 2000, 2001; Reid et al., 2001; Testi et al., 2001). Comparatively, little work has been done beyond 2.5 $`\mu \mathrm{m}`$, where to date only six L dwarfs have published L-band magnitudes (Jones et al., 1996; Leggett, Allard, & Hauschildt, 1998; Tokunaga & Kobayashi, 1999; Leggett et al., 2001). Since brown dwarfs emit most of their flux in the near infrared, observations beyond 2.5 $`\mu \mathrm{m}`$ are useful to the determination of the bolometric luminosity of these objects. In addition water ($`\mathrm{H}_2\mathrm{O}`$), carbon monoxide (CO), and methane ($`\mathrm{CH}_4`$) which are among the most abundant atmospheric species after hydrogen (H2) are important opacity sources at these wavelengths.

Indeed $`\mathrm{H}_2\mathrm{O}`$ and $`\mathrm{CH}_4`$ opacities dominate the L-band (2.5-4.1 $`\mu \mathrm{m}`$) radiative transfer. In particular the $`L_s`$ band (2.5-3.5 $`\mu \mathrm{m}`$) overlaps the strong $`\nu _3`$ fundamental band of $`\mathrm{CH}_4`$ at 3.3 $`\mu \mathrm{m}`$, which is observable in some of the late L dwarfs and all of the T dwarfs (Noll et al., 2000). Because chemical equilibrium favors $`\mathrm{CH}_4`$ over CO with falling atmospheric temperature (Fegley & Lodders, 1996), the abundance of $`\mathrm{CH}_4`$ helps to constrain the effective temperature ($`\mathrm{T}_{\mathrm{eff}}`$). In the L band (3.5-4.1 $`\mu \mathrm{m}`$) $`\mathrm{H}_2\mathrm{O}`$ alone is the dominant opacity source. Atmosphere models (Marley et al., 2001) suggest that the (K-L) color of L dwarfs is relatively insensitive to the details of cloud models, or even the absence of clouds, making this color a potentially interesting diagnostic of $`\mathrm{T}_{\mathrm{eff}}`$. When used together, the $`L_s`$ and L bands can probe for the appearance of $`\mathrm{CH}_4`$ in the latest L dwarfs, identify for trends in color which are a function of spectral type and temperature, and provide flux measurements for use in estimating the bolometric luminosity of these objects. Here we report L-band observations for a sample of L and T dwarfs, their observed trends, and comment on the L-dwarf $`\mathrm{T}_{\mathrm{eff}}`$ scale.

## 2 Observational Technique

The data were acquired using the Near Infrared Camera, NIRC (Matthews & Soifer, 1994) on Keck I at the W.M. Keck Observatory located in Mauna Kea, Hawaii on April 13th and 14th, UT 2000. Photometry was obtained using K, $`K_s`$, $`L_s`$, L and a 3.3 $`\mu \mathrm{m}`$ narrow band filter111The narrow band results will be presented elsewhere. for 18 L dwarfs and 5 T dwarfs which were observable at a low airmass, and chosen to span the L-dwarf spectral sequence. Where multiple objects of a given spectral type were available to observe, preference was given to those objects which were most typical in their (J-K) and (H-K) colors respective to spectral type (Kirkpatrick et al., 2000). The resulting data set includes objects with spectral types ranging from L0 to L8 in the Kirkpatrick et al. (1999) classification as well as five T dwarfs ranging from T0 to T6 in the A. Burgasser (2001, private communication) classification. Table 1 lists the targets, their coordinates, observed magnitudes and colors. Both nights were photometric with seeing ranging between 0.35″and 0.6″. Since the radiation background from the sky is considerably larger at 3 $`\mu \mathrm{m}`$ and beyond, integration times were limited to 0.012 and 0.03 seconds in the L and $`L_s`$ bands. To achieve these integration times, a 64 $`\times `$ 64 subsection of the chip was read out using the pat4xa timing pattern. Standard chopping and nodding procedures were also employed to minimize errors resulting from temporal variations in the sky radiation background.

To place our photometry onto a standard infrared photometric system, standard stars were chosen from the UKIRT Fundamental and Extended List (Hawarden et al., 2001) and were observed at the same airmass directly before or after the target objects. In the case where a UKIRT standard of the same airmass was not available to observe, standard stars were then chosen from the Persson et al. (1998) NICMOS standard list and their magnitudes were converted to the UKIRT system using the formula given by Hawarden et al. (2001). These lists provided standard K-band magnitudes, but not L or $`L_s`$ magnitudes. To determine a standard magnitude for our L-band observations we observed only main sequence standard stars with spectral types between A0V and G6V, and estimated a standard L or $`L_s`$ magnitude for each star using the intrinsic colors given by Bessell & Brett (1988) for a main sequence star of that spectral type. Since the (K-L) and (K-L) values of these stars vary between 0.0 for A0V to 0.05 for G6V, we estimate that this technique produces no more than a 5% error in our magnitude calculation and this error has been included in the final L and $`L_s`$ magnitude tabulation (Table 1). As an additional check standard stars that deviated largely from the mean L and $`L_s`$ zeropoint correction determined from all of the standard stars were discarded. This resulted in the elimination of no more than one standard star each night and we are confident that this technique produces valid L and $`L_s`$ magnitudes for our objects.

The standard stars used to calibrate the photometry in Table 1 are on the UKIRT IRCAM3 natural system. Convolutions of the filter transmission profiles for the NIRC and IRCAM3 filters with model spectra reveals a constant zeropoint offset between the two systems with no dependence on spectral type. Thus the photometry presented in this paper is on the UKIRT IRCAM3 system. It is important to note that this is not the same photometric system as the Mauna Kea Observatory near-infrared system (MKO-NIR). While the K band magnitudes for our L dwarfs can be placed on the MKO-NIR system using the transformation equations provided by (Hawarden et al., 2001), there is not a conversion equation for the L magnitudes. A comparison of the L filters for both the NIRC and MKO-NIR systems reveals that the MKO-NIR filter extends farther into the blue from 3.4 to 4.1 $`\mu \mathrm{m}`$ while the NIRC filter only extends from 3.5 to 4.1 $`\mu \mathrm{m}`$. This results in the MKO-NIR filter being more sensitive to the appearance of $`\mathrm{CH}_4`$ absorption at 3.3 $`\mu \mathrm{m}`$. The net effect is that for the hottest L dwarfs, our L magnitudes are almost identical to the MKO-NIR magnitudes. For the latest L dwarfs where $`\mathrm{CH}_4`$ is just beginning to affect the spectrum, our L magnitudes are 10% brighter, and for the T dwarfs where $`\mathrm{CH}_4`$ absorption is strongly established, our L magnitudes are 20% brighter (differences of 0.1 and 0.2 mag respectively). Thus any comparison of the L magnitudes from these two systems needs to take this offset into account.

As a final comment on the magnitudes presented here, we need to mention that although observations were acquired using both the K and $`K_s`$ filters, we found that in general the K-band magnitude differed from the $`K_s`$-band magnitude by no more than 0.05 for the L dwarfs and 0.1 for the T dwarfs. In addition since transformation equations from the UKIRT IRCAM3 system to the other standard infrared systems use the UKIRT K-band magnitudes for conversion and not $`K_s`$, we present only the K-band magnitudes.

## 3 Infrared Colors

Figure 1 presents the (K-L) color of our objects as a function of L- and T-dwarf spectral class. The closed squares represent the objects from Table 1, and the open circles represent (K-L) colors for six L dwarfs by Leggett et al. (2001). In general the (K-L) color gradually increases with later L- and T-dwarf spectral type. A linear least squares fit of the data to the L dwarf index shows that the (K-L) colors are reasonably well fit (correlation coefficient, $`r=0.80`$) by a line with a slope of 0.092 and an intercept point of 0.87. For L dwarfs the increase in (K-L) color is due to increasing molecular absorption in the K band. As the $`\mathrm{T}_{\mathrm{eff}}`$ decreases with later spectral type, the wings of the $`\mathrm{H}_2\mathrm{O}`$ band at 1.9 $`\mu \mathrm{m}`$ and pressure induced $`\mathrm{H}_2`$-$`\mathrm{H}_2`$ absorption together suppress the K-band flux. In the L band the principle opacity source is $`\mathrm{H}_2\mathrm{O}`$ alone. As the temperature decreases, the mean opacity increases more slowly in the L band than in the K band. The net result is that the total flux emitted in the K band falls faster than the total flux emitted in the L band with decreasing temperature. This produces the redder (K-L) colors with later L-dwarf spectral type.

A similar mechanism is responsible for the gradual reddening of the T dwarfs. For these objects $`\mathrm{CH}_4`$ joins $`\mathrm{H}_2\mathrm{O}`$ and $`\mathrm{H}_2`$-$`\mathrm{H}_2`$ as the principle opacity sources. By definition, the transition from the L to T dwarfs occurs when the $`\nu _2+\nu _3`$ combination band of $`\mathrm{CH}_4`$ first appears at 2.2 $`\mu \mathrm{m}`$ (Kirkpatrick et al., 1999). A linear least squares fit ($`r=0.97`$) of the T dwarf data produces a slightly steeper line with a slope of 0.161 and an intercept point of -0.19. The change in slope from the L to T dwarfs suggests that the additional appearance of $`\mathrm{CH}_4`$ opacity further suppresses the flux in the K band, resulting in a steeper (K-L) color increase throughout the T dwarfs. Since the abundance of $`\mathrm{CH}_4`$ varies with $`\mathrm{T}_{\mathrm{eff}}`$ and it is the dominant opacity source for T dwarfs at K band, but not L, the (K-L) color should prove to be a good diagnostic for estimating the $`\mathrm{T}_{\mathrm{eff}}`$ of T dwarfs.

Figure 2 presents the (K-$`L_s`$) color of our objects as a function of L- and T-dwarf spectral class. Like the (K-L) color, the increase in (K-$`L_s`$) through the L-dwarf spectral sequence is due to increasing $`\mathrm{H}_2\mathrm{O}`$ and $`\mathrm{H}_2`$-$`\mathrm{H}_2`$ absorption in the K band. A linear least squares fit applied to this data ($`r=0.50`$) produces a line with a slope of 0.06 and an intercept point of 0.14. It is interesting to note that there appears to be substantially more scatter in the (K-$`L_s`$) color for the later L dwarfs (and early T dwarfs) as compared to the early L dwarfs. This scatter may arise from differences in cloud cover (Gelino et al., 2001), or possibly atmospheric dynamics or metallicity producing greater variability in $`\mathrm{CH}_4`$ abundance in these objects.

For the T dwarfs $`\mathrm{CH}_4`$ absorption definitely occurs both at 3.3 $`\mu \mathrm{m}`$ and 2.2 $`\mu \mathrm{m}`$ which results in a suppression of the flux in the $`L_s`$ and K bands. Since the $`\nu _3`$ fundamental band of $`\mathrm{CH}_4`$ at 3.3 $`\mu \mathrm{m}`$ is a factor of 100 times stronger than the $`\nu _2+\nu _3`$ combination band at 2.2 $`\mu \mathrm{m}`$, there is a greater suppression of the flux in the $`L_s`$ band than in the K band which results in a bluer color for the T dwarfs. If $`\mathrm{CH}_4`$ absorption were to gradually appear in the latest L dwarfs with falling $`\mathrm{T}_{\mathrm{eff}}`$, we would expect the (K-$`L_s`$) color of the late L dwarfs to gradually transition into the bluer color of the early T dwarfs. However, the abrupt transition in (K-$`L_s`$) from the latest L dwarfs to the earliest T dwarfs suggests that perhaps very early T dwarfs, or very late L dwarfs have yet to be discovered, or that substantial $`L_s`$ band $`\mathrm{CH}_4`$ absorption appears rapidly over a very small $`\mathrm{T}_{\mathrm{eff}}`$ range.

Note that the two L1 dwarfs both have very blue (K-$`L_s`$) colors, but unexceptional (K-L) colors. These two objects were also unusually faint in narrow band 3.3 $`\mu \mathrm{m}`$ (NIRC pahl filter) images obtained the same night. This suggests that another opacity source is affecting the flux in the $`L_s`$ band for these two objects. Clearly 3 to 4 $`\mu \mathrm{m}`$ spectroscopy is needed to constrain what is suppressing the flux in the $`L_s`$ band. We are confident that the colors of the L1 dwarfs are not due to photometric errors since these objects were observed at different times during the night and with different standard stars. One of the L1 dwarfs, DENIS J1441-0945, is a binary system and if the two components of the system have different $`L_s`$ magnitudes this may be the cause of the color discrepancy. Indeed it does appear that the pair is of unequal brightness although they are too close for individual photometry222Magnitudes for this object in Table 1 are for the combined pair.. However this would imply that the other L1 dwarf is also a binary system and we see no evidence of multiplicity in our 0.5″images of 2MASSW J1035+2507.

Finally Figure 3 presents a (K-$`L_s`$) verses (K-L) color-color diagram with the data from Table 1. The progressive reddening in (K-L) is apparent as is the greater variability in (K-$`L_s`$) at the L to T transition. Note that the L7 object, 2MASSI J1526+2043, has the same (K-L) and (K-$`L_s`$) colors as an L4 dwarf. This L7 also has 2MASS (J-$`K_s`$), (J-H) and (H -$`K_s`$) colors which are most similar to an L3 or L4 dwarf (Kirkpatrick et al., 2000). It is possible that this L7 dwarf has been misclassified and this should be kept in mind when studying its location in Figures 1, 2 and 3.

## 4 Comparison to Models

A complete analysis of the L-dwarf $`\mathrm{T}_{\mathrm{eff}}`$ scale requires consideration of infrared and optical spectra as well as photometry. However, we can use the results presented here to draw some tentative conclusions by comparing to atmospheric models of Marley et al. (2001). These models account for sedimentation of atmospheric condensates as described in Ackerman & Marley (2001). Marley et al. (2001) found that the models with the precipitation efficiency factor $`f_{\mathrm{rain}}=5`$ best fit the (J-K) colors of the L dwarfs and we use those models here. The models suggest that the (K-L) color of dusty L-dwarf models is insensitive to details of cloud structure or gravity. For 1300 K $`\mathrm{T}_{\mathrm{eff}}`$ 2100 K a remarkably accurate linear expression ($`r=0.99`$) reproduces the (K-L) model colors as a function of $`\mathrm{T}_{\mathrm{eff}}`$. Equating this expression with the regression line derived in Figure 1 for the (K-L) data gives the following estimate of L dwarf $`\mathrm{T}_{\mathrm{eff}}`$ as a function of Kirkpatrick et al. (1999) spectral type $`\mathrm{L}_\mathrm{K}`$:

$$\mathrm{T}_{\mathrm{eff}}=2220100\mathrm{L}_\mathrm{K}$$
$`(1)`$
where $`\mathrm{T}_{\mathrm{eff}}`$ is the effective temperature in Kelvin, and $`\mathrm{L}_\mathrm{K}`$ = 0 to 8. This equation was fit using only models appropriate for $`M30\mathrm{M}_\mathrm{J}`$, however the colors depend only slightly on mass for this temperature range, and introduce only a few hundreths variation in the color difference at a fixed $`\mathrm{T}_{\mathrm{eff}}`$.

This equation is presented for its heuristic value as an indicator of the $`\mathrm{T}_{\mathrm{eff}}`$ range of L-dwarfs and not as a definitive analysis. Nevertheless it is interesting to make some comparisons. The earliest L dwarfs appear to have $`\mathrm{T}_{\mathrm{eff}}2100`$ to $`2200\mathrm{K}`$ and the L-dwarf regime extends down to about $`\mathrm{T}_{\mathrm{eff}}1400\mathrm{K}`$. On average each Kirkpatrick et al. (1999) spectral type spans about 100 K in $`\mathrm{T}_{\mathrm{eff}}`$.

The $`\mathrm{T}_{\mathrm{eff}}`$ derived here for L dwarfs can be compared to the results predicted by other research groups. Leggett et al. (2001) compares low resolution spectra and colors of a selection of L dwarfs to models of Allard, Hauschildt, & Schweitzer (2000) to derive $`\mathrm{T}_{\mathrm{eff}}`$ for several objects. Their spectral synthesis fitting gave an estimated $`\mathrm{T}_{\mathrm{eff}}`$ of 1800-1900 K for the entire L dwarf sequence (L0 to L8). Comparison of observed fluxes to evolution models by Chabrier et al. (2000) resulted in $`\mathrm{T}_{\mathrm{eff}}`$ ranging from 2050 - 2350 K for L0 to 1400 to 1600 K for L7. Basri et al. (2000) proposed that $`\mathrm{T}_{\mathrm{eff}}`$ on the Martín et al. (1999) scale varies from about 2200 K (L0) to 1600 K (L6). The above $`\mathrm{T}_{\mathrm{eff}}`$ estimates are all based on fits to the “cleared dust” models of Allard et al. (2000) which do not include dust in the radiative transfer. Since the Marley et al. (2001) models do include dust opacity, accounting for sedimentation, the atmosphere models are warmer at fixed $`\mathrm{T}_{\mathrm{eff}}`$ than comparable cleared-dust models by several hundred degrees. Thus the $`\mathrm{T}_{\mathrm{eff}}`$ range implied by Eq. 1 is cooler and may be more reasonable.

Our estimated effective temperatures are somewhat higher, particularly for the later-type objects, then those of Pavlenko, Zapatero Osorio, & Rebolo (2000) who derive temperatures of 2200 K for M9.5 to 1200 K for L7 based on spectral synthesis fitting of far optical spectra using dusty atmosphere models. We note that their fitted $`\mathrm{T}_{\mathrm{eff}}`$ depends sensitively upon the treatment of dust and alkali gas opacity as well as the inclusion of an unknown “additional opacity”. In contrast the models of Marley et al. (2001) fit both the optical and near-infrared colors of L dwarfs by including a more physically-based description of condensate sedimentation and properly accounting for chemical equilibrium with cloud formation. Both dust sedimentation and the chemical equilibrium treatment of alkali metal opacity become more important with later L spectral type (Marley et al., 2001). These effects are likely responsible for the discrepancies between our results.

## 5 Conclusions

We find that the (K-L) and (K-$`L_s`$) colors become redder with later L-dwarf spectral type. Similarly the (K-L) color of the T dwarfs steadily reddens by a full magnitude over the range of objects observed, and this color may prove useful in the future classification and determination of effective temperature in T dwarfs. We find that (K-$`L_s`$) for the early T dwarfs is bluer than (K-$`L_s`$) for the latest L dwarfs due to $`\mathrm{CH}_4`$ absorption but that a gradual transition between the two spectral classes is not seen. This suggests that perhaps very early T dwarfs or very late L dwarfs have yet to be found or that $`\mathrm{CH}_4`$ absorption appears rapidly over a very small effective temperature range. By using the observed (K-L) colors and theoretical models which include dust opacity and account for cloud sedimentation (Marley et al., 2001), we estimate an $`\mathrm{T}_{\mathrm{eff}}`$ range of 2200 K to 1400 K for the entire L-dwarf spectral sequence (L0 to L8).

We thank Sandy Leggett for helpful comments and advice. D.S. and M.M. acknowledge support from NASA grants NAG2-6007 and NAG5-8919 and NSF grants AST-9624878 and AST-0086288. The data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. We would also like to thank instrument specialists Bob Goodrich and Randy Campbell for their support and assistance in using the NIRC instrument.


format: text

text:
# A Simple Model of Liquid-Liquid Phase Transitions

\[

## Abstract

In recent years, a second fluid-fluid phase transition has been reported in several materials at pressures far above the usual liquid-gas phase transition. In this paper, we introduce a new model of this behavior based on the Lennard-Jones interaction with a modification to mimic the different kinds of short-range orientational order in complex materials. We have done Monte Carlo studies of this model that clearly demonstrate the existence of a second first-order fluid-fluid phase transition between high- and low-density liquid phases.

PACS number: 64.70.Ja

\]

The most common example of a first-order phase transition is that between a liquid and a gas, such as boiling water. On the other hand, while many transitions between different solid phases of homogeneous materials are also well known, it is only relatively recently that evidence of a second fluid-fluid phase transition has been found. In fact, liquid-liquid phase transitions (LLPT) have been suggested in liquid S, Ga, Se, Te, I2, Cs, and Bi .

Stell and Hemmer  showed the existence of LLPT in a one dimensional model with a softened hard core potential and a long-range negative attraction. Their work was later studied in more detail by Franzese, Malescio, et. al.  and Sadr-Lahijany, Scala, et. al. .

Mitus, Patashinskii and Shumilo  proposed LLPT in molten salt at high pressure based on a phenomenological model. Ferraz and March suggested a similar LLPT in carbon, with indirect experimental evidence being found by Togaya . Glosli and Ree  published results of a first-order liquid-liquid phase transition in molten Carbon between two thermodynamically stable liquid phases. Extensive computer simulations on models of water have supported the existence of a LLPT in the metastable region . Experimental results supporting the evidence of liquid-liquid phase transitions in water have also been found . Katayama and Mizutani et. al.  found a liquid-liquid phase transition in molten phosphorus using x-ray diffraction.

Our objective is to investigate the general phenomenon instead of studying LLPT for any particular substance. We have developed a simple model that exhibits a transition between high- and low-density liquids at high pressure.

The behavior of our model is constructed to be similar to behavior seen in simulations of real substances such as water and carbon but without introducing the complexity of having to simulate molecular orientations as in water and carbon.

To mimic the effects of local ordering, we have represented the different relative local orientations of the molecules with a spin-one-half variable. The interactions between particles with the same spin are given by the original Lennard-Jones expression, while the interactions between particles with opposite spin are purely repulsive.

$`\varphi _{}(r)`$ $`=\varphi _{}(r)`$ $`=4ϵ({\displaystyle \frac{\sigma _l^{12}}{r^{12}}}{\displaystyle \frac{\sigma _l^6}{r^6}})`$ (1)
$`\varphi _{}(r)`$ $`=\varphi _{}(r)`$ $`=4ϵ({\displaystyle \frac{\sigma _u^{12}}{r^{12}}})`$ (2)
Note that there are different values of $`\sigma `$ for like and unlike spins. This is an important feature of the model, and some properties, including the symmetries of the solid phases, are sensitive to the relative values of $`\sigma _l`$ and $`\sigma _u`$. The LLPT occurs when $`\sigma _u`$ is smaller than $`\sigma _l`$, so that by re-orienting the spins, the particles are capable of forming different co-ordination numbers and local structures. We have performed most of our investigations for the case in which the ratio is $`1/2`$.

Our model can be modified to couple to a fictitious external magnetic field. If we define the magnetic moment as the sum of all spins, then the Hamiltonian is given by,

$$H=\frac{1}{2}\underset{ij}{}\varphi _{ij}(r_{ij})h\underset{i}{}\sigma _i$$
(3)
$`h`$ is the fictitious magnetic field which is set to zero in our simulations and $`\sigma _i=\pm 1`$ is the direction of the spins of individual particles.

We used the units $`T^{}=k_BT/ϵ`$, $`P^{}=\sigma _l^2k_BP/ϵ`$. We have set the potential cutoffs at 3$`\sigma _l`$.

We performed Monte Carlo simulations of the two-dimensional version of our model in various ensembles. Although our simulations are done in two dimensions, our model is not limited to two dimensions. The simulations turned out to be rather difficult, and it was necessary to extend the usual techniques to improve efficiency. However, we did find clear evidence of an LLPT at high pressures. We are able to map out both PT and P$`\rho `$ phase diagrams.

For our simulations of fluids, we confined the particles to a square with periodic boundary conditions. For those simulations that were extended to include solid phases, we used parallelograms to allow us to vary the angle of the boundary conditions, as well as the volume of the container. The optimization of the Monte Carlo simulation step sizes for individual particle motion, changes of the volume of the box, and changes of the angle of the parallelogram were dynamically optimized using the Acceptance Ratio Method . Metropolis flips to maintain equilibrium for the spins associated with each particle were also carried out.

Volume changes on the dense fluid turned out to be rather inefficient because $`\varphi _{}`$ increases very rapidly at short distances. We solved this problem by introducing a cluster Monte Carlo move. The clusters are formed by creating bonds between particles with probability,

$$P(r)=\{\begin{array}{cc}1\hfill & r<r_{min}\hfill \\ (r_{max}r)/(r_{max}r_{min})\hfill & r[r_{min},r_{max}]\hfill \\ 0\hfill & \text{otheriwse}\hfill \end{array}$$
(4)
and then changing the size of the box by rescaling the locations of the centers of mass of the clusters. This cluster move is extremely effective, with improvements in the acceptance ratio by factors of up to $`10^{10}`$.

The probability of accepting a proposed cluster move is given as follows. Let $``$ be the set of all edges which joins any two particles which are less than 3$`\sigma _l`$ apart, and let $``$ be the subset of $``$ consisting of bonds formed with the probability P(r), where $`r`$ is the edge length joining the two particles. The probability of forming $``$ is,

$$P()=\underset{i}{}P(r_i)\underset{j}{}1P(r_j)$$
(5)
Bond configurations need to remain invariant in order to satisfy detailed balance. Therefore $`=^{}`$, $`r_i=r_i^{}`$ for all $`r_i`$ and $`r_i^{}^{}`$, but $`r_jr_j^{}`$ for $`r_j`$ and $`r_j^{}^{}^{}`$. Therefore,

$$\underset{i}{}P(r_i)/\underset{i^{}^{}}{}P(r_i^{})=1$$
(6)
and,

$$\frac{P(^{})}{P()}=\underset{j^{}^{}^{}}{}(1P(r_j^{}))/\underset{j}{}(1P(r_j))$$
(7)
Hence the acceptance probability is,

$`A(\varpi \varpi ^{})=`$ (9)
$`\text{min}(1,e^{\beta (\mathrm{\Delta }E+p\mathrm{\Delta }V)}{\displaystyle \frac{V^N}{V^N}}{\displaystyle \frac{_i^{}^{}[1P(r_i)]}{_i[1P(r_i)]}})`$
Another new move that has proven extremely effective is to form clusters of nearby particles that are less than $`1/\sqrt{\rho sin(\pi /3)}`$ apart and attempt to flip the spins of all particles in a cluster.

For canonical ensemble simulations using two boxes in equilibrium with each other, the total volume is conserved . An additional Monte Carlo move is introduced to transfer particles between the boxes.

Thermodynamic quantities were calculated from simulations in the constant pressure ensemble over relatively large pressure and temperature domains using the Multiple Histogram Method . The coexistence curve is mapped out by tracing the ridge line of the isothermal compressibility. The position of the tricritical point was estimated using finite size scaling. The P$`\rho `$ diagram was determined by combining other data with the results of simulations using the canonical ensemble .

The phase diagram for our model is shown in figure (1). At low pressures its behavior is virtually identical to that of the usual Lennard-Jones model. At higher pressures, the coexistence curve for the LLPT is shown projecting into the fluid region of the phase diagram. This line of first-order phase transitions actually ends in a tricritical point, which separates it from a line of critical points associated with the ordering of the “spins” used to define the model. This line joins the liquid-gas and liquid-liquid first-order coexistence curves at the tricritical points. However, it is only our access to the “magnetic” degrees of freedom in the model that make it obvious that the LLPT curves end in tricritical rather than critical points. If we did not have access to the magnetic degrees of freedom, rather careful measurement of the exponents characterizing the divergences at the end of the line of first-order transitions would be necessary to distinguish the two cases. It is not completely clear whether a fictitious ordering field in a more elaborate model would not have the same consequence of producing tricritical points.

The apparent discontinuity of the coexistence curves between the high-density solid, the low-density solid and the low-density liquid is an artifact of the accuracy of our simulations. These lines are obtained from the Multiple Histogram extrapolations in opposite directions. The fact that these two lines do not meet exactly, but terminated very close to each other shows the accuracy of our simulations.

To determine the liquid-liquid coexistence curve in the PT diagram, simulations were done with 160 particles at constant pressure, using larger systems to check the results. Sets of histograms were collected at $`T^{}`$=0.53, 0.55, and 0.565 over a range of pressures covering across the coexistence region. The system was equilibrated for 30,000 MCS/P before beginning to take data for another 300,000 MCS/P.

To determine the tricritical point for the liquid-liquid phase transition, additional simulations for 160, 240, 320, and 480 particles were done at constant pressures at $`T^{}`$=0.55 and 0.565. The Monte Carlo simulations used a total of 300,000 MCS/P for the 160 particle systems and typical equilibration runs of 30,000 MCS/P. For larger systems, longer runs were used. For the 480 particle system, a total of 1,200,000 MCS/P were used with typical equilibration runs of 100,000 MCS/P.

To determine the liquid-solid coexistence curve, simulations were performed at $`P^{}`$=0.5, 0.6, 0.78 and 0.9 over a range of temperatures covering the solid-liquid coexistence region. The system was equilibrated for 60,000 steps per particle before collecting data for another 600,000 steps per particle.

The phase transition between high- and low-density phases is strongly first-order. Figure (3a,b) shows snap shots of high- and low-density liquids near the coexistence curve for a system with 240 particles. The low-density liquid is dominated by parallel nearest-neighbor interactions with a stronger core repulsion and with co-ordination number six. The high-density liquid is dominated by an anti-parallel nearest-neighbor interaction that allows the particles to come closer together and with co-ordination number three. The differences in local orderings for low- and high-density liquids stabilize the liquids and create the possibility of liquid-liquid phase transitions.

The volume-energy probability density function (figure (4)) has a saddle point that is typical of a first-order transition. The coexistence curve for the LLPT has a negative slope, as shown in the PT diagram in figure (1), reflecting the higher entropy of the high-density liquid. The sharp bend of the liquid-solid coexistence curve at the triple point found from the multiple-histogram analysis is consistent with our observation of the liquid-liquid phase transition.

To find the tricritical point, we use the fact that peak values of the isothermal compressibility grow as $`𝒪(N)`$ at the coexistence curve and as $`𝒪(N^{0.925})`$  at the tricritical point, assuming that the transition is in the expected two-dimensional Ising class. Simulations were done with N=160, 240, 320 and 480 near the tricritical point. Peak values of the isothermal compressibility are then plotted against system sizes for various temperatures. Figure (6) shows the size dependence of the isothermal compressibility. Although it is hard to determine the exact location of the tricritical point, the data suggests that it is located at $`T^{}=0.56\pm 0.01`$.

Figure 2 shows the growing divergence of the magnetic susceptibility for 160, 240 and 320 particles. These magnetic susceptibility are plotted at $`P^{}=0.5`$ which is far away from the liquid-liquid coexistence curve and the gas-liquid coexistence curve.

The $`\rho `$-T diagram (figure (5)) was mapped out using simulations in the canonical ensemble using dual boxes to simulate coexistence without the inconvenience of an interface between the phases . The low-density liquid-gas coexistence region is only slightly lower than that of pure Lennard-Jones  due to the small effect of the differences in the models at low densities.

Figure(7a) shows a snapshot of the high-density crystalline state, which has three-fold symmetry and zero magnetic moment. The low-density crystalline state (figure(7b)) is hexagonal close packed with uniform spin. At zero temperature, we can calculate the location of the transition between the two solid phases to arbitrary precision. Its location is determined to be $`P^{}=0.534819\pm 10^7`$, with $`\rho ^{}=1.82646`$ and $`\rho ^{}=0.943122`$ for the high- and low-density solids respectively. For comparison, $`\rho ^{}=0.934721`$ at $`T^{}=0`$ and $`P^{}=0`$, so the low-density solid phase changes its density very little up to the boundary of the high-density phase.

It was found that symmetry of the solid phase is highly dependent on the ratio $`\tau =\sigma _u/\sigma _l`$. Figure (8) shows a solid with four-fold symmetry obtained from a simulation with $`\tau =0.6`$. We believe that the full range of solid phases is quite rich for this model.

We have developed a relatively simple model that demonstrates a liquid-liquid phase transition between high- and low-density phases. Although our simulations are two dimensional, our model is not restricted to two dimensions and we see no reasion to believe that a three dimensional simulation with our model would produce significantly different behavior. By comparing the behavior of our model with the properties of real systems, we hope to learn which properties are generic and which depend on details of a particular material. One immediate point of interest is the negative slope of the liquid-liquid coexistence line in our model, which reflects the high entropy in the high-density phase. This feature is, indeed, found in most materials that exhibit an LLPT. However, this does not appear to be universal, since the coexistence curve of molten carbon reported by Glosli and Ree  has a positive slope. At the present time, this difference in materials properties is not understood.

We would like to thank Prof. Robert Griffiths for his suggestions and helpful comments. We also would like to acknowledge the support from the Pittsburgh Supercomputing Center.


text:
# Gaussian Process Regression with Mismatched Models

## 1 Introduction

There has in the last few years been a good deal of excitement about the use of Gaussian processes (GPs) as an alternative to feedforward networks . GPs make prior assumptions about the problem to be learned very transparent, and even though they are non-parametric models, inference—at least in the case of regression considered below—is relatively straightforward. One crucial question for applications is then how ‘fast’ GPs learn, i.e., how many training examples are needed to achieve a certain level of generalization performance. The typical (as opposed to worst case) behaviour is captured in the learning curve, which gives the average generalization error $`ϵ`$ as a function of the number of training examples $`n`$. Good bounds and approximations for $`ϵ(n)`$ are now available , but these are all restricted to the case where the ‘student’ model exactly matches the true ‘teacher’ generating the data. In practice, such a match is unlikely, and so it is important to understand how GPs learn if there is some model mismatch. This is the aim of this paper.

In its simplest form, the regression problem is this: We are trying to learn a function $`\theta _{}`$ which maps inputs $`x`$ (real-valued vectors) to (real-valued scalar) outputs $`\theta _{}(x)`$. We are given a set of training data $`D`$, consisting of $`n`$ input-output pairs $`(x^l,y^l)`$; the training outputs $`y^l`$ may differ from the ‘clean’ teacher outputs $`\theta _{}(x^l)`$ due to corruption by noise. Given a test input $`x`$, we are then asked to come up with a prediction $`\widehat{\theta }(x)`$, plus error bar, for the corresponding output $`\theta (x)`$. In a Bayesian setting, we do this by specifying a prior $`P(\theta )`$ over hypothesis functions, and a likelihood $`P(D|\theta )`$ with which each $`\theta `$ could have generated the training data; from this we deduce the posterior distribution $`P(\theta |D)P(D|\theta )P(\theta )`$. For a GP, the prior is defined directly over input-output functions $`\theta `$; this is simpler than for a Bayesian feedforward net since no weights are involved which would have to be integrated out. Any $`\theta `$ is uniquely determined by its output values $`\theta (x)`$ for all $`x`$ from the input domain, and for a GP, these are assumed to have a joint Gaussian distribution (hence the name). If we set the means to zero (as is commonly done), this distribution is fully specified by the covariance function $`\theta (x)\theta (x^{})_\theta =C(x,x^{})`$. The latter transparently encodes prior assumptions about the function to be learned. Smoothness, for example, is controlled by the behaviour of $`C(x,x^{})`$ for $`x^{}x`$: The Ornstein-Uhlenbeck (OU) covariance function $`C(x,x^{})=\mathrm{exp}(|xx^{}|/l)`$ produces very rough (non-differentiable) functions, while functions sampled from the radial basis function (RBF) prior with $`C(x,x^{})=\mathrm{exp}[|xx^{}|^2/(2l^2)]`$ are infinitely differentiable. Here $`l`$ is a lengthscale parameter, corresponding directly to the distance in input space over which we expect significant variation in the function values.

There are good reviews on how inference with GPs works , so I only give a brief summary here. The student assumes that outputs $`y`$ are generated from the ‘clean’ values of a hypothesis function $`\theta (x)`$ by adding Gaussian noise of $`x`$-independent variance $`\sigma ^2`$. The joint distribution of a set of training outputs $`\{y^l\}`$ and the function values $`\theta (x)`$ is then also Gaussian, with covariances given (under the student model) by

$$y^ly^m=C(x^l,x^m)+\sigma ^2\delta _{lm}=(𝐊)_{lm},y^l\theta (x)=C(x^l,x)=(𝐤(x))_l$$

Here I have defined an $`n\times n`$ matrix $`𝐊`$ and an $`x`$-dependent $`n`$-component vector $`𝐤(x)`$. The posterior distribution $`P(\theta |D)`$ is then obtained by conditioning on the $`\{y^l\}`$; it is again Gaussian and has mean and variance

$`\theta (x)_{\theta |D}\widehat{\theta }(x|D)`$ $`=`$ $`𝐤(x)^\mathrm{T}𝐊^1𝐲`$ (1)
$`(\theta (x)\widehat{\theta }(x))^2_{\theta |D}`$ $`=`$ $`C(x,x)𝐤(x)^\mathrm{T}𝐊^1𝐤(x)`$ (2)
From the student’s point of view, this solves the inference problem: The best prediction for $`\theta (x)`$ on the basis of the data $`D`$ is $`\widehat{\theta }(x|D)`$, with a (squared) error bar given by (2). The squared deviation between the prediction and the teacher is $`[\widehat{\theta }(x|D)\theta _{}(x)]^2`$; the generalization error (which, as a function of $`n`$, defines the learning curve) is obtained by averaging this over the posterior distribution of teachers, all datasets, and the test input $`x`$:

$$ϵ=[\widehat{\theta }(x|D)\theta _{}(x)]^2_{\theta _{}|D}_D_x$$
(3)
Now of course the student does not know the true posterior of the teacher; to estimate $`ϵ`$, she must assume that it is identical to the student posterior, giving from (2)

$$\widehat{ϵ}=[\widehat{\theta }(x|D)\theta (x)]^2_{\theta |D}_D_x=C(x,x)𝐤(x)^\mathrm{T}𝐊^1𝐤(x)_{\{x^l\}}_x$$
(4)
where in the last expression I have replaced the average over $`D`$ by one over the training inputs since the outputs no longer appear. If the student model matches the true teacher model, $`ϵ`$ and $`\widehat{ϵ}`$ coincide and give the Bayes error, i.e. the best achievable (average) generalization performance for the given teacher.

I assume in what follows that the teacher is also a GP, but with a possibly different covariance function $`C_{}(x,x^{})`$ and noise level $`\sigma _{}^2`$. This allows eq. (3) for $`ϵ`$ to be simplified, since by exact analogy with the argument for the student posterior

$$\theta _{}(x)_{\theta _{}|D}=𝐤_{}(x)^\mathrm{T}𝐊_{}^1𝐲,\theta _{}^2(x)_{\theta _{}|D}=\theta _{}(x)_{\theta _{}|D}^2+C_{}(x,x)𝐤_{}(x)^\mathrm{T}𝐊_{}^1𝐤_{}(x)$$

and thus, abbreviating $`𝐚(x)=𝐊^1𝐤(x)𝐊_{}^1𝐤_{}(x)`$,

$$ϵ=𝐚(x)^\mathrm{T}\mathrm{𝐲𝐲}^\mathrm{T}𝐚(x)+C_{}(x,x)𝐤_{}(x)^\mathrm{T}𝐊_{}^1𝐤_{}(x)_D_x$$

Conditional on the training inputs, the training outputs have a Gaussian distribution given by the true (teacher) model; hence $`\mathrm{𝐲𝐲}^\mathrm{T}_{\{y^l\}|\{x^l\}}=𝐊_{}`$, giving

$$ϵ=C_{}(x,x)2𝐤_{}(x)^\mathrm{T}𝐊^1𝐤(x)+𝐤(x)^\mathrm{T}𝐊^1𝐊_{}𝐊^1𝐤(x)_{\{x^l\}}_x$$
(5)
## 2 Calculating the learning curves

An exact calculation of the learning curve $`ϵ(n)`$ is difficult because of the joint average in (5) over the training inputs $`X`$ and the test input $`x`$. A more convenient starting point is obtained if (using Mercer’s theorem) we decompose the covariance function into its eigenfunctions $`\varphi _i(x)`$ and eigenvalues $`\mathrm{\Lambda }_i`$, defined w.r.t. the input distribution so that $`C(x,x^{})\varphi _i(x^{})_x^{}=\mathrm{\Lambda }_i\varphi _i(x)`$ with the corresponding normalization $`\varphi _i(x)\varphi _j(x)_x=\delta _{ij}`$. Then

$$C(x,x^{})=\underset{i=1}{\overset{\mathrm{}}{}}\mathrm{\Lambda }_i\varphi _i(x)\varphi _i(x^{}),\text{and similarly}C_{}(x,x^{})=\underset{i=1}{\overset{\mathrm{}}{}}\mathrm{\Lambda }_i^{}\varphi _i(x)\varphi _i(x^{})$$
(6)
For simplicity I assume here that the student and teacher covariance functions have the same eigenfunctions (but different eigenvalues). This is not as restrictive as it may seem; several examples are given below. The averages over the test input $`x`$ in (5) are now easily carried out: E.g. for the last term we need

$$(𝐤(x)𝐤(x)^\mathrm{T})_{lm}_x=\underset{ij}{}\mathrm{\Lambda }_i\mathrm{\Lambda }_j\varphi _i(x^l)\varphi _i(x)\varphi _j(x)_x\varphi _j(x^m)=\underset{i}{}\mathrm{\Lambda }_i^2\varphi _i(x^l)\varphi _i(x^m)$$

Introducing the diagonal eigenvalue matrix $`(𝚲)_{ij}=\mathrm{\Lambda }_i\delta _{ij}`$ and the ‘design matrix’ $`(\mathrm{\Phi })_{li}=\varphi _i(x^l)`$, this reads $`𝐤(x)𝐤(x)^\mathrm{T}_x=\mathrm{\Phi }𝚲^2\mathrm{\Phi }^\mathrm{T}`$. Similarly, for the second term in (5), $`𝐤(x)𝐤_{}^\mathrm{T}(x)_x=\mathrm{\Phi }𝚲𝚲_{}\mathrm{\Phi }^\mathrm{T}`$, and $`C_{}(x,x)_x=\mathrm{tr}𝚲_{}`$. This gives, dropping the training inputs subscript from the remaining average,

$$ϵ=\mathrm{tr}𝚲_{}2\mathrm{tr}\mathrm{\Phi }𝚲𝚲_{}\mathrm{\Phi }^\mathrm{T}𝐊^1+\mathrm{tr}\mathrm{\Phi }𝚲^2\mathrm{\Phi }^\mathrm{T}𝐊^1𝐊_{}𝐊^1$$

In this new representation we also have $`𝐊=\sigma ^2𝐈+\mathrm{\Phi }𝚲\mathrm{\Phi }^\mathrm{T}`$ and similarly for $`𝐊_{}`$; for the inverse of $`𝐊`$ we can use the Woodbury formula to write $`𝐊^1=\sigma ^2[𝐈\sigma ^2\mathrm{\Phi }𝒢\mathrm{\Phi }^\mathrm{T}]`$, where $`𝒢=(𝚲^1+\sigma ^2\mathrm{\Phi }^\mathrm{T}\mathrm{\Phi })^1`$. Inserting these results, one finds after some algebra that

$$ϵ=\sigma _{}^2\sigma ^2\left[\mathrm{tr}𝒢\mathrm{tr}𝒢𝚲^1𝒢\right]+\mathrm{tr}𝒢𝚲_{}𝚲^2𝒢$$
(7)
which for the matched case reduces to the known result for the Bayes error

$$\widehat{ϵ}=\mathrm{tr}𝒢$$
(8)
Eqs. (7,8) are still exact. We now need to tackle the remaining averages over training inputs. Two of these are of the form $`\mathrm{tr}𝒢𝐌𝒢`$; if we generalize the definition of $`𝒢`$ to $`𝒢=(𝚲^1+v𝐈+w𝐌+\sigma ^2\mathrm{\Phi }^\mathrm{T}\mathrm{\Phi })^1`$ and define $`g=\mathrm{tr}𝒢`$, then they reduce to $`\mathrm{tr}𝒢𝐌𝒢=g/w`$. (The derivative is taken at $`v=w=0`$; the idea behind introducing $`v`$ will become clear shortly.) So it is sufficient to calculate $`g`$. To do this, consider how $`𝒢`$ changes when a new example is added to the training set. One has

$$𝒢(n+1)𝒢(n)=\left[𝒢^1(n)+\sigma ^2\psi \psi ^\mathrm{T}\right]^1𝒢(n)=\frac{𝒢(n)\psi \psi ^\mathrm{T}𝒢(n)}{\sigma ^2+\psi ^\mathrm{T}𝒢(n)\psi }$$
(9)
in terms of the vector $`\psi `$ with elements $`(\psi )_i=\varphi _i(x_{n+1})`$, using again the Woodbury formula. To obtain the change in $`g`$ we need the average of (9) over both the new training input $`x_{n+1}`$ and all previous ones. This cannot be done exactly, but we can approximate by averaging numerator and denominator separately; taking the trace then gives $`g(n+1)g(n)=\mathrm{tr}𝒢^2(n)/[\sigma ^2+g(n)]`$. Now, using our auxiliary parameter $`v`$, $`\mathrm{tr}𝒢^2=g/v`$; if we also approximate $`n`$ as continuous, we get the simple PDE $`g/n=(g/v)/(\sigma ^2+g)`$ with the initial condition $`g|_{n=0}=\mathrm{tr}(𝚲^1+v𝐈+w𝐌)^1`$. Solving this using the method of characteristics  gives a self-consistent equation for $`g`$,

$$g=\mathrm{tr}\left[𝚲^1+\left(v+\frac{n}{\sigma ^2+g}\right)𝐈+w𝐌\right]^1$$
(10)
The Bayes error (8) is $`\widehat{ϵ}=g|_{v=w=0}`$ and therefore obeys

$$\widehat{ϵ}=\mathrm{tr}𝐆,𝐆^1=𝚲^1+\frac{n}{\sigma ^2+\widehat{ϵ}}𝐈$$
(11)
within our approximation (called ‘LC’ in ). To obtain $`ϵ`$, we differentiate both sides of (10) w.r.t. $`w`$, set $`v=w=0`$ and rearrange to give

$$\mathrm{tr}𝒢𝐌𝒢=g/w=(\mathrm{tr}\mathrm{𝐌𝐆}^2)/[1(\mathrm{tr}𝐆^2)n/(\sigma ^2+\widehat{ϵ})^2]$$

Using this result in (7), with $`𝐌=𝚲^1`$ and $`𝐌=𝚲^1𝚲_{}𝚲^1`$, we find after some further simplifications the final (approximate) result for the learning curve:

$$ϵ=\widehat{ϵ}\frac{\sigma _{}^2\mathrm{tr}𝐆^2+n^1(\sigma ^2+\widehat{ϵ})^2\mathrm{tr}𝚲_{}𝚲^2𝐆^2}{\sigma ^2\mathrm{tr}𝐆^2+n^1(\sigma ^2+\widehat{ϵ})^2\mathrm{tr}𝚲^1𝐆^2}$$
(12)
which transparently shows how in the matched case $`ϵ`$ and $`\widehat{ϵ}`$ become identical.

## 3 Examples

I now apply the result for the learning curve (11,12) to some exemplary learning scenarios. First, consider inputs $`x`$ which are binary vectors111This scenario may seem strange, but simplifies the determination of the eigenfunctions and eigenvalues; for large $`d`$, however, one expects other distributions with continuously varying $`x`$ to give similar results .  with $`d`$ components $`x_a\{1,1\}`$, and assume that the input distribution is uniform. We consider covariance functions for student and teacher which depend on the product $`xx^{}`$ only; this includes the standard choices (e.g. OU and RBF) which depend on the Euclidean distance $`|xx^{}|`$, since $`|xx^{}|^2=2d2xx^{}`$. All these have the same eigenfunctions , so our above assumption is satisfied. The eigenfunctions are indexed by subsets $`\rho `$ of $`\{1,2\mathrm{}d\}`$ and given explicitly by $`\varphi _\rho (x)=_{a\rho }x_a`$. The corresponding eigenvalues depend only on the size $`s=|\rho |`$ of the subsets and are therefore $`(\begin{array}{c}d\\ s\end{array})`$-fold degenerate; letting $`e=(1,1\mathrm{}1)`$ be the ‘all ones’ input vector, they have the values $`\mathrm{\Lambda }_s=C(x,e)\varphi _\rho (x)_x`$ (which can easily be evaluated as an average over two binomially distributed variables, counting the number of $`+1`$’s in $`x`$ overall and among the $`x_a`$ with $`a\rho `$). With the $`\mathrm{\Lambda }_s`$ and $`\mathrm{\Lambda }_s^{}`$ determined, it is then a simple matter to evaluate the predicted learning curve (11,12) numerically. First, though, focus on the limit of large $`d`$, where much more can be said. If we write $`C(x,x^{})=f(xx^{}/d)`$, the eigenvalues become, for $`d\mathrm{}`$, $`\mathrm{\Lambda }_s=d^sf^{(s)}(0)`$ and the contribution to $`C(x,x)=f(1)`$ from the $`s`$-th eigenvalue block is $`\lambda _s(\begin{array}{c}d\\ s\end{array})\mathrm{\Lambda }_sf^{(s)}(0)/s!`$, consistent with $`f(1)=_{s=0}^{\mathrm{}}f^{(s)}(0)/s!`$ The $`\mathrm{\Lambda }_s`$, because of their scaling with $`d`$, become infinitely separated for $`d\mathrm{}`$. For training sets of size $`n=𝒪(d^L)`$, we then see from (11) that eigenvalues with $`s>L`$ contribute as if $`n=0`$, since $`\mathrm{\Lambda }_sn/(\sigma ^2+\widehat{ϵ})`$; they have effectively not yet been learned. On the other hand, eigenvalues with $`s<L`$ are completely suppressed and have been learnt perfectly. We thus have a hierarchical learning scenario, with different scalings of $`n`$ with $`d`$—as defined by $`L`$—corresponding to different ‘learning stages’. Formally, we can analyse the stages separately by letting $`d\mathrm{}`$ at a constant ratio $`\alpha =n/(\begin{array}{c}d\\ L\end{array})`$ of the number of examples to the number of parameters to be learned at stage $`L`$ (note $`(\begin{array}{c}d\\ L\end{array})=𝒪(d^L)`$ for large $`d`$). An independent (replica) calculation shows that our approximation for the learning curve actually becomes exact in this limit. The resulting $`\alpha `$-dependence of $`ϵ`$ can be determined explicitly: Set $`f_L=_{sL}\lambda _s`$ (so that $`f_0=f(1)`$) and similarly for $`f_L^{}`$. Then for large $`\alpha `$,

$$ϵ=f_{L+1}^{}+(f_{L+1}^{}+\sigma _{}^2)\alpha ^1+𝒪(\alpha ^2)$$
(13)
This implies that, during successive learning stages, (teacher) eigenvalues are learnt one by one and their contribution eliminated from the generalization error, giving plateaux in the learning curve at $`ϵ=f_1^{}`$, $`f_2^{}`$, …. These plateaux, as well as the asymptotic decay (13) towards them, are universal , i.e. student-independent. The (non-universal) behaviour for smaller $`\alpha `$ can also be fully characterized: Consider first the simple case of linear perceptron learning (see e.g. ), which corresponds to both student and teacher having simple dot-product covariance functions $`C(x,x^{})=C_{}(x,x^{})=xx^{}/d`$. In this case there is only a single learning stage (only $`\lambda _1=\lambda _1^{}=1`$ are nonzero), and $`ϵ=r(\alpha )`$ decays from $`r(0)=1`$ to $`r(\mathrm{})=0`$, with an over-fitting maximum around $`\alpha =1`$ if $`\sigma ^2`$ is sufficiently small compared to $`\sigma _{}^2`$. In terms of this function $`r(\alpha )`$, the learning curve at stage $`L`$ for general covariance functions is then exactly given by $`ϵ=f_{L+1}^{}+\lambda _L^{}r(\alpha )`$ if in the evaluation of $`r(\alpha )`$ the effective noise levels $`\stackrel{~}{\sigma }^2=(f_{L+1}+\sigma ^2)/\lambda _L`$ and $`\stackrel{~}{\sigma }_{}^2=(f_{L+1}^{}+\sigma _{}^2)/\lambda _L^{}`$ are used. Note how in $`\stackrel{~}{\sigma }_{}^2`$, the contribution $`f_{L+1}^{}`$ from the not-yet-learned eigenvalues acts as effective noise, and is normalized by the amount of ‘signal’ $`\lambda _L^{}=f_L^{}f_{L+1}^{}`$ available at learning stage $`L`$. The analogous definition of $`\stackrel{~}{\sigma }^2`$ implies that, for small $`\sigma ^2`$ and depending on the choice of student covariance function, there can be arbitrarily many learning stages $`L`$ where $`\stackrel{~}{\sigma }^2\stackrel{~}{\sigma }_{}^2`$, and therefore arbitrarily many over-fitting maxima in the resulting learning curves. Fig. 1(left) demonstrates that this conclusion holds not just for $`d\mathrm{}`$; even for the cases shown, with $`d=10`$, up to three over-fitting maxima are apparent. Our theory provides a very good description of the numerically simulated learning curves even though, at such small $`d`$, the predictions are still significantly different from those for $`d\mathrm{}`$ (see Fig. 1(right)) and therefore not guaranteed to be exact.

In the second example scenario, I consider continuous-valued input vectors, uniformly distributed over the unit interval $`[0,1]`$; generalization to $`d`$ dimensions ($`x[0,1]^d`$) is straightforward. For covariance functions which are stationary, i.e. dependent on $`x`$ and $`x^{}`$ only through $`xx^{}`$, and assuming periodic boundary conditions (see  for details), one then again has covariance function-independent eigenfunctions. They are indexed by integers222Since $`\mathrm{\Lambda }_q=\mathrm{\Lambda }_q`$, one can assume $`q0`$ if all $`\mathrm{\Lambda }_q`$ for $`q>0`$ are taken as doubly degenerate.  $`q`$, with $`\varphi _q(x)=e^{2\pi iqx}`$; the corresponding eigenvalues are $`\mathrm{\Lambda }_q=𝑑xC(0,x)e^{2\pi iqx}`$. For the (‘periodified’) RBF covariance function $`C(x,x^{})=\mathrm{exp}[(xx^{})^2/(2l^2)]`$, for example, one has $`\mathrm{\Lambda }_q\mathrm{exp}(\stackrel{~}{q}^2/2)`$, where $`\stackrel{~}{q}=2\pi lq`$. The OU case $`C(x,x^{})=\mathrm{exp}(|xx^{}|/l)`$, on the other hand, gives $`\mathrm{\Lambda }_q(1+\stackrel{~}{q}^2)^1`$, thus $`\mathrm{\Lambda }_qq^2`$ for large $`q`$. I also consider below covariance functions which interpolate in smoothness between the OU and RBF limits: E.g. the MB2 (modified Bessel) covariance $`C(x,x^{})=e^a(1+a)`$, with $`a=|xx^{}|/l`$, yields functions which are once differentiable ; its eigenvalues $`\mathrm{\Lambda }_q(1+\stackrel{~}{q}^2)^2`$ show a faster asymptotic power law decay, $`\mathrm{\Lambda }_qq^4`$. To subsume all these cases I assume in the following analysis of the general shape of the learning curves that $`\mathrm{\Lambda }_qq^r`$ (and similarly $`\mathrm{\Lambda }_q^{}q^r_{}`$). Here $`r=2`$ for OU, $`r=4`$ for MB2, and (due to its faster-than-power law decay of eigenvalues) effectively $`r=\mathrm{}`$ for RBF.

From (11,12), it is clear that the $`n`$-dependence of the Bayes error $`\widehat{ϵ}`$ has a strong effect on the true generalization error $`ϵ`$. From previous work , we know that $`\widehat{ϵ}(n)`$ has two regimes: For small $`n`$, where $`\widehat{ϵ}\sigma ^2`$, $`\widehat{ϵ}`$ is dominated by regions in input space which are too far from the training examples to have significant correlation with them, and one finds $`\widehat{ϵ}n^{(r1)}`$. For much larger $`n`$, learning is essentially against noise, and one has a slower decay $`\widehat{ϵ}n^{(r1)/r}`$. These power laws can be derived from (11) by approximating factors such as $`[\mathrm{\Lambda }_q^1+n/(\sigma ^2+\widehat{ϵ})]^1`$ as equal to either $`\mathrm{\Lambda }_q`$ or to $`0`$, depending on whether $`n/(\sigma ^2+\widehat{ϵ})<`$ or $`>\mathrm{\Lambda }_q^1`$. With the same technique, one can estimate the behaviour of $`ϵ`$ from (12). In the small $`n`$-regime, one finds $`ϵc_1\sigma _{}^2+c_2n^{(r_{}1)}`$, with prefactors $`c_1`$, $`c_2`$ depending on the student. Note that the contribution proportional to $`\sigma _{}^2`$ is automatically negligible in the matched case (since then $`ϵ=\widehat{ϵ}\sigma ^2=\sigma _{}^2`$ for small $`n`$); if there is a model mismatch, however, and if the small-$`n`$ regime extends far enough, it will become significant. This is the case for small $`\sigma ^2`$; indeed, for $`\sigma ^20`$, the ‘small $`n`$’ criterion $`\widehat{ϵ}\sigma ^2`$ is satisfied for any $`n`$. Our theory thus predicts the appearance of plateaux in the learning curves, becoming more pronounced as $`\sigma ^2`$ decreases; Fig. 2(left) confirms this333If $`\sigma ^2=0`$ exactly, the plateau will extend to $`n\mathrm{}`$. With hindsight, this is clear: a GP with an infinite number of nonzero eigenvalues has no limit on the number of its ‘degrees of freedom’ and can fit perfectly any amount of noisy training data, without ever learning the true teacher function.. Numerical evaluation also shows that for small $`\sigma ^2`$, over-fitting maxima may occur before the plateau is reached, consistent with simulations; see inset in Fig. 2(right). In the large $`n`$-regime ($`\widehat{ϵ}\sigma ^2`$), our theory predicts that the generalization error decays as a power law. If the student assumes a rougher function than the teacher provides ($`r<r_{}`$), the asymptotic power law exponent $`ϵn^{(r1)/r}`$ is determined by the student alone. In the converse case, the asymptotic decay is $`ϵn^{(r_{}1)/r}`$ and can be very slow, actually becoming logarithmic for an RBF student ($`r\mathrm{}`$). For $`r=r_{}`$, the fastest decay for given $`r_{}`$ is obtained, as expected from the properties of the Bayes error. The simulation data in Fig. 2 are compatible with these predictions (though the simulations cover too small a range of $`n`$ to allow exponents to be determined precisely).

In summary, the above approximate theory makes a number of non-trivial predictions for GP learning with mismatched models, all borne out by simulations: for large input space dimensions, the occurrence of multiple over-fitting maxima; in lower dimensions, the generic presence of plateaux in the learning curve if the student assumes too small a noise level $`\sigma ^2`$, and strong effects of model mismatch on the asymptotic learning curve decay. The behaviour is much richer than for the matched case, and could guide the choice of (student) priors in real-world applications of GP regression; RBF students, for example, run the risk of very slow (logarithmic) decay of the learning curve if the target (teacher) is less smooth than assumed. An important issue for future work is to analyse to which extent hyperparameter tuning (e.g. via evidence maximization) can make GP learning robust against some forms of model mismatch, e.g. a misspecified functional form of the covariance function.


text:
# Public Release of 2dF data from the Fornax Cluster Spectroscopic Survey

## 1 Spectroscopy of Everything

Almost 20 years ago Morton, Krug & Tritton (1985) made the first complete spectroscopic survey of a given region of sky. They surveyed stellar objects in a 0.3 deg2 region of sky to a limit of $`B=20`$, measuring 600 objective prism spectra with the UKST and 100 slit spectra with the AAT. All the objects measured were normal stars except for a small number of QSOs and white dwarfs, the latter being used to put new constraints on the white dwarf luminosity function.

Thanks to the arrival of the 2dF spectrograph on the AAT, we have recently completed the first stage of a complete spectroscopic survey more than one order of magnitude larger than the Morton et al. study, measuring 7000 spectra in a $`2\pi `$ deg2 area as part of our study of the Fornax Cluster. In this article we describe the public release of 3600 spectra from our first field. We hope that this public release will encourage colleagues making surveys for rare objects to choose these fields, as much of the follow-up spectroscopy that might be required is available from our data.

Our 2dF Fornax Cluster Spectroscopic Survey (FCSS; see Drinkwater et al. 2000a) was designed to make the most complete census possible of low-luminosity galaxies in the Fornax Cluster. As well as the conventional low-surface brightness dwarfs, we already had evidence (Drinkwater & Gregg 1998) that very compact, high surface brightness dwarf cluster galaxies had been missed in previous work, so we took the unusual step of observing all objects in each 2dF field, both resolved (“galaxies”) and unresolved (“stars”). In this way we avoided any morphological bias as to what a cluster galaxy should look like. Observing all the “stars” in each 2dF field leads to a large increase in the number of targets observed: within our magnitude limits ($`16.5<b_J<19.8`$) over 50% of the objects in each 2dF field are stars. Thanks to the flexibility of 2dF, however, this only leads to a small increase in total observing time as we schedule the stars at times when we could not usefully observe galaxies, such as twilight or at high airmass or through cloud.

## 2 Status of the Survey

As of 2001 January we have completed observations of two of our four planned 2dF fields centred on the Fornax Cluster. We have analysed the data from our first field and have released it for public access as described below. The properties of the sample are summarised in Fig. 1, a plot of surface brightness against magnitude for all the objects we successfully observed. The main all-object sample was selected in the magnitude range $`16.5<b_J<19.8`$. We extended our observations of the unresolved objects (“stars”) to a slightly fainter magnitude limit of $`b_J20.2.`$ There is some incompleteness for galaxies with surface brightness lower than about 23.5 $`b_J`$ mag arcsec-2 as we were unable to measure their redshifts (see Drinkwater et al. 2000a).

The diversity of objects in our sample is illustrated by Fig. 2, a cone diagram of the sample in which the distance axis is scaled (as $`z^{1/4}`$) so as to display both Galactic stars and the most distant QSOs. Note that although we avoid the use of conventional morphological classifications in our analysis (we use redshift, spectral signatures and luminosity instead), we do show the morphological classifications in both Fig. 1 and Fig. 2. The presence of unresolved objects (“stars”) in both the Fornax Cluster and among the background field galaxies demonstrates the importance of our all-object strategy. These objects would have been missed by conventional galaxy surveys only looking at objects that are resolved on UKST sky survey plates.

The scientific highlight of our results to date has been the discovery of a new population of Fornax Cluster dwarf galaxies so compact they were previously mistaken for Galactic stars (Drinkwater et al. 2000b). They can be seen among the more normal, resolved Fornax Cluster dwarf galaxies in Fig. 2. These ”ultra-compact dwarf” (UCD) galaxies are unlike any known type of stellar system. They are smaller and more concentrated than any known dwarf galaxy, but are 2-3 magnitudes more luminous than the largest Galactic globular clusters. Numerical simulations have shown that the UCDs could have been formed by tidal stripping of nucleated dwarf galaxies in close orbits around the central galaxy of the cluster, NGC 1399 (Bekki, Couch & Drinkwater 2001). We are making detailed follow-up observations of the UCDs with Hubble Space Telescope imaging to measure their sizes along with high resolution spectroscopy on the VLT to measure their velocity dispersions and hence masses.

Other results from the early stages of our survey include the first large QSO sample selected purely from spectroscopic identification without any pre-selection bias (Fig. 3; see Meyer et al. 2001) and the identification of a population of $`L_{}`$ compact field galaxies also unresolved from the ground (Drinkwater et al. 1999). We have also shown that some low surface brightness galaxies previously assumed to be members of the Fornax Cluster are in fact background objects (Jones et al., in preparation).

## 3 Public Web Access

The 3600 2dF spectra from our first field are now available for public access at our web site
http://astro.ph.unimelb.edu.au/data/ in Melbourne. The main function of the web site is to allow searches of our database by position on the sky and to return our 2dF spectra and redshift measurements of the selected objects. The complete database of measurements is available from the website as a plain text file as well as additional help material and references. The database consists of all the objects for which we have measured reliable redshifts (about 95% of the sample at present). Our 2dF spectra were all obtained with the same observing configuration as for the 2dF Galaxy Redshift Survey: a wavelength coverage of 3600–8010Å and a resolution of 9Å (dispersion of 4.3Å per pixel). The spectra have typical signal-to-noise ratios of 10 or more per pixel. The redshifts were measured by cross-correlation with star and galaxy templates; this gives velocities accurate to about 65$`\mathrm{km}\mathrm{s}^1`$ for galaxies (see Drinkwater et al. 2000a for details).

The database is searched by entering a co-ordinate and search radius, with a table listing the parameters of any objects observed within that radius being returned. The parameters listed include magnitudes and classifications from the APM Catalogues (Irwin, Maddox & McMahon, 1994) as well as our own redshift measurements. Following our philosophy of avoiding morphological bias in classifications, we do not specifically classify objects in the database (although the APM classifications are listed for reference). Instead we recommend the use of our spectroscopic redshifts (in conjunction with the magnitudes if necessary) to make classifications on a more physical basis, e.g. classifying objects with redshifts less than 700$`\mathrm{km}\mathrm{s}^1`$ as Galactic stars.

To obtain the 2dF spectrum of any object returned in the output table, click on the “plot” button for that object. When you click on the plot button, the spectrum will appear in a new window. The spectrum can then be replotted with different axis limits if necessary and downloaded in text, GIF, PostScript or FITS format. An example of the PostScript output is shown in Fig. 4.

## References

Bekki, K., Couch, W.J., Drinkwater M.J., 2001, ApJ, in press

Drinkwater, M.J., Gregg, M.D., 1998, MNRAS, 296, L15

Drinkwater, M.J., et al., 1999, ApJ, 511, L97

Drinkwater, M.J., et al., 2000a, A&A,

Drinkwater M.J., Jones, J.B., Gregg, M.D., Phillipps, S., 2000b, PASA, 17, 227

Irwin M.J., Maddox S., McMahon R., 1994, Spectrum, 2, 14

Meyer M.J., Drinkwater M.J., Phillipps S., Couch W.J., 2001, MNRAS, in press

Morton D.C., Krug P.A., Tritton K.P., 1985, MNRAS, 212, 325


text:
# Disorder induced hexagonal-orthorhombic transition in Y{_{1-𝑥}}{^limit-from3+}Gd{_𝑥}{^limit-from3+}MnO3

## I Introduction

In the search for new composition-properties relations ABO3 compounds have attracted a lot of attention. The perovskite materials, ABO3, have been researched extensively because this structure forms the basis for interesting physical properties such as high $`T_c`$ superconductivity and colossal magnetoresistance. Non-perovskite AMnO3, with A = Y, Ho,…,Lu, attracted renewed interest, due to their ferroelectric properties. These hexagonal AMnO3 have a basically different structure than most ABO3 compounds, that are distorted perovskites. These properties arise due to the strong correlation of the $`3d`$ electrons with the O $`2p`$ orbitals.

In this paper, we report the transition from the orthorhombic perovskite (o) to the hexagonal (h) phase by changing the ionic radius of the A ion and by high pressure synthesis. The effect of the ionic radius on the transition is studied by partially replacing Y by Gd ions in YMnO3. The resulting phase diagram leads us to discuss the effect of disorder, in terms of the ionic radius variance, on the stability of the hexagonal and orthorhombic phases.

The basic building block of the perovskite is an oxygen octahedron with a transition metal, B, in its centre. The A ions, usually lanthanides or alkaline earth metal ions, occupy the holes between the octahedra, that form a 3D corner shared network. In this picture B is sixfold and A is 12-fold coordinated. Most perovskites have a distorted structure, derived from this building block. The distortions have various origins, including a ferroelectric transition for B a $`d^0`$ transition metal ion like Ti4+ . The most common distortion originates from the relative small radius of the A ions compared with the holes between the octahedra. This results in a cooperative rotation of the octahedra known as the GdFeO3 distortion. While the structure is interesting in its own right, it has also large effects on the physical properties. It is well documented that the physical properties depend strongly on the magnitude of the structural distortions. An overview for the manganites is given in Ref.’s .

The magnitude of the GdFeO3 distortion depends strongly on the tolerance factor, $`t`$:

$$t=\frac{r_{A^{3+}}+r_{O^2}}{\sqrt{2}(r_{B^{3+}}+r_{O^2})},$$
(1)
where $`r_X`$ is the radius of the X ion. The tolerance factor gives the relation between the radii of ions A, B and O in a ideal cubic perovskite. For $`t=1`$ the size of the lanthanide is exactly right to compose the cubic perovskite system. For Mn3+, $`r_{Mn^{3+}}=0.645`$ Å and $`r_{O^2}=1.42`$ Å this yields a ionic radius $`r_{A^{3+}}=1.50`$ Å, where the largest lanthanide, La, has a radius of 1.22 Å. The corresponding tolerance factor $`t=0.90`$ indicates a large distortion for LaMnO3. With increasing atomic number, the lanthanide radius decreases and thereby the distortion increases. For the manganites, the tolerance factor is conventionally regarded as the factor controlling the boundary between the hexagonal and orthorhombic structures. The orthorhombic perovskite phase is stable for $`t>0.855`$, corresponding to $`r_Ar_{Dy}`$. For $`t<0.855`$, $`r_Ar_{Ho}`$, the hexagonal phase prevails . Yttrium, although not in the lanthanide series, behaves chemically identical and its radius falls between dysprosium and holmium. An overview on the ionic radii and tolerance factors of some relevant compounds is given in Table I.

The high temperature phase of the hexagonal AMnO3 consists of 8-fold coordinated A ions in bicapped antiprisms. Trigonal bipyramidal holes are formed between two layers of face-sharing antiprisms by the edges of the capping oxygens with the antiprisms. The capping oxygens of two adjacent layers are located on the same $`ab`$ plane. Half of the bipyramidal holes are occupied by Mn. The apical oxygens of the MnO5 bipyramid are also the oxygens that make up the antiprism. The two polyhedra are sketched in Fig. 1, where the shared edge is shown. The Mn-Oap distance is thus equal to the distance between the antiprism oxygen layer and the capping oxygen layer. The steric hindrance of the Mn restricts this layer separation and therefore increases the A-Oap bond length. Thus, the eightfold co-ordination is not uniform. The two apical oxygens have slightly larger bond lengths. Furthermore, the structure is unstable against a ferroelectric distortion at lower temperatures. The apical oxygens move in such a way that one bond becomes ’normally’ short, while the other becomes about 1 Å larger. The asymmetric A environment is the main reason for the ferroelectric behaviour. As we have two lanthanide positions in $`P6_3cm`$, we have two non-equivalent, although similar, dipole moments. Four out of the six moments per unit cell point upwards, the other two downwards.

Although the hexagonal phase of YMnO3 at ambient conditions is the thermodynamically stable phase, there are several ways to obtain orthorhombic YMnO3. Using thin film growth, an appropriate substrate will force the coherent growth of the orthorhombic phase . Synthesis routes via organic precursors and low reaction temperatures yield the orthorhombic phase . And last, high pressure synthesis favours the orthorhombic phase, because it has a higher density .

## II Experimental

Polycrystalline ceramic samples of AMnO3, where A is a mixture of Y and Gd, have been synthesised using regular solid state synthesis at ambient pressure. Starting materials were Y2O3, Gd2O3 and MnO2. Stoichiometric amounts corresponding to formulae which range from pure YMnO3 to pure GdMnO3 were weighted and wet-mixed using acetone as liquid medium. The pressed pellets were sintered for 24 hours at 1250 and for 24 hours at 1400.

High pressure experiments were carried out both on the mixture of oxides and on as-prepared samples. X-ray powder diffraction patterns were identical for both methods. The high pressure high temperature piston cylinder apparatus is a Depth of the Earth Quickpress 3.0, with an experimental range up to 25 kbar and 2100C . The lower pressure limit is $`1`$ a 2 kbar.

The sample environment is a complex set-up, including a graphite furnace and the pressure medium. Care has to be taken to prevent contamination from the graphite resistance furnace or any of the other materials, *e.g.* Al2O3 or NaCl, in the sample assembly. Therefore, a small amount of powder $`0.5`$ g was encapsulated in a Pt capsule. The capsule consisted of a tube (diameter 4 mm and height 6 mm) and two pre-shaped lids, which were welded together using a small welding apparatus. Recovery of the pellet from the sample assembly is improved by the Pt container.

X-ray diffraction patterns were recorded using a Bruker-AXS D8 powder diffractometer, with primary and secondary monochromator, using Cu $`K_\alpha `$ radiation. Patterns were analysed for phase determination using the evaluation software EVA , the Powder Diffraction File  and the Inorganic Crystal Structure Database . Patterns of non-contaminated samples, containing only the hexagonal and orthorhombic phases, were used in the Rietveld refinement using TOPAS R . Rietveld refinements included lattice parameters, zero point correction and the ratio between the two phases. Atomic positions were assumed to be constant using the positions determined by single crystal x-ray diffraction on YMnO3 . The ratio Y:Gd was fixed at the nominal composition.

## III Results and discussion

In Fig. 2, we present a phase diagram of Y1-xGdxMnO3 as a function of $`r_A`$ and its variance $`\sigma ^2`$ of samples given in Table II. The variance $`\sigma ^2`$ is given by,

$$\sigma ^2=\underset{1}{\overset{n}{}}x_i(r_ir_A)^2$$
(2)
The phase diagram can be divided in three regions. Low $`r_A`$ compounds, $`r_A<1.078`$, are hexagonal. Large $`r_A`$ and small $`\sigma ^2`$ compounds are orthorhombic. The intermediate region shows both phases. The data for HoMnO3 and DyMnO3 have been taken from the literature .

By changing the value of $`x`$ we substitute Y by Gd, whereby $`r_A`$ increases linearly. The tolerance factor of orthorhombic DyMnO3 is equal to that of Y1-xGdxMnO3, with $`x=0.23`$. Thus for $`x0.23`$ we expect to see the orthorhombic phase. For smaller $`x`$ and $`t`$ the hexagonal phase is expected. For $`x=0`$ and $`x=0.06`$ we indeed found the hexagonal structure. However, for $`0.19x0.38`$ we do not observe a sharp transition to the orthorhombic structure, but a mixture of the hexagonal and orthorhombic phases. Only for $`x=0.5`$ an (almost) pure orthorhombic compound is found. The anomalous behaviour of the mixed Y,Gd samples is best illustrated by focussing on the sample with $`x=0.25`$, which has an almost identical tolerance factor as DyMnO3. Where the latter sample is orthorhombic, the former segregates in both phases. The only difference between the two compounds is that one is undoped and the other has a mixed lanthanide composition.

The relative amounts of orthorhombic and hexagonal fractions are determined by Rietveld refinement of the powder diffraction data, see Fig. 4. We can rule out element segregation of Y and Gd, since the continuous increase in the lattice parameter $`a`$ of the hexagonal phase with increasing Gd concentration indicates perfect mixing of the A ions. A linear increase in $`a`$ with $`r_A`$ is also observed for the single A cation h-AMnO3 series . Fig. 4 shows that the orthorhombic phase fraction increases linearly with the Gd fraction. This allows us to apply the lever rule on the h-o transition

$$\frac{xx_{hexa}}{xx_{ortho}}=\frac{y_{ortho}}{y_{hexa}}$$
(3)
where $`x_{hexa}`$ and $`x_{ortho}`$ are the boundary values for the respective phases and $`y_{ortho}/y_{hexa}`$ the ratio of the two fractions. From the observed ratios $`y_{ortho}/y_{hexa}`$ as a function of $`x`$, the boundary values are derived. The boundary values derived from all five mixed phase samples are plotted in Fig. 2 as inverted triangles.

We construct a preliminary phase diagram by assuming that the h-o transition at $`\sigma ^2=0`$ occurs halfway between hexagonal YMnO3 and orthorhombic DyMnO3. The phase boundaries are drawn in Fig. 2 as straight lines through the two calculated boundary values and the assumed $`\sigma ^2=0`$ midpoint. This phase diagram can be described as follows. The phase line associated with the upper $`t`$ limit of the hexagonal phase does not depend on the variance. Slightly increasing $`\sigma ^2`$ and $`t`$ results in the appearance of a two phase region, consisting of both the hexagonal and the orthorhombic phase. With increasing tolerance factor, the fraction of the orthorhombic phase increases until the lower boundary limit for the orthorhombic phase is crossed. This limit strongly depends on $`\sigma ^2`$. The lower limit for the orthorhombic phase increases from $`r_A=1.078`$ Å at $`\sigma ^2=0`$ to $`r_A=1.093`$ Å at $`\sigma ^2=\sigma _{max}^2`$. Note that $`r_A=1.093`$ Å corresponds with the ionic radius of terbium, the second smallest lanthanide to form the perovskite structure.

We have shown the dependence of the h-o transition on the average radius $`r_A`$ and $`\sigma ^2`$. In the next section, the effect of high pressure experiments on the h-o transition will be discussed.

We applied high pressure and high temperature to convert or synthesise some of the conventionally hexagonal samples in the orthorhombic state. Pressure generally stabilises the most dense phase, in this case the orthorhombic structure. First, we consider samples with $`\sigma ^2=0`$. YMnO3 is still hexagonal at 5 kbar, but the orthorhombic phase is found using a pressure of 15 kbar. The data is shown in Fig. 5. The necessary pressure for the h-o transition, less than 15 kbar, is much less than reported previously in the literature . The h-o transition phase line has been sketched in Fig. 5 by using midpoints. Extrapolating the pressure dependence of the h-o transition yields a critical pressure of $`27`$ kbar for HoMnO3. The error bar on this value, $`10`$ kbar, is large because of the sparse data points.

We also carried out high pressure experiments on the Y1-xGdxMnO3 compounds. The necessary pressure to induce the orthorhombic phase from the phase mixture decreased with increasing $`r_A`$.

We note the following observations of the two-phase region:

1. The phase mixture is only observed for ambient presurre synthesis for samples with $`\sigma ^20`$.
2. We do not observe the mixed phase for any of the experiments at high pressure ($`5`$ kbar). We cannot exclude the presence of the mixed phase at low pressures, but our experimental set-up is not well-suited for those pressures.
3. The lattice parameters of both phases in the two phase region indicate no segregation into Y-rich and Gd-rich phases, see Fig. 3, which is unconventional.
4. Literature reports that synthesis via organic precursors can result in a mixture of hexagonal and orthorhombic phases for AMnO3, $`\sigma ^2=0`$, compounds .

These observations lead to the following conclusions. For $`\sigma ^2=0`$ compounds either the hexagonal or the orthorhombic structure is stable. Whereas low temperature synthesis may yield mixed phase samples, a high temperature anneal will convert the unstable phase. The hexagonal or orthorhombic phase will be stable depending on the tolerance factor. However for $`\sigma ^20`$ mixed phase samples can be obtained for a broad range of tolerance factors. Even a high temperature anneal retains the phase segregated state. Surprisingly, the phase segregation is not accompanied by two limiting compositions, *e.g* Y1-xGdxMnO3 with $`x=0.1`$ and $`x=0.5`$. The continuous increase of the lattice parameters of both the hexagonal and the orthorhombic phase throughout the two phase region indicates that the composition of the hexagonal and the orthorhombic state are the same in the mixed state. We have no explanation for this unconventional form of phase segregation.

## IV Conclusions

We have constructed phase diagrams for the h-o phases of AMnO3, including the effects of average ionic radius, hydrostatic pressure and variance. For compounds with $`\sigma ^20`$ a phase separation in the orthorhombic and the hexagonal phase is found. The mixed region exists only at low pressures. We have shown that at ambient pressure, this region expands towards higher values of the average radius with increasing variance. The upper limit for the hexagonal phase is not affected by an increase in the variance. We speculate that disorder, introduced by a large variance or soft chemical synthesis routes, allows the occurrence of the phase separation in the absence of other driving forces. Suppressing the disorder by applying external pressure or annealing at high temperatures prevents the existence of the phase separation. Pressure favours the denser, orthorhombic phase, whereas thermal annealing promotes the hexagonal phase.

This work is supported by the Netherlands Foundation for the Fundamental Research on Matter. (FOM).