text:
# NACO Polarimetric Differential Imaging of TW Hya:

## 1 Introduction

The simulations of viscous disks with an embedded Jupiter-mass planet (see, e.g. D’Angelo et al. dan03 (2003)) predict dramatic changes in the disk structure: gaps are being opened and shock waves follow the path of the orbiting body. Detecting these signatures is an important step in confirming that our models for planet formation and disk evolution are right and would provide a unique possibility to study a newly formed exoplanet and its interaction with the circumstellar material.

However, the investigations of protoplanetary disks require high-contrast, high-spatial-resolution imaging very close to the bright central star. Although in recent years numerous groups studied disks around nearby stars at their planet-forming age with the HST, instrumental limitations do not allow to probe the disk structure at scales similar to our inner Solar System. The technique of polarimetric differential imaging (PDI) as described by Kuhn et al. (kuh01 (2001)) is capable of tracing scattered (i.e. polarized) light from a circumstellar disk very close to the central star. This promising method enhances the contrast between the disk and the star and, thereby, might be important to detect the detailed inner structure of protoplanetary disks.

The basic idea of the PDI is to take the difference of two orthogonally polarized, simultaeneously acquired images of the same object in order to remove all non-polarized light components. The non-polarized light consists mainly of the central star’s light and the speckle noise, which is – close to the star – the dominant noise source of ground-based observations at optical and near-infrared wavelengths (Racine et al. rac99 (1999)). After this subtraction only the polarized light, such as scattered light from the disk, is left over.

To explore the potential of this technique for planet formation studies as well as to demonstrate the NACO/VLT platform capabilities we targeted the 8 Myr old (Webb et al. web99 (1999)) and very close ($`56\pm 7`$ pc, Wichmann et al. wic98 (1998)) classical T Tauri star TW Hya (RA: $`11^\mathrm{h}01^\mathrm{m}51.9^\mathrm{s}`$ Dec:$`34\mathrm{°}42^{}17^{\prime \prime }`$ J2000). At the visual wavelength regime TW Hya shows a polarization variability between 0-3% (Mekkaden mek98 (1998)). Being an excellent target for planet formation studies, during the past few years several attempts have been made to image its disk as close as possible to the star and to identify possible disk structures (Krist et al. kri00 (2000), Weinberger et al. wei02 (2002)). In addition, the presence of a giant planet was predicted (Calvet et al. cal02 (2002)) from modeling the spectral energy distribution (SED) (see also Steinacker & Henning ste03 (2003)). The first application of the PDI to TW Hya by Potter (2003a , 2003b ) made use of the 36 element curvature sensing Adaptive Optics (AO) system Hokupa’a (Graves et al. gra00 (2000)), of the Gemini North telescope. This H-band data set probed the circumstellar material closer than 0.5″to the star.

In the following we describe the first results obtained with the PDI technique and the 188 element Shack-Hartmann AO system NACO, which is attached to the Nasmyth B port of the VLT UT4.

## 2 Observations and Data Reduction

The observations were carried out in April 2002 during the commissioning of the NACO/VLT system (Lenzen et al. len98 (1998), Rousset et al. rou98 (2003), Hartung et al. har00 (2000)) in the Ks band. In order to simultaneously measure two orthogonal components of the polarized light, a Wollaston prism was introduced in the light path. The Wollaston prism splits the light into an ordinary and an extraordinary beam (o\- and e-beam), separated by 3.5″ in the Ks-band. To eliminate the instrumental polarization and to increase the signal-to-noise ratio (SNR) of the data, a redundant data set with Wollaston angles of 0$`\mathrm{°}`$, 45$`\mathrm{°}`$, 90$`\mathrm{°}`$ and 135$`\mathrm{°}`$ was acquired. Each of these images has a field of view (FOV) of 3$`\mathrm{}\times 29.4\mathrm{}`$ sampled with a pixel size of 0.027$`\mathrm{}`$/pixel. At every Wollaston position a 3-point dithering was applied with steps of about 9$`\mathrm{}`$ along the field of view to allow sky subtraction and to reduce the influence of bad pixels. No additional off-source sky frames were taken.

To ensure the high dynamic range needed for disk analysis, we repeated the complete set of these observations with defined integration times of 0.4 s and 30 s. The total time spent on the source was 24 s and 1800 s in the short and long exposure series, respectively. As the TW Hya saturated the detector, the AO-performance can not be accurately measured. However, based on simulations with the NAOS Prepartion Software v1.62 we estimate a typical Strehl-ratio of 50% for these observations.

The data reduction was carried out using self–developed IDL routines. First, every frame was manually inspected and those showing reflections or electronic ghosts were excluded from further reduction. Hot pixels were efficiently removed by a 3.5-sigma filtering process. The raw images were sky subtracted and flat field corrected. The sky frame was calculated individually for each group of images taken with the same polarization angle. Following these basic corrections, the ordinary and extraordinary peaks of the individual exposures were centered and extracted. To remove all of the non-polarized intensity (mainly from the central star) we subtracted images of ordinary polarization from those of extraordinary polarization.

The fine alignment of the individual images before subtraction and co-addition was carried out by a two-level gaussian fitting procedure (see also Section 4).

The subtraction process gives the orthogonal $`Q_i`$ and $`U_i`$ Stokes components, $`i`$ standing for the number of dithering positions. For the 4 angles of the Wollaston prism we thus derived the values for $`Q_i^0`$, $`U_i^{45}`$, $`Q_i^{90}`$ and $`U_i^{135}`$ (see, e.g. Huard hua97 (1997)). Using the redundancy of our data set we derived the $`Q`$ and $`U`$ mean polarization vectors:

$$Q=\underset{i=1}{\overset{3}{}}\frac{(Q_i^0Q_i^{90})}{2}\mathrm{and}U=\underset{i=1}{\overset{3}{}}\frac{(U_i^{45}U_i^{135})}{2}$$

and the polarized intensity:

$$PI=\sqrt{\frac{Q^2+U^2}{2}}.$$

This reduction procedure was repeated independently on both the long and the short exposure series, resulting in the $`PI_{\mathrm{Short}}`$ and $`PI_{\mathrm{Long}}`$ polarized intensity maps.

Due to the brightness of TW Hya, the inner regions of the raw frames exceeded the regime of linear detector response. These data points ($`r<0.06\mathrm{}`$ in the short exposure series and $`r<0.45\mathrm{}`$ in the long exposure series) were excluded from the data reduction and further analysis. The field of view (FOV) was limited to $`3\mathrm{}`$ by the field mask. During the dithering process slight movements perpendicular to the dithering directions occured. Since the resulting FOV is the intersection of the FOVs of the individual frames, the final field of view was reduced to $`0.54\mathrm{}`$ and $`1.43\mathrm{}`$ for the short and long exposure series, respectively.

To characterize the surface brightness distribution of the disk we measured the radial profile of the polarized intensity, derived from 3-pixel-wide apertures, roughly equivalent to the full width at half maximum (FWHM) of the point spread function (PSF).

## 3 Signal Analysis

For detecting the polarization signal as close as possible to the star a very careful data inspection is essential – residuals from the subtractions as well as sub-pixel misalignments between the individual frames are the dominant noise sources. While the scalar PI image contains no information on the direction of the polarization, the Q- and U-components are projections of the polarization vector and thus show a sinusoidal modulation with the angle in respect to the polarization direction of the Wollaston-prism (i.e. a maximum at the angles parallel to the polarization direction of the Wollaston prism and a minimum perpendicular to it). In the case of axisymmetrically distributed dust this modulation leads to the butterfly pattern as shown in the right panel of Fig. 1 and described in Kuhn et al. (kuh01 (2001)). By using the information about the polarization orientation and by comparing this pattern to what is expected from an axi-symmetric light scatterer, we can better discriminate between noise and signal.

The right panel of Fig. 1 shows the Q-component of the polarization signal radially integrated (over 5-degree bins) in an annulus between 0.75″ and 1.0″ from the long exposure series. A very similar, but $`\pi `$/2-shifted modulation is evident in the U-images. The clearly sinusoidal modulation arises from the strong polarization signal. When probing the extent of the polarized circumstellar material, we looked for the presence of a similar pattern (phase- and frequency-correct sinusoidal modulation) as a function of radial distance from TW Hya.

We found an extended polarization pattern with a high signal-to-noise ratio between 0.5 and 1.4″ distance from the star in both the Q and U long exposure images. In the short exposure images we identified light-scattering material between 0.2″ and 0.4″ from TW Hya in both the Q- and U-components. Inspection of the higher signal-to-noise Q-component (short exposure) showed this pattern also between 0.1″and 0.2″ (see the left panel of Fig. 2).

In order to estimate the contrast ratio between the polarized and the unpolarized fluxes we compared the fluxes of the direct images (after sky subtraction and bad pixel removal) to the amplitude of the fitted cosine. The typical contrast ratios are $`F_{P/U}`$1% and $`F_{P/U}`$3.5–4% for radii of $`r=0.1\mathrm{}`$ – $`0.3\mathrm{}`$ and $`r1\mathrm{}`$, respectively.

## 4 Reliability Tests and Confirmation

To confirm our results, we analyzed the Q- and U-components of the PDI data set of the PSF-comparison star GSC 07208-00319, which was observed and reduced in an identical way to the TW Hya data set. In this case neither polarized emission in the PI image nor any sign of sinusoidal modulation in the Q- and U-band images could be found (see right panel of Fig. 3).

We also conducted extensive tests to check the reliablity of our data reduction and analysis toolkit. During these test runs the effects of the following errors were investigated: sub-pixel shifts between the individual frames, different centering, interpolation, and flat field normalization methods. For each test run the complete data reduction and analysis has been executed both for TW Hya and the PSF comparison star.

As the possible sub-pixel misalignments might produce a bipolar, positive-negative pattern mimicking the butterfly pattern of the extended polarization signal, special care was taken to investigate their effects on our results. In order to simulate centering problems, we have introduced normally distributed random shifts with FWHMs of 0.2, 0.5, 1.0 and 1.5 pixels after the centering step. The sinusoidal modulation of the TW Hya short exposure frames gets strongly distorted and asymmetric even for the sub-pixel misalignments (0.2 and 0.5) and the pattern is completely disrupted for misalignments as large as 1 or 1.5 pixels. The short exposure images of the PSF comparison star start to show a sinusoidal modulation, but its frequency is one half of that of the polarization signal. These images show, in fact, offset negative and positive peaks which can easily be discriminated from the butterfly pattern shown in Fig. 1. The long exposure images are less sensitive to sub-pixel misalignment but are noticeably distorted in the case of large (1 or 1.5 pixel) centering errors.

In order to explore the influence of systematic errors in the centering procedure we have introduced fixed shifts of 0.1 and 0.2 pixels after our centering procedure. We found that systematic centering errors as small as 0.1 pixels lead to a strong positive-negative pattern and the disruption of the polarization signal. These systematic errors produced similar pattern for the PSF comparison star as reported for the case of random shifts and can, therefore, not resemble the correct polarization pattern.

Additional tests included the comparison of centering methods based on cross-correlation and on Gauss fitting, bilinear and cubic interpolations at the rotation of the images, as well as the local and global normalization of the flat field images. We found that the Gauss-fitting centering, the bilinear interpolation and the individual normalization of each stripe in the Wollaston-masked flat field leads to the highest signal-to-noise ratios.

These test runs proved that the strongest noise sources in the reduction procedure are the subpixel misalignments and the flat field errors. We found no combination of the artificially introduced errors that could create a false (phase- and frequency-correct) sinusoidal modulation and thus could lead to the misinterpretation of our data set.

Chromatic effects from the Wollaston prism could also lead to slight differences between the o- and e-beams and thus to an imperfect subtraction. However, the MgF2 prism used in the CONICA camera introduces a relative displacement of only 86 mas between the blue and red edge of the broad K-band (CONICA Manual). In the narrower Ks-band this effect is even smaller and is therefore not expected to influence the PDI observations.

As TW Hya and the PSF star were observed in very similar airmass ranges (1.02 – 1.08 and 1.09 – 1.22) and seeing conditions (0.6″), we expect that any atmospheric effect leading to a false polarization signal would effect the PSF comparison images in the same way. Since there is no evidence for such an atmospheric influence in the comparison images (see right panel of Fig. 3), we exclude atmospheric chromatism as a possbility to mimick PDI signals.

The examination of the data reduction, possible instrumental and atmospheric effects led us to conclude that the extended polarization signal in our images indeed originate from the light scattering on the circumstellar dust around TW Hya.

## 5 Results

Based on the NACO observations, we obtained Ks-band images of the circumstellar disk around the classical T Tauri star TW Hya, which probe the disk structure closer to the star than any previous observation. Our diffraction limited PDI of TW Hya shows an extended butterfly pattern (see Figs. 1, 2) characteristic for a spatially resolved axisymmetric light scattering source. This pattern is present between 0.5″ and 1.4″ on the long exposure Q- and U-images and between 0.1″ and 0.4″ in the short exposure images. The observations of the regions closer than 0.1″ to the star and between 0.4″ and 0.5″ do not have sufficient SNR for a reliable analysis. The SNR of the data set is sufficiently high to plot the azimuthally integrated Ks-band surface brightness density distribution between 0.5″ and 1.4″ (see Fig. 4).

## 6 Discussion

### 6.1 Scattering

The presence of the extended polarization signal across our FOV confirms the existence of light scattering material. The sinusoidal modulation of the Q- and U-images confirms the centrosymmetric polarization pattern. The fact that we detected such emission even as close as $``$0.1″ ($``$6 AU) to TW Hya is the first direct evidence that its dust disk extends so far in. This detection does not support the model of Krist et al. (kri00 (2000)) who proposed the presence of an inner ’dark’ zone ($`r18`$ AU = 0.3″) with suppressed flaring angle, necessary to explain the exceptional brightness of the TW Hya disk.

Based on the modeling of the SED, Calvet et al. (cal02 (2002)) predicted the presence of a mostly dust-evacuated inner gap at $`r<34`$ AU ($`r<0.07^{\prime \prime }`$). Such a gap could be a tell-tale signpost of a massive planet in close orbit around TW Hya and would be of high importance. However, our commissioning data do not allow to draw any definite conclusions on the existence of this feature.

### 6.2 Radial Profile and Slope Change

In Fig. 4 we show the comparison of the radial polarized intensity derived from our PDI observation and (arbitrarily shifted) radial surface brightness profiles measured by the HST (Weinberger wei02 (2002)). Although these measurements probe different quantities at different wavelengths (Ks and J, H, respectively) a remarkable similarity is evident between the radial slopes. Based on this similarity, we conclude that the polarization degree is nearly independent of the radius and thus the polarized intensity characterizes the surface brightness. Thus, the PDI method can lead to contrast enhancement without the loss of information.

The radial surface brightness profile between $`0.9\mathrm{}`$ and $`1.4\mathrm{}`$ can be well fitted by a power-law function $`I_{\mathrm{Ks},1}r^{3.1\pm 0.3}`$. This is slightly steeper than the $`I_{\mathrm{J},\mathrm{H}}r^{2.6\pm 0.1}`$ behaviour found by Weinberger et al. (wei02 (2002)) at 1.1 and 1.6 $`\mu `$m, for the region between 0.8″ and 2.7″.

To place our results in context, we compare them to surface brightness slopes predicted by simple analytical approximations (e.g. Whitney & Hartmann whi92 (1992)). These approximations consider only isotropic scattering and simple geometric effects and do not include wavelength-dependent terms. The predicted surface brightness power-law exponent for an optically thick, geometrically thin disk is $`3`$, while for a flaring disk an exponent of -2 is expected. Our Ks-band results in the $`0.9\mathrm{}<r<1.4\mathrm{}`$ region support the flat disk geometry or a very small flaring. However, the surface brightness slope is determined by an intricate interplay between many disk and dust parameters and is not a conclusive measure of the disk geometry.

As to underline this statement, a strong, gradual change in the brightness profile around 0.8″ from TW Hya is evident. The annulus between $`0.5\mathrm{}<r<0.7\mathrm{}`$ is characterized by the very gentle surface brightness behaviour of $`I_{\mathrm{Ks},2}r^{0.9\pm 0.2}`$. The slope of the curve is even less steep than that expected from a flared disk (Whitney & Hartmann whi92 (1992)).

Although the radial brightness slope is strongly influenced by the disk geometry, slope changes can have different origins. Even though a detailed analysis of the reasons of the slope changes is beyond the scope of the current paper, we refer to some of them. The most obvious one is an actual change in the surface density slope of the disk, translating into a slope change in the surface brightness. Such clear-cut changes might be the results of external or internal perturbation of the disk by a companion or a massive planet as seen in the case of HR 4796A (see, e.g. Boccaletti et al. boc03 (2003)).

A second cause might be the change of the flaring angle and through this the change in the illumination of the disk’s surface. Alternatively, a radial dependence of the scattering cross section of the dust grains might introduce modulations in the surface brightness profile. Disentangling the effects of these processes requires high-angular resolution observations at longer wavelengths.

We note, that this change in the slope at 0.8″ has been detected previously at shorter wavelengths by HST direct imaging and coronographic studies (Krist et al. kri00 (2000), Weinberger et al. wei02 (2002)), together with other slope changes in the outer parts ($`>`$1.4″) of the disk.

Our data were acquired during the commissioning of the NACO instrument and suffer from high detector noise and limited AO performance. The current performance of NACO is expected to result in roughly five times better SNR.

The first results of a PDI survey of the TW Hya Association are presented by Huélamo & Brandner (hue03 (2003)).

For future PDI observers it might be important to point out some of the limiting factors in our data set. First, the dynamic range of the data set has to be increased by applying series of different exposure times. Furthermore, the flat field calibration plays a critical role and should be carried out with great care. Finally, the neccessarily saturated star will pose limitations to the accurate alignment of the individual frames, which will in practice determine the inner radius until which the polarization signal can be detected.

## 7 Summary

We presented the first high-contrast observations carried out with NACO/VLT using the PDI technique. The extended scattering pattern of TW Hya provides the first direct proof that the disk extends closer than 0.5″, up to $``$ 0.1″ ($``$ 6 AU) from the star.

We derive the first Ks-band radial polarized intensity profile between 0.5″ and 1.4″. This profile strongly resembles the surface brightness profile seen at shorter wavelengths and shows as well a strong, gradual slope change around 0.8″. The polarized intensity profile between 0.9″ and 1.4″ stands close to that expected from a flat disk.

These results demonstrate the potential of the PDI technique at the NACO/VLT platform for future planet formation studies. This technique is capable of imaging a new, beforehand unaccessible regime of circumstellar disks at very high resolution. Such images can help understanding the evolution of protoplanetary disks at the scales comparable to our inner Solar System.

###### Acknowledgements.

We acknowledge the outstanding support of the ESO staff both in Garching and at Paranal to the construction and commissioning of NACO. We thank O. Schütz for his useful comments on the manuscript. The suggestions of the referees (R. Racine and P. Bastien) helped to improve the clarity of the paper.


format: text

text:
# Paramagnetic anisotropic magnetoresistance in thin films of SrRuO3

## I Introduction

The phenomenon of anisotropic magnetoresistance (AMR) in magnetic conductors expresses the dependence of the resistivity $`\rho `$ on the angle $`\delta `$ between the current $`𝐉`$ and the magnetization $`𝐌`$. In polycrystals the AMR effect is commonly found to follow: $`\rho =\rho _{}+(\rho _{}\rho _{})\mathrm{cos}^2\delta `$ where $`\rho _{}`$ is the resistivity when $`𝐉𝐌`$, and $`\rho _{}`$ is the resistivity when $`𝐉𝐌`$ AMR . This simple relation is not expected to hold in crystalline samples where both the current orientation relative to the lattice as well as the magnetization orientation relative to the lattice play an important role.

Here we present AMR measurements of thin films of the $`4`$*d* itinerant ferromagnet SrRuO3 above $`T_c`$ ($`T_c150`$ K). Those films are epitaxial and characterized by large uniaxial magnetocrystalline anisotropy (MCA) easy axis . In our measurements, we study the AMR in uncommon conditions: (a) while in most AMR measurements the orientation of $`𝐌`$ is changed without changing its magnitude, here, because of the large MCA both orientation and magnitude of $`𝐌`$ are changing; and (b) while in most AMR measurement the applied field $`𝐇`$ is parallel to $`𝐌`$, here, because of the large MCA $`𝐌𝐇`$ except for the cases where $`𝐇`$ is along the easy or hard axes. Therefore, to explore the AMR in SrRuO3 it is not sufficient to measure magnetoresistance (MR) as a function of angle, but we need to independently determine both the magnitude and orientation of $`𝐌`$.

Our films are grown on miscut ($`2^{}`$) SrTiO3 substrates using reactive electron beam epitaxy. The films have orthorhombic structure ($`a=5.53`$ Å, $`b=5.57`$ Å, $`c=7.85`$ Å) and they grow uniformly (without twinning) with the $`c`$ axis in the film plane and the $`a`$ and $`b`$ axes at $`45^{}`$ out of the plane (see Fig. 1) Marshall . The MCA is uniaxial with the easy axis along the $`b`$ axis and thus the uniform growth of the measured films can be confirmed by MR measurements in the ferromagnetic phase angle . The films were patterned by photolithography to allow Hall effect (HE) and MR measurements. The sample whose results are presented here is $`300`$ Å thick with $`T_c`$ $`147`$ K and resistivity ratio of $`13`$.

## II Measurements and discussion

Our measurements have two parts: (a) extraordinary Hall effect (EHE) EHE  measurements from which we extract both the magnitude and the orientation of $`𝐌`$, and (b) MR measurements for currents in the $`[001]`$ and $`[1\overline{1}0]`$ directions. Combining the two measurements we show that the MR can be fit very well with a second order magnetization expansion.

Both HE and MR measurements were performed as a function of the angle $`\varphi `$ between the applied field $`𝐇`$ and the easy axis, where $`𝐇`$ is rotating in the $`(001)`$ plane. Each film has two kinds of patterns: a pattern with current along the $`[1\overline{1}0]`$ direction (denoted Pab) and a pattern with current along the $`[001]`$ direction (denoted Pc). While in Pab the angle between $`𝐉`$ and $`𝐇`$ varies with $`\varphi `$, in Pc the field $`𝐇`$ is always perpendicular to $`𝐉`$. The measurement configuration is illustrated in Fig. 1.

The Hall field in magnetic conductors has two contributions:

$$𝐄_H=R_0𝐉\times 𝐁R_s𝐉\times \mu _0𝐌$$

where $`𝐁`$ is the magnetic field, and $`R_0`$ and $`R_s`$ are the ordinary and the extraordinary Hall coefficients, respectively. By measuring the HE in our films at a temperature where $`R_s`$ vanishes rs  we determined $`R_0`$, which enabled to extract $`\mu _0R_sM_{}`$ (where $`M_{}`$ is the component of $`𝐌`$ which is perpendicular to the film plane) at all temperatures. This is not sufficient, however, since we need to determine both components of $`𝐌`$. For that we note that based on symmetry considerations we may assume that if a field $`𝐇`$, that is applied in the $`(001)`$ plane at an angle $`\varphi `$ relative to the easy axis, creates a magnetization $`𝐌`$ pointing at an angle $`\alpha `$ relative to the easy axis, then applying the same field at an angle $`\varphi `$ will create the same magnetization, but at an angle $`\alpha `$ demag . In our case the easy axis is at $`45^{}`$ out of the plane thus symmetry considerations yield that $`M_{}(\varphi )=M_{}(\varphi )`$ where $`M_{}`$ and $`M_{}`$ are the in-plane and perpendicular components of $`𝐌`$, respectively. Consequently, by measuring the Hall resistivity ($`\rho _{EHE}`$) at $`\varphi `$ and $`\varphi `$ we obtain: $`\rho _{EHE}\left(\varphi \right)=\mu _0R_sM_{}(\varphi )`$ and $`\rho _{EHE}\left(\varphi \right)=\mu _0R_sM_{}(\varphi )`$, which allows to determine $`𝐌`$ (multiplied by $`\mu _0R_s`$). Figures 2a and 2b show the change in the magnitude and direction of $`𝐌`$ as a function of $`\varphi `$ at $`T=170`$ K determined with that method. As expected, $`𝐌`$ obtains its maximum value at $`\varphi =0`$ and lags behind $`𝐇`$ except for $`𝐇`$ along $`a`$ and $`b`$. It seems that this is the first time that EHE which is sensitive only to the perpendicular component of the magnetization is used for extracting the full magnetization vector based on symmetry consideration. This is possible only due to the tilted easy axis. Where the easy axis is perpendicular or parallel to the film this scheme is not applicable.

Figures 2(c) and 2(d) present the MR measured at $`T=170`$ K for $`H=6`$ T and $`8`$ T in Pab and Pc. To fit the MR data we expand MR in $`𝐌`$, noting that due to the MR symmetry under field inversion the lowest order expansion is of second order. Since in our experiment $`𝐌`$ remains in the $`(001)`$ plane it is sufficient to use two components of $`𝐌`$. We use the freedom of choosing the principle axes to take them in the crystallographic directions of $`a`$ and $`b`$. Therefore, the general expansion of the MR to lowest order is:

$$MR=(\rho (H)\rho (0))/\rho (0)=A(M_b^2+\beta M_a^2+\gamma M_bM_a)$$
(1)
where $`M_b`$ and $`M_a`$ are the components of $`𝐌`$ along the easy axis ($`b`$) and hard axis ($`a`$), respectively. The lines in Figs. 2(c) and 2(d) are fits of the MR data based on the measured $`𝐌`$ and Eq. 1.

We obtained the following fitting parameters at $`T=170`$ K: $`\beta =1.4\pm 0.2`$, $`\gamma =0.6\pm 0.1`$ for the Pab pattern, and $`\beta =1.7\pm 0.2`$, $`\gamma =0.1\pm 0.1`$ for the Pc pattern. The error limits indicate the evaluated changes in the fitting parameters in case there is difference between the instrumental $`\varphi `$ and the actual $`\varphi `$ of up to $`2^{}`$.

The fits of the MR data based on Eq. 1 allow us to determine the ”clean” AMR effect; namely, how would the resistivity change if we could rotate $`𝐌`$ in the $`(001)`$ plane without changing its magnitude. Figure 3 shows the expected behavior at $`T=170`$ K for Pab and Pc for a value of magnetization obtained with $`H=8`$ T along the easy axis.

For Pc we note that although $`𝐉𝐌`$ there is a significant AMR ($`\beta 1`$), and that $`\gamma =0`$ is within our experimental accuracy, as could be expected from symmetry considerations. These results exhibit strong dependence of the MR not only on the angle between $`𝐌`$ and $`𝐉`$ (which remains constant) but also on the direction of $`𝐌`$ relative to the crystal. For Pab we note that the extremum values are not obtained for $`𝐉𝐌`$ or $`𝐉𝐌`$ but at intermediate angles. In fact, the extremum values for Pab are in between those obtained in Pc (along the $`a`$ and $`b`$ axes) and those observed in polycrystals (parallel and perpendicular to $`𝐉`$). This shows that in our case the AMR related to the orientation of $`𝐌`$ with respect to the lattice is of the same order of magnitude as the AMR due to the relative orientation of $`𝐌`$ and $`𝐉`$. We also note that for $`\alpha =45^{}`$, which corresponds to $`𝐌`$ perpendicular to the plane the MR is different in Pab and Pc despite the fact that in both cases $`𝐉𝐌`$. This illustrates the dependence of MR on the direction of $`𝐉`$ relative to the crystal.

In conclusion, we have presented AMR investigation of SrRuO3 with simultaneous measurements of $`𝐌`$ and MR in the same pattern thus enabling accurate determination of its AMR behavior despite the change in the magnitude of $`𝐌`$ and in its relative angle with $`𝐇`$. The results indicate significant AMR even in the paramagnetic state, where $`𝐌`$ is relatively small, and large effect of the orientation of $`𝐌`$ and $`𝐉`$ relative to the crystal axes.

We acknowledge support by the Israel Science Foundation founded by the Israel Academy of Sciences and Humanities.


text:
# An Automorphism of A Free Metabelian Group Without Fixed Points

## 1. Introduction

The study of free metabelian groups began in the 30-es with the work of Magnus . Later Bachmuth wrote a series of papers  devoted to automorphisms of the free metabelian groups. In  Shpilrain studied the fixed points of the endomorphisms of free metabelian groups. There, he raised the question whether any IA-automorphism (i.e., an automorphism which acts trivially on the abelianization) of the a metabelian group has a nontrivial fixed point. In this paper we construct explicitly an IA-automorphism without fixed points when the number of generators is at least $`3`$, thus giving negative answer to Shpilrain’s question. Such an automorphism does not exist in the two-generator case by a result of Bachmuth.

In the next section we have summarized several known results about free metabelian groups. Section 3 describes an automorphism without fixed points.

The author thanks his adviser E. Zelmanov for useful discussions during the preparation of the manuscript, and the referee for suggesting significant refinements of remark 3.8.

## 2. Free Metabelian Group

Let $`F_n`$ be the free group on the generators $`g_1,\mathrm{},g_n`$. Let $`F_n^{}=[F_n,F_n]`$ be its commutator subgroup, and let $`F_n^{\prime \prime }=[F_n^{},F_n^{}]`$ be the second commutator subgroup. The group $`M_n=F_n/F_n^{\prime \prime }`$ is called the free metabelian group of rank $`n`$. We will denote by $`g_i`$ the generators of $`M_n`$.

Let $`\varphi `$ be the two dimensional representation of the group $`M_n`$ defined by

$$\varphi (g_i)=\left(\begin{array}{cc}s_i& t_i\\ 0& 1\end{array}\right).$$

The coefficients of these matrices lie in the ring $`[s_i^{\pm 1},t_i]`$ obtained from $``$ by adjoining the commuting variables $`s_i`$ and $`t_i`$ together with the inverses of $`s_i`$. This representation is known as the Magnus representation. Magnus  showed that:

###### Theorem 2.1.

The representation $`\varphi `$ is faithful.

The following lemma from  describes the image of the representation $`\varphi `$.

###### Lemma 2.2.

Let $`S=s_1^{j_1}\mathrm{}s_n^{j_n}`$, where $`j_i`$ are arbitrary integers. Let $`\gamma _i`$ be polynomials in $`s_i^{\pm 1}`$. Then the matrix

$$\left(\begin{array}{cc}S& \gamma _it_i\\ 0& 1\end{array}\right)$$

lies in the image of $`\varphi `$ if and only if $`\gamma _i`$ satisfy the identity

$$\gamma _i(1s_i)=1S.$$

Also, every element in $`\mathrm{Im}(\varphi )`$ is a matrix of the above type.

The following lemma (also from ) gives a representation of the IA-endomorphisms by matrices over the ring $`[s_i^{\pm 1}]`$:

###### Lemma 2.3.

Let $`\alpha `$ be an IA-endomorphism of the group $`M_n`$. There exist polynomials $`\alpha _{i,j}`$ on $`s_i^{\pm 1}`$ such that

$$\varphi (\alpha (g))=\overline{\alpha }(\varphi (g)),$$

where $`\overline{\alpha }`$ is the endomorphism of the coefficient ring $`[s_i^{\pm 1},t_i]`$ of the representation $`\varphi `$ defined by

$$\begin{array}{cc}\hfill \overline{\alpha }(t_i)=& \alpha _{i,j}t_j\hfill \\ \hfill \overline{\alpha }(s_i)=& s_i.\hfill \end{array}$$

## 3. Construction of an Automorphism

In this section we give an explicit example of an automorphism of the group $`M_n`$, which has no fixed points.

###### Lemma 3.1.

Let $`\alpha _n`$ be an endomorphism of $`M_n`$ for $`n3`$ defined by

$$\begin{array}{cc}\alpha _n(g_i)=[[g_1,g_2]g_n,g_i]g_i\hfill & \text{if }i<n\hfill \\ \alpha _n(g_n)=[g_1,g_2]g_n.\hfill & \end{array}$$

Then $`\alpha _n`$ is an IA-automorphism of $`M_n`$. Here, $`[g,h]=ghg^1h^1`$ denotes the commutator of the elements $`g`$ and $`h`$.

###### Proof.

Let $`\beta _1`$ be the automorphism on $`M_n`$ which is the conjugation by the element $`g_n`$, i.e., $`\beta _1:gg_ngg_n^1`$; and let $`\beta _2`$ be the automorphism defined by $`\beta _2(g_i)=g_i`$ for all $`i<n`$ and $`\beta _2(g_n)=[g_1,g_2]g_n`$. The composition $`\beta _2\beta _1`$ acts on the generators as follows:

$$\begin{array}{cc}\beta _2\beta _1(g_i)=\beta _2([g_n,g_i]g_i)=[[g_1,g_2]g_n,g_i]g_i\hfill & \text{if }i<n\hfill \\ \beta _2\beta _1(g_n)=\beta _2(g_n)=[g_1,g_2]g_n\hfill & \text{for }i=n.\hfill \end{array}$$

We have $`\alpha _n=\beta _2\beta _1`$, and therefore $`\alpha _n`$ is an automorphism. The automorphism $`\alpha _n`$ acts trivially on the abelianization of $`M_n`$ because both $`\beta _1`$ and $`\beta _2`$ do. ∎

###### Lemma 3.2.

The automorphism $`\alpha _n`$ has no nontrivial fixed points in $`M_n`$.

###### Proof.

Let us give a direct proof using the Magnus representation.

Assume the contrary: let $`g1`$ be a fixed point of the automorphism $`\alpha _n`$. First, let us compute the action of the endomorphism $`\overline{\alpha }_n`$ (defined in lemma 2.3) corresponding to the automorphism $`\alpha _n`$.

###### Lemma 3.3.

The endomorphism $`\overline{\alpha }_n`$ acts on $`t_i`$ as follows:

$$\begin{array}{cc}\hfill \overline{\alpha }(t_i)=& s_nt_i+(1s_i)\left((1s_2)t_1(1s_1)t_2+t_n\right),\text{ for }in\hfill \\ \hfill \overline{\alpha }(t_n)=& (1s_2)t_1(1s_1)t_2+t_n.\hfill \end{array}$$

###### Proof.

We have that

$$\varphi ([g_1,g_2])=(\begin{array}{cc}s_1& t_1\\ 0& 1\end{array})(\begin{array}{cc}s_2& t_2\\ 0& 1\end{array})(\begin{array}{cc}s_1^1& s_1^1t_1\\ 0& 1\end{array})\left(\begin{array}{cc}s_2^1& s_2^1t_2\\ 0& 1\end{array}\right)$$

$$=\left(\begin{array}{cc}1& (1s_2)t_1(1s_1)t_2\\ 0& 1\end{array}\right).$$

Therefore,

$$\varphi (\alpha _n(g_n))=\left(\begin{array}{cc}s_n& (1s_2)t_1(1s_1)t_2+t_n\\ 0& 1\end{array}\right)$$

and for $`in`$

$$\varphi (\alpha _n(g_i))=\varphi (\alpha _n(g_n))\left(\begin{array}{cc}s_i& t_i\\ 0& 1\end{array}\right)\varphi (\alpha _n(g_n))^1=$$

$$=\left(\begin{array}{cc}s_i& s_nt_i+(1s_i)\left((1s_2)t_1(1s_1)t_2+t_n\right)\\ 0& 1\end{array}\right).$$

This, combined with the definition of $`\overline{\alpha }_n`$, proves the lemma. ∎

Let

$$\varphi (g)=\left(\begin{array}{cc}S& \gamma _it_i\\ 0& 1\end{array}\right),$$

where $`\gamma _i`$ are polynomials in $`s_i^{\pm 1}`$.

###### Lemma 3.4.

If $`\varphi (g)=\overline{\alpha }_n(\varphi (g))`$, then $`\gamma _i=0`$ for $`3i<n`$, and there exists $`A[s_i^{\pm 1}]`$, such that

$$\begin{array}{cc}\hfill \gamma _1=& (1s_2)s_nA,\hfill \\ \hfill \gamma _2=& (1s_1)s_nA,\hfill \\ \hfill \gamma _n=& (1s_n)A.\hfill \end{array}$$

###### Proof.

We have that

$$\begin{array}{c}\overline{\alpha }_n(\gamma _it_i)=\gamma _i\overline{\alpha }_n(t_i)=\hfill \\ =_{in}\gamma _i\left(s_nt_i+(1s_i)\left((1s_2)t_1(1s_1)t_2+t_n\right)\right)+\hfill \\ +\gamma _n\left((1s_2)s_nt_1(1s_1)s_nt_2+t_n\right).\hfill \end{array}$$

Equating the coefficients of $`t_i`$ in $`\gamma _it_i`$ and $`\overline{\alpha }_n(\gamma _it_i)`$, we have:

$$\begin{array}{ccc}\text{ coefficient at }t_1:\hfill & \hfill \gamma _1=& s_n\gamma _1+_{j<n}(1s_i)(1s_2)\gamma _k+(1s_2)s_n\gamma _n\hfill \\ \text{ coefficient at }t_2:\hfill & \hfill \gamma _2=& s_n\gamma _2_{j<n}(1s_i)(1s_1)\gamma _k(1s_1)s_n\gamma _n\hfill \\ \text{ coefficient at }t_i:\hfill & \hfill \gamma _i=& s_n\gamma _i,\text{ for }3i<n\hfill \\ \text{ coefficient at }t_n:\hfill & \hfill \gamma _n=& _{j<n}(1s_j)\gamma _j+\gamma _n.\hfill \end{array}$$

The element $`1s_n`$ is not a zero divisor in the ring $`Z[s_i^{\pm 1}]`$. Therefore, we have $`\gamma _i=0`$, for $`3i<n`$. Thus, the system is reduced to

$$\begin{array}{cc}\hfill (1s_n)\gamma _1=& (1s_2)\left(\left((1s_1)\gamma _1+(1s_2)\gamma _2\right)+s_n\gamma _n\right)\hfill \\ \hfill (1s_n)\gamma _2=& (1s_1)\left(\left((1s_1)\gamma _1+(1s_2)\gamma _2\right)+s_n\gamma _n\right)\hfill \\ \hfill 0=& (1s_1)\gamma _1+(1s_2)\gamma _2.\hfill \end{array}$$

The expression in the middle brackets of the first two equations is zero. Therefore, they are reduced to

$$\begin{array}{cc}\hfill (1s_n)\gamma _1=& (1s_2)s_n\gamma _n\hfill \\ \hfill (1s_n)\gamma _2=& (1s_1)s_n\gamma _n.\hfill \end{array}$$

The elements $`1s_i`$ are prime in the ring $`Z[s_i^{\pm 1}]`$. Therefore, $`(1s_n)`$ divides $`\gamma _n`$, $`(1s_1)`$ divides $`\gamma _2`$ and $`(1s_2)`$ divides $`\gamma _1`$, i.e., there exist elements $`A_i`$ such that $`\gamma _1=(1s_2)A_1`$, $`\gamma _2=(1s_1)A_2`$ and $`\gamma _n=(1s_n)A_n`$. From the system above it follows that $`A_1=A_2=s_nA_n`$, which completes the proof. ∎

###### Lemma 3.5.

If $`\varphi (g)=\overline{\alpha }_n(\varphi (g))`$, then $`\varphi (g)=1`$.

###### Proof.

The matrix $`\left(\begin{array}{cc}S& \gamma _it_i\\ 0& 1\end{array}\right)`$ is in the image of $`\varphi `$. By lemma 2.2 we have that $`1S=\gamma _i(1s_i)`$. Applying the previous lemma, we have

$$1S=(1s_n)^2A.$$

If $`S1`$, then the point $`(1,\mathrm{},1)`$ is a simple zero of the left hand side, but it is a double of the right hand side, a contradiction. Therefore, the only solution to the above equation is $`S=1`$ and $`A=0`$, i.e., $`\varphi (g)=1`$. ∎

This shows that any fixed point of the automorphism $`\alpha _n`$ lies in the kernel of $`\varphi `$. By lemma 2.1 the representation $`\varphi `$ is faithful. Therefore, the automorphism $`\alpha _n`$ has no nontrivial fixed points. ∎

Thus, we have shown:

###### Theorem 3.6.

There exists an automorphism of the free metabelian group on more than $`2`$ generators which acts trivially on the abelianization, and the only fixed point is the identity.

###### Remark 3.7.

Every IA-automorphism of the group $`M_2`$ has a fixed point. Bachmuth  has shown that every such automorphism is inner, and therefore, has a fixed point.

###### Remark 3.8.

For the construction of the automorphism $`\alpha _n`$ we were motivated by the Lie correspondence. It is well known that many problems regarding a Lie group $`G`$ can be transferred to similar problems about its Lie algebra.

In general, these two problems are not equivalent, because the exponential map is not a bijection. In the case of connected unipotent Lie groups, the exponential map is a bijection, and one can use the Lie correspondence to relate the question of existence of IA-automorphisms without fixed points to the existence of derivations of their Lie algebras with trivial kernel.

We could not apply this correspondence directly to our problem because $`M_n`$ is not a unipotent Lie group. One way to go around this problem is to embed $`M_n`$ into a pro-unipotent Lie group $`\stackrel{~}{M}_n`$, where $`\stackrel{~}{M}_n`$ is constructed as follows: Let $`A=x_1,\mathrm{},x_n`$ be the algebra of formal power series of the non-commuting variables $`x_i`$ and let $`A^{\prime \prime }`$ be its second commutator subalgebra. Finally, let $`\stackrel{~}{M}_n`$ be the closure of the subgroup of the multiplicative group of $`1+A/A^{\prime \prime }`$ generated by the images of $`1+x_i`$ together with the operations $`gg^\alpha `$ for any real number $`\alpha `$, where the operation is defined as a formal power series.

The Lie algebra corresponding to the pro-unipotent group $`\stackrel{~}{M}_n`$ is the free metabelian Lie algebra $`ML_n`$ on $`n`$ generators, i.e., the quotient of the free Lie algebra by its second commutator algebra. Our first step was to construct a derivation $`D_n`$ of $`ML_n`$ with a trivial kernel — we took $`D_n(x_i)=[x_i,x_n]`$ for $`i<n`$ and $`D_n(x_n)=[x_1,x_2]`$ and proved that it has a trivial kernel using a well known basis of the free metabelian algebra $`ML_n`$. The next step was to lift this derivation to the automorphisms of $`\stackrel{~}{M}_n`$ and $`M_n`$, i.e., to construct an automorphism $`\alpha _n`$ of $`M_n`$ such that its extension to $`\stackrel{~}{M}_n`$ is almost the same as $`\mathrm{exp}(D_n)`$.

In the case of the free metabelian group, the Lie correspondence is relatively weak — for example there exist derivations of $`ML_2`$ without kernel but they can not be used to construct automorphisms of the group $`M_2`$.

One way of showing that the automorphism $`\alpha _n`$ moves any nontrivial element $`gM_n`$ is to pass to a suitable nilpotent quotient of $`M_n`$ and use the Lie correspondence there. Although such a proof would be conceptually easier, the are many additional details that have to be proved, which would make it more technical. Therefore, we chose the more direct proof involving the Magnus representation of the free metabelian group.


text:
# A Revised Ephemeris and FUSE Observations of the Supersoft X-ray Source CAL 831footnote 11footnote 1This paper utilizes public domain data obtained by the MACHO Project, jointly funded by the US Department of Energy through the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48, by the National Science Foundation through the Center for Particle Astrophysics of the University of California under coopertative agreement AST-8809616, and by the Mount Stromlo and Siding Spring Observatory, part of the Australian National University. 2footnote 22footnote 2Based on observations made with the NASA-CNES-CSA Far Ultraviolet Spectroscopic Explorer. FUSE is operated for NASA by the Johns Hopkins University under NASA contract NAS5-3298.

## 1 Introduction

The X-ray binary CAL 83 is one of a small number of supersoft sources which are characterized by an extremely soft X-ray spectrum with almost no radiation above $`0.5`$ keV. This group includes a handful of sources in the Magellanic Clouds and the Galaxy (e.g. Greiner 1996, Hasinger 1996, Cowley et al. 1998), plus some less well studied objects in other nearby galaxies (e.g. Kahabka 1999; Swartz et al. 2002; Kong & Di Stefano 2003). CAL 83 itself is a member of the LMC, based on its position and radial velocity.

From radial velocity variations of the optical emission lines (especially He II and H), Crampton et al. (1987) discovered that CAL 83 is a binary system with a period near one day. The emission lines are formed in an accretion disk that is viewed at a low inclination angle, so the resulting velocity amplitude is small (K$`30`$ km s-1). There is no spectroscopic evidence of either star, but the period and velocities suggest that they have masses well under 2M$``$ (see Crampton et al. 1987, Smale et al. 1988).

CAL 83 is known to exhibit irregular brightness variations of more than a magnitude over timescales of months (e.g. Smale et al. 1988, Alcock et al. 1997, Greiner & Di Stefano 2002). Although CAL 83 is considered a nearly steady X-ray source, occasionally it drops into an X-ray-off state (e.g. Kahabka 1996, 1998). Using MACHO data, Alcock et al. and Greiner & Di Stefano have shown that these low X-ray states are followed within $`3050`$ days by a low optical state. The most likely explanation is that there is a change in the accretion flow and hence in the accretion disk which provides the optical luminosity.

In spite of these large changes in the mean brightness of CAL 83, Smale et al. were able to detect small superimposed orbital variations with a full amplitude of $``$0.2 mag and a period of P = 1.0436 days. During the orbital period the light curve shows a sinusoidal variation with one maximum and one minimum per cycle. Cowley et al. (1991) later refined the orbital period to be P = 1.0475 days by combining all published and new spectroscopic data as well as the photometry of Smale et al.

Given the low inclination of the system, it is unlikely that the minimum is caused by an occulatation by the secondary star. Furthermore, the phasing of the emission line velocities implies that the light minimum does not coincide with conjunction of the stars. The variation is likely to arise from a partial obscuration of the bright inner accretion disk. Possible occulting regions might be a splash region raised by impact of the accretion stream or a gas stream which is significantly above the orbital plane. Although we are uncertain what causes the orbital light curve, we have been able to use it to refine the orbital period to allow intercomparison of two decades of spectroscopic data.

## 2 Analysis of Optical Data

### 2.1 Calibration of the MACHO data

CAL 83 was observed for $``$7 years by the MACHO Project. We extracted both the “blue” and “red” data for this star using the interactive online web browser333See http://www.macho.mcmaster.ca/Data/MachoData.html. After transformation of the instrumental magnitudes, using the method described by Alcock et al. (1999), the derived values are close to $`V`$ and $`R`$ band magnitudes. We refer to these as $`V_{MAC}`$ and $`R_{MAC}`$, as we have done in previous papers. The linear transformation equations adopted here differ slightly from those given by Greiner & Di Stefano (2002), in the sense that our $`V_{MAC}`$ and $`R_{MAC}`$ values are 0.012 mag brighter and 0.013 mag fainter, respectively, than their “generic” photometry. The combination of offsets makes our $`V_{MAC}R_{MAC}`$ colors 0.025 magnitudes bluer. All of these systematic differences, however, are small compared to the estimated zero-point calibration accuracies ($`\pm `$0.035 mag) given by Alcock et al. Our primary concern is with the internal errors of the MACHO data. Hence we propogate, in all of our calculations, only the photometric uncertainties of the instrumental magnitudes. In applying the transformations it was assumed that the $`V_{MAC}R_{MAC}`$ color for CAL 83 is constant, but small variations are present in the MACHO data (see $`\mathrm{\S }`$2.2).

Figure 1 shows the $`V_{MAC}`$ light curve for CAL 83. The source is highly variable with rapid changes in brightness between high, intermediate, and low optical states (see Greiner & Di Stefano). There were two epochs when the mean magnitude remained relatively stable (MJD 48917–49640 and 50980–51541). CAL 83 was in an intermediate state during Segment 1 and a high state during Segment 2.

### 2.2 Color Variations Based on MACHO Data

In order to understand what causes the light variations, it is of interest to determine whether or not the color changes as the system brightens. Greiner & Di Stefano (2002), in their study of the long-term on/off states, found no relationship between color and brightness. They pointed out that the reddest $`VR`$ data have very large errors and argued that the mean $`VR`$ color varies by no more than 0.03$`\pm `$0.05 mag for $`>`$1 mag changes in system brightness (see their Fig. 2). However, we find that the mean $`V_{MAC}R_{MAC}`$ depends on CAL 83’s optical state. In Figure 2 we plot $`V_{MAC}`$ as a function of $`V_{MAC}R_{MAC}`$. The dividing lines between high, intermediate, and low states roughly correspond to the limits set by Greiner & Di Stefano. One sees that the color is bluest when the system is in its high state, although the change is small.

Removing the 17 points in the figure with the largest uncertainties (i.e., error estimates $``$0.048 mag), the remaining points have a distribution of errors that approximates the expectation for a Gaussian distribution with $`\sigma 0.018`$ mag. That is, $``$68% of the remaining points have errors $``$1$`\sigma `$ and $``$97% have errors $``$2$`\sigma `$. This close agreement with the expected percentages for a normal distribution (i.e., 68%, 95%, etc.) validates our selection of points to remove from the analysis. We find the mean $`V_{MAC}R_{MAC}`$ colors are $`0.065\pm 0.019`$ (56 points), $`0.000\pm 0.038`$ (66 points), and $`0.024\pm 0.027`$ (31 points) for the high, intermediate, and low states, respectively. The observed dispersion within each sample is comparable to the internal precision of the entire MACHO dataset ($`\sigma _{VR}=0.028`$ mag for $`V18`$ mag) as demonstrated by Alcock et al. (1999). Simple T-tests show that the samples are significantly different compared to each other. However, systematic errors are present in the means because we have assumed constant color terms in the transformation equations. Hence, we have slightly overestimated the color differences between states. The largest mean difference (between high and intermediate states) is $``$0.065 mag, but after compensating for the transformation color terms, the revised differential is $``$0.048 mag. The corrected value is still much greater than the expected internal precision of MACHO data (i.e., 0.028 mag) and therefore judged to a meaningful difference. All of the other corrections are smaller than this example. After making the adjustments, the color difference of even the smallest differential (between intermediate and low states) is still significant at a 98% level.

### 2.3 Determination of a Revised Orbital Ephemeris from MACHO Data

After removal of a few errant (i.e., spuriously faint) points, data from the two time segments marked in Figure 1 were searched for orbital variations. Since the data are more complete for $`V_{MAC}`$, observations in only that bandpass were used to redetermine the orbital period. The $`V_{MAC}`$ magnitudes for each segment were analyzed separately to derive periodograms, using the technique described by Horne & Baliunas (1986). Both segments contain one significant peak with very similar values near P = 1.04 days. (Due to the daily sampling window, several aliases also are present.) Therefore, we combined the data from the two segments after subtracting the mean linear fit to each segment. This resulted in a single set of relative magnitudes, $`V^{}`$, from which the combined power spectrum, shown in Figure 3, was determined. Fitting a sinewave to the $`V^{}`$ data yielded the following ephemeris:

T(min) = JD 2451500.953$`\pm `$0.004 + 1.047568$`\pm `$0.000003E days

For the epoch of the FUSE observations presented in this paper (see $`\mathrm{\S }`$3), the new ephemeris gives phases 0.15 earlier than those predicted by Smale et al. (1988), so this revision of the orbital period is important for correctly phasing the FUSE data.

Plots of the $`V^{}`$ and $`R^{}`$ light curves, folded on the new ephemeris, are shown in Figure 4, with points from the two time segments indicated by different symbols. We note that the amplitude of the orbital variation appears to depend on the mean brightness level in both bandpasses. The $`V_{MAC}`$ data show that when the system was at $`V_{MAC}`$ = 17.34 during the first segment, the semiamplitude of the variation was 0.053$`\pm `$0.001 mag, but during the second segment (mean magnitude $`V_{MAC}`$ = 16.97) the semiamplitude was 0.075$`\pm `$0.001 mag. Similarly, the $`R_{MAC}`$ data show that the amplitude in the first (fainter) time segment was 0.050$`\pm `$0.003, while in the last (brighter) segment the amplitude was 0.082$`\pm `$0.002. In addition, the data of Smale et al. have a semiamplitude of 0.11 mag when the system had a mean magnitude of $`V`$ = 16.87. Hence, the amplitude of the orbital modulation appears to increase as the system brightens.

Smale et al. (1988) do not tabulate their photometric data, so we are unable to plot their data on the MACHO ephemeris. However, using Smale et al.’s ephemeris, we estimate minimum light occurred near HJD 2446051.585, which corresponds to our phase 0.07. The discrepancy is small but reasonable, given that their assumed time of minimum is based on an analysis of data which are folded on a somewhat different orbital period and may contain fluctuations in mean brightness.

In the bottom panel of Figure 4 we plot $`V_{MAC}R_{MAC}`$ versus orbital phase. Although there is no significant variation of color with phase, the figure confirms that the mean color in the second segment (when the source was in a high optical state) is bluer than in the first segment (intermediate state), as found in $`\mathrm{\S }`$2.2.

### 2.4 Previously Unpublished Photometry of CAL 83

$`B`$ and $`V`$ photometry was obtained on 8 nights in 1999 March using the TEK2K #3 CCD on the CTIO 0.9 m telescope. The images were reduced with DAOPHOT (Stetson 1987) and calibrated via observations of Landolt’s (1992) standard stars. Differential photometry was then calculated using a set of reference stars in the field of view of CAL 83. The same comparison stars were used in our previous photometric studies (Crampton et al. 1987, Cowley et al. 1998), hence all of our measurements are on a common scale. The new photometric measures are listed in Table 1 and shown in Figure 5. For a short observing run, the $``$1 day orbital period results in poor phase coverage, so it is not possible to use these data to determine an independent value for the period. Hence, the $`V`$ light curve in the top panel of the figure is fitted using a sinewave with P and T(min) fixed at the new values derived in $`\mathrm{\S }`$2.3. The resulting semiamplitude (0.062$`\pm `$0.012 mag) and mean magnitude ($`V`$ = 16.986$`\pm `$0.009) are consistent with the MACHO reductions for the same epoch. The color, plotted in the bottom panel, shows little variation and has a mean value of $`BV`$ = $``$0.059$`\pm `$0.015.

## 3 Re-analysis of Published Spectroscopic Data Using the New Ephemeris

We have gone back to the spectroscopic data given by both Crampton et al. (1987) and by Smale et al. (1988). Using the He II (4686Å) velocities from Crampton et al., we have computed revised orbital elements for CAL 83 which are presented in Table 2. A plot of the He II velocities (from Crampton et al.) versus phase is shown in the lower panel of Figure 6, where it is compared to the FUSE O VI data. Although the He II velocities tabulated by Smale et al. are not plotted, their phasing and amplitude are very similar. We note that their velocities appear to be systematically higher by $``$30 km s-1. It is likely that this is related to the longterm variations in the shape of the He II line pointed out by Crampton et al. Over a timescale of months, the $`\lambda `$4686 line shows wings extending several thousand km s-1, alternately to longer and then shorter wavelengths. Thus, depending on the overall line shape during the epoch when observations were made, the mean He II wavelength might be measured somewhat differently, thus accounting for apparent shifts in the systemic velocity.

From HST data observed 1996 November 10 (our orbital phase $`\mathrm{\Phi }`$ = 0.69), Gänsicke et al. (1998) measured a velocity of +300 km s-1 for the N V doublet at $`\lambda \lambda `$1239 & 1243. This value consistent with our He II velocity curve. They note that the lines are asymmetrical with a redward wing extending to $``$+800 km s-1.

The emission lines are typical of those from an accretion disk which surrounds the compact star. Provided there is no strong contribution to these lines from the mass transfer stream, then they should trace the orbital motion of this star. The velocity curve shows that minimum light occurs about 0.12P before the compact star is at superior conjunction.

Data in the literature show that the mean He II equivalent widths (EW) changes over the long term, and probably with phase, although this information is incomplete for most epochs. Our 1984 January spectra from Las Campanas have a mean EW = $`5`$Å, while spectra from 1985 December when the system was faint ($`V17.3`$, Crampton et al. 1987) show nearly double that value. This suggests that when the system is faint the emission line EW may be greater, but further data are needed to confirm this.

## 4 FUSE Observations and Measurements

In 2001 September, CAL 83 was observed using FUSE. This instrument is sensitive over the wavelength range 912–1187Å, at a spectral resolution of R = 20000 or more. FUSE has four independent optical channels, feeding two detectors via Rowland circle spectrographs. Two optical trains use LiF-coated optics (990–1187Å) and two have SiC-coated optics for wavelengths down to 912Å. As each detector is in two sections, there are 8 channels of data, designated Lif1a, Lif1b, etc., through Sic2b. A more complete description of the instrument is given by Moos et al. (2000).

The observations of CAL 83 were taken during 11 consecutive FUSE orbits, with earth-occultation interruptions. The star was located in the large aperture, and the alignment of the FUSE channels should have been stable over the whole observation. The visibility times were generally of order 4000 ksec. The whole observing program covered about 0.6P of the binary orbit. The basic observational details are given in Table 3. The data were extracted and processed with CALFUSE, current in April 2002. Because the source is faint, the airglow spectrum constitutes a major contaminant and effectively obliterates a significant part of the FUSE spectral range.

The only stellar line feature seen in the FUSE spectrum is O VI emission at 1031.9Å. The second line of the O VI doublet, 1038Å, is overlaid by airglow lines with peak emission some ten times that of the O VI line. The 1032Å line is weak, so in order to measure the continuum and O VI with more reliability we have binned the observations into 5 spectra, using consecutive pairs except for the last 3 observations which were summed to form a single spectrum. Figure 7 shows the binned spectra in the region of the 1032Å O VI emission line, together with some nearby absorption features. These spectra combined the LiF channels only, as the SiC channels were very noisy.

The continuum flux was measured in the airglow-free wavelength region 1045-1080Å, in both the LiF1a and LiF2b channels. Similar changes were seen in both, but the LiF2b showed a much greater range, consistent with a poorly-determined zero level in the extraction. As the LiF1a channel is that used for guiding, we have adoped continuum values from this channel only. (The profiles shown in Figure 7 include the LiF2b to maximize the signal, but the O VI profile is not affected by zero point uncertainties.) The continuum flux values, given in Table 3, have a mean value of $`1.2\times 10^{14}`$ erg s-1 cm-2 Å-1 and are almost flat across the FUSE bandpass. Since the changes in the continuum are similar for the LiF1a and LiF2b channels, we consider them to be real and note that they have a range and phasing similar to the MACHO data. Figure 6 shows the FUSE continuum changes compared with the mean optical light curve.

The O VI line was measured for flux and radial velocity, and these values are also given in Table 3. We estimate the uncertainty in measurement of the line flux is no larger than $``$20% of its value. The O VI flux appears to be greatest between phases 0.5-0.8 and weakest near phase 0.1, although the phase coverage is incomplete. The O VI flux variation appears to be similar to that of the optical He II line, for the overlapping phases. Additional FUSE observations are needed to see if the observed O VI flux variations are phase related.

The O VI line was measured in several ways, first by profile fitting of line and then by using the centroid of the flux above the local continuum. In addition, the velocities were measured by cross-correlation against the mean CAL 83 spectrum. There is close agreement (in all cases $`<`$15 km s-1) between the profile fitting technique and the cross-correlation velocities. The adopted values given in Table 3 are the mean of the two methods. We estimate the error of each measure to be no more than $`\pm `$20 km s-1. Both sets of measured velocities are plotted versus phase in the bottom panel of Figure 6, with different symbols for each method to show the close agreement. In this figure we also show the He II velocities from Crampton et al. (1987).

As a check on the stability of the wavelength scale, the double absorption feature of C II at 1036–7Å was also measured by cross-correlation of only this part of the spectrum against the mean for all CAL 83 FUSE spectra. Although the C II lines are weak and hence difficult to measure, their velocities show no correlation with each other or with the O VI velocities. Even though in Figure 6 it may appear that their positions are not constant, noise in these very weak features makes the profiles and/or the doublet separation appear to change. The measurements of C II are consistent with no velocity variation, with a scatter of about $`\pm `$12 km s-1.

Finally, the wavelength scale itself was checked by measuring the airglow emissions near 1038Å. These airglow lines were found to have the same wavelength for all binned spectra, to within $`\pm `$0.01Å, but their average position is too high by 15 km s-1. This wavelength scale error is typical of large aperture observations, but may partly reflect the average FUSE orbital motion through the visibility windows to the target. We have applied a correction of +15 km s-1 to the O VI velocity measures.

We note that the mean O VI velocity of +434 km s-1 is much more positive than the systemic velocity of +257 km s-1 given here from He II, 4686Å emission. It is possible that the shortward edge of the line has P Cygni-like absorption (see Figure 7), which may account for some of the velocity offset. However, the velocity difference is more than 0.5Å, so this explanation would imply a much broader and stronger intrinsic emission profile and a large absorption in a disk wind. As noted above, the longterm changes in the emission wings of He II were reported by Crampton et al. (1987; also see discussion in Cowley et al. 1998). It is possible that the offset seen in the velocities from the FUSE data may be partially due to the long wavelength extension of the line profile, which can be seen in Figure 7.

Fitting the O VI FUSE velocities to a sinewave results in a semiamplitude K = 23$`\pm `$6 km s-1 with maximum velocity occuring at photometric phase $`\mathrm{\Phi }=0.23\pm 0.03`$. Thus, the phasing of the O VI emission appears to differ from He II, whose maximum velocity is at $`\mathrm{\Phi }=0.88`$ with a semiamplitude of K = 27 km s-1 (see Table 2). If further observations in the far ultraviolet confirm this difference, then the O VI emission may arise in a different location from the He II.

## 5 Summary

To summarize, we have used MACHO photometry to refine the orbital period of the LMC supersoft X-ray source CAL 83. With our new ephemeris, we have been able to compare previously published spectroscopic and photometric data using consistent phases. In a high optical state, its color becomes bluer and the amplitude of the orbital variations increases. New FUSE observations show that orbital variations in brightness in the far ultraviolet are similar to those in the optical region. However, the velocity of the O VI emission in the FUSE spectra indicates its behavior may differ from the optical He II line.

We acknowledge help from the CTIO staff when our photometric observations were made.


text:
# Alloy effects in Ga_{1-𝑥}In_𝑥N/GaN heterostructures

## Abstract

We show that the large band offsets between GaN and InN and the heavy carrier effective masses preclude the use of the Virtual Crystal Approximation to describe the electronic structure of $`\mathrm{Ga}_{1x}\mathrm{In}_x\mathrm{N}/\mathrm{GaN}`$ heterostructures while this approximation works very well for the $`\mathrm{Ga}_{1x}\mathrm{In}_x\mathrm{As}/\mathrm{GaAs}`$ heterostructures.

The use of semiconductor alloys proves to be necessary to adjust the electronic parameters to the design of specific devices. $`\mathrm{Ga}_{1x}\mathrm{In}_x\mathrm{N}`$ alloys inserted between GaN Adelmann  would be ideally suited to cover the entire near infrared/ultra-violet spectrum Nakamura1 . Despite the lack for translation invariance, one may very often describe the electronic behavior of an alloy by means of a band structure, the so-called Virtual Crystal Approximation Harrison  (VCA). Close to it, the Coherent Potential Approximation allows to introduce part of the disorder effects in the VCA scheme, which leads to a damping of the VCA Bloch states. An example of alloys which is well described by the VCA/CPA is $`\mathrm{Ga}_{1x}\mathrm{In}_x\mathrm{As}`$. The VCA has also been used to describe the electronic properties of $`\mathrm{Ga}(\mathrm{In})\mathrm{N}`$ heterostructures Hangleiter

In this letter, we shall show that the $`\mathrm{Ga}_{1x}\mathrm{In}_x\mathrm{N}`$ system behaves in a radically different way. A numerical computation of the electronic states of $`\mathrm{Ga}_{1x}\mathrm{In}_x\mathrm{N}/\mathrm{GaN}`$ quantum wells and quantum dots will show very large differences from the VCA predictions. Instead, we shall show that $`\mathrm{Ga}_{1x}\mathrm{In}_x\mathrm{As}/\mathrm{GaAs}`$ electronic states are very close to the VCA predictions. The difference between the two alloy families arises from the much larger band offsets between GaN and InN than those between GaAs and InAs and from the heavier effective masses in the nitride system than in the arsenide system. We note that the large potential offset between $`\mathrm{In}`$ and $`\mathrm{Ga}`$ has led Kent and Zunger Kent  to predict that nitride based alloys cannot be described by models which neglect fluctuations.

In the following, we use one band effective mass Hamiltonians to describe the electron and hole kinematics. The effective mass is slightly anisotropic in the conduction band ($`m_{ez}=0.184m_0,m_{e//}=0.166m_0`$) and anisotropic in the valence band ($`m_{hz}=1.1m_0,m_{h//}=0.5045m_0`$) Morel1 ; Morel2 . The conduction $`\mathrm{\Delta }E_c`$ (valence $`\mathrm{\Delta }E_v`$) band offsets between GaN and InN is taken equal to 1.8 eV (0.9 eV) while those between GaAs and InAs are taken equal to 0.41 eV and 0.29 eV respectively Stier . Thus, there is a much larger energy fluctuation in the nitrides than in the arsenides when in a given unit cell a Ga atom replaces an In atom. The two heterostructures we shall be dealing with are a 3.2 nm thick $`\mathrm{Ga}_{0.83}\mathrm{In}_{0.17}\mathrm{N}/\mathrm{GaN}`$ quantum well and a $`\mathrm{Ga}_{0.83}\mathrm{In}_{0.17}\mathrm{N}`$ truncated pyramid with hexagonal basis (6 nm side), 2.6 nm height and a basis angle of $`30^{}`$ embedded into a GaN matrix and floating on a 1.1 nm thick wetting layer (see Fig.(1)). A pyramid with the same geometrical parameters made of $`\mathrm{Ga}_{0.5}\mathrm{In}_{0.5}\mathrm{As}`$ and embedded into GaAs will be considered for comparison (50% of In instead of 17% was considered to get bound states to the pyramid). It is well known that nitride heterostructures contain huge internal fields Morel1 ; Morel2 ; Bernardini . This will be modeled by assuming that there exist piecewise constant electric fields which are oriented along the $`Z`$ direction and have a magnitude of 2.45 MV/cm in $`\mathrm{Ga}_{0.83}\mathrm{In}_{0.17}\mathrm{N}`$ and - 0.1 MV/cm in GaN Morel1 ; Morel2 .

The quantum well problem will be reduced to a two dimensional (2D) calculation by considering an effective in - plane potential which is the average of the actual 3D potential weighted by $`\left|\mathrm{\Psi }(Z)\right|^2`$ where $`\mathrm{\Psi }\left(Z\right)`$ is the envelope function solution of the 1D problem in the presence of the electric field and a square well confining potential with a depth of 0.31 eV and 0.15 eV ($`=0.17\mathrm{\Delta }E_c`$ or $`0.17\mathrm{\Delta }E_v`$) (see Fig.(1)). This effective 2D potential is then supplemented by a hard box potential which confines the carrier in the layer plane into a square with dimensions $`L\times L=36\mathrm{nm}\times 36\mathrm{nm}`$. A given sample is then randomly generated by filling each unit cell either by In or by Ga corresponding to a potential energy of either 0 (In) or $`\mathrm{\Delta }E_c`$ (Ga) for electrons and either 0 (In) or $`\mathrm{\Delta }E_v`$ (Ga) for holes. For the quantum dots a complete 3D diagonalization was undertaken enclosing the pyramid in a box with a square basis with a side of 14.9 nm and a height of 7.9 nm. The box is then subdivided into smaller boxes of atomic size where the potential is supposed to be constant. The wave functions are then developed on the sinus solutions leading to a Hilbert space of dimension 75000.

For the QDs, two situations have been considered : either uncorrelated sites where the probability that a given cell is occupied by In is $`x`$ ( = 17%) irrespective of the occupancy of the other sites or site correlations which favor the In cluster formation Narakawa ; Donnell  by allowing a probability $`p>x`$ for a cell to be occupied by an In if its first neighborings are already occupied by at least one In. The only samples retained in the analysis are those such that the In fraction in the whole QD lies in the interval $`[0.17,0.18]`$.

The first few eigenstates of those 2D (40 states) and 3D (6 states) problems are evaluated by means of exact diagonalizations using the Lánczos algorithm. Once the eigenstates are known one can calculate various averages like the mean position, the mean square deviation to it, etc..

While the properties of a single sample have virtually no chance to be compared to any experiment, the average properties, obtained by taking the arithmetic mean of the outcomes of a given physical quantity in a series of samples, have more physical substance since they can (in principle) be compared to the results obtained on ensembles of quantum dots or quantum wells. This is for these average properties that one can discuss the effectiveness or failure of the VCA. One important quantity is $`D_e`$, the average electron density of states which is :

$`D_e(ϵ)={\displaystyle \frac{1}{\pi N}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{i}{}}{\displaystyle \frac{\mathrm{\Gamma }/2}{\left(ϵϵ_{ei}\left(j\right)\right)^2+\left(\mathrm{\Gamma }/2\right)^2}}`$ (1)
where N is the number of samples (which differ from one another by the locations of the In atoms) while $`ϵ_{ei}(j)`$ is the $`i^{th}`$ eigenvalue of the $`j^{th}`$ sample and $`\mathrm{\Gamma }`$ has been taken equal to 1 meV. Here we replace the usual delta peaks by Lorentzian distributions for clarity. Equally important (and measurable) is the average electron - hole optical density of states defined as :

$`D_{eh}(ϵ)={\displaystyle \frac{1}{\pi N}}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{k,l}{}}{\displaystyle \frac{\left(\mathrm{\Gamma }/2\right)\left|\mathrm{\Psi }_{ek}|\mathrm{\Psi }_{hl}\right|^2}{\left(ϵϵ_{ek}(j)ϵ_{hl}(j)\right)^2+\left(\mathrm{\Gamma }/2\right)^2}}`$ (2)
In the finite QW structure the VCA results in peaks related to the confinement box which fulfill :

$`ϵ_{np}=V_0+ϵ_{1z}+{\displaystyle \frac{\mathrm{}^2\pi ^2}{2m_{//}L^2}}\left(n^2+p^2\right)`$ ; $`n,p=1,2,\mathrm{}`$ (3)
where $`V_0=x\mathrm{\Delta }E_c`$ is the VCA average potential, $`ϵ_{1z}`$ is the lowest confining energy in the Z direction and L = 36 nm.

Fig.(2) shows for the QW structure the mean square deviation $`\delta \sigma `$ plotted versus the mean location $`\mathrm{X}`$ for the electron (upper panel) and the hole (lower panel) along the X axis (the results are equivalent for the Y axis) for 50 samples, together with the VCA results. The In distribution is at random ($`\mathrm{x}=0.17`$). It is obvious that $`\mathrm{X}`$ and $`\delta \sigma `$ for holes differ strongly from the VCA predictions while the electrons with their lighter mass are less readily confined by In rich fluctuations and thus display $`\delta \sigma (\mathrm{X})`$ which are closer from the VCA results than the holes. Fig.(3) shows the electron - hole optical density of states $`D_{eh}(ϵ)`$ versus the energy $`ϵ`$ averaged over $`N=50`$ samples. $`D_{eh}`$ is proportional to the light absorption coefficient if we neglect electron-hole interaction, which is justified here due to the strong e-h separation by the electric field (770 meV compared to the Coulomb interaction energy of the order of 44 meV). The VCA results in a rounded (by $`\mathrm{\Gamma }`$) staircase - like $`D_{eh}`$ for an infinite QW while the large box produces VCA peaks whose position follows eq. (3). We see that the random alloy results in a broadening of these peaks which is quite important. Note for instance the significant ($`D_{eh}>0.05`$) band tail which develops down to $`15`$ meV below the nominal edge of the ground peak. This bandtail corresponds to increasingly localized states around In clusters, in particular for the holes. The high energy decrease of the optical density of states is unphysical. It reflects the high energy cut-off of the electron and hole eigenvalues. Notice that the photoluminescence line is found experimentally at $``$2.6 eV in $`\mathrm{Ga}_{0.83}\mathrm{In}_{0.17}\mathrm{N}/\mathrm{GaN}`$ Morel1  with which our results are coherent.

Fig.(4) shows the calculated electron density of states for $`\mathrm{Ga}_{0.83}\mathrm{In}_{0.17}\mathrm{N}/\mathrm{GaN}`$ pyramid (100 samples); the optical density of states displays similar results although the hole states can be unbound in the pyramid. The VCA results shows a twofold degenerate excited state located some 90 meV above the ground state. It lies some 80 meV below a narrow triplet. The random alloy case shows both a redshift and a considerable broadening. Note in particular the tendency towards the closing of the gaps between the remnants of the VCA peaks. This effect is more pronounced when In segregation is introduced in the calculations. In addition, increasing the segregation produces densities of states which bear less and less resemblance with the VCA result. In contrast (fig.(5)) the $`\mathrm{Ga}_{0.5}\mathrm{In}_{0.5}\mathrm{As}/\mathrm{GaAs}`$ pyramids (50 samples) display densities of states which are very close to the VCA results. Even though the 50% alloy should display the maximum of disorder, we see that the calculated DOS features are very well separated from each others, and much narrower than the ones found in the 17% nitride alloy. Note in particular that the bandtail extends only a few meV’s below the VCA peaks. Thus, if one knows the shape of a $`\mathrm{Ga}_{0.5}\mathrm{In}_{0.5}\mathrm{As}`$ pyramid, fig.(5) shows that it does make sense to attempt to identify the peaks, to fit their energy difference etc. On the other hand, the huge deviations from the VCA results in the 17% nitride dot demonstrate that any peak attribution in this system is rather elusive.

###### Acknowledgements.

We would like to thank Pr. Runge and Pr. Lefèbvre for fruitful discussions and insights. The LPA-ENS is Laboratoire associé au CNRS et aux Universités Paris 6 et Paris 7 .


text:
# A stochastic cellular automaton model for traffic flow with multiple metastable states

## 1 Introduction

Traffic problems have been attracting not only engineers but also physicists . Especially it has been widely accepted that the phase transition from free to congested traffic flow can be understood using methods from statistical physics . In order to study the transition in detail, we need a realistic model of traffic flow which should be minimal to clarify the underlying mechanisms. In recent years cellular automata (CA)  have been used extensively to study traffic flow in this context. Due to their simplicity, CA models have also been applied by engineers, e.g. for the simulation of complex traffic systems with junctions and traffic signals .

Many traffic CA models have been proposed so far , and among these CA, the deterministic rule-184 CA model (R184), which is one of the elementary CA classified by Wolfram , is the prototype of all traffic CA models. R184 is known to represent the minimum movement of vehicles in one lane and shows a simple phase transition from free to congested state of traffic flow. In a previous paper , using the ultra-discrete method , the Burgers CA (BCA) has been derived from the Burgers equation

$$v_t=2vv_x+v_{xx},$$
(1)
which was interpreted as a macroscopic traffic model . The BCA is written using the minimum function $`\mathrm{min}`$ by

$`U_j^{t+1}=U_j^t`$ $`+`$ $`\mathrm{min}\{U_{j1}^t,LU_j^t\}\mathrm{min}\{U_j^t,LU_{j+1}^t\},`$ (2)
where $`U_j^t`$ denotes the number of vehicles at the site $`j`$ and time $`t`$. If we put the restriction $`L=1`$, it can be easily shown that the BCA is equivalent to R184. Thus we have clarified the connection between the Burgers equation and R184, which offers better understanding of the relation between macroscopic and microscopic models.

The BCA given above is considered as the Euler representation of traffic flow. As in hydrodynamics there is an another representation, called Lagrange representation , which is specifically used for car-following models. The Lagrange version of the BCA is given by

$$x_i^{t+1}=x_i^t+\mathrm{min}\{V_{max},x_{i+S}^tx_i^tS\},$$
(3)
where $`V_{max}=S=L`$ and $`x_i^t`$ is the position of $`i`$-th car at time $`t`$. Note that in (3) $`S`$ corresponds a “perspective” or anticipation parameter  which represents the number of cars that a driver sees in front, and $`V_{max}`$ is the maximum velocity of cars. (3) is derived from the BCA mathematically by using an Euler-Lagrange (EL) transformation  which is a discrete version of the well-known EL transformation in hydrodynamics.

In this paper, we will develop the BCA (3) to a more realistic model by introducing slow-to-start (s2s) effects  and a driver’s perspective $`S`$. Moreover, a stochastic generalization is also considered by combining it with the Nagel-Schreckenberg (NS) model .

## 2 Traffic models in Lagrange form

First, let us extend (3) to the case $`V_{max}S`$ and combine it with the s2s model. The s2s model  is written in Lagrange form as

$`x_i^{t+1}`$ $`=`$ $`x_i^t+\mathrm{min}\{1,x_{i+1}^tx_i^t1,x_{i+1}^{t1}x_i^{t1}1\}.`$ (4)
Note that the inertia effect of cars is taken into account in this model. Comparing (4) and (3), we see that, in the s2s model, the velocity of a car depends not only on the present headway $`d_i^t=x_{i+1}^tx_i^t1`$, but also on the past headway $`d_i^{t1}=x_{i+1}^{t1}x_i^{t1}1`$. This rule has only a nontrivial effect if $`d_i^{t1}=0`$ and $`d_i^t=1`$, i.e. if the leading car has started to move in the previous time step. In this case the following car is not allowed to move immediately (s2s).

Before combining (3) and (4), it is worth pointing out that we can choose the perspective parameter as $`S=2`$ in the model according to observed data. We define the size of a cell as 7.5 m and $`V_{max}=5`$ in our model according to the NS model. Since $`V_{max}`$ corresponds to about 100 km/hour in reality, then one time step in the CA model becomes 1.3 s. Moreover, the gradient of the free line and jamming line in the fundamental diagram, which is the dependence of the traffic flow $`Q`$ on density $`\rho `$, is known to be about 100 km/hour and $`15`$ km/hour  according to many observed data (see Fig. 1) . These values correspond to the typical free velocity and the jam velocity, respectively.

Thus, considering the fact that the positive and negative gradient of each line is given by $`V_{max}=5`$ and $`S/2`$, respectively, in the CA model , we should choose $`S=1.5`$ in the CA model. Since only integer numbers for $`S`$ are allowed in this model, we will simply choose $`S=2`$ for studying the effect of the perspective of drivers. It is noted that other possibilities, such as velocity-dependent $`S`$ or stochastic choice of $`S`$, are also possible.

Now by combining (3) and (4) we propose a new Lagrange model with $`S=2`$ which is defined by the rules listed below: Let $`v_i^{(0)}`$ be the velocity of the $`i`$-th car at a time $`t`$. The update procedure from $`t`$ to $`t+1`$ is divided into five stages:

1. Accerelation

$$v_i^{(1)}=\mathrm{min}\{V_{max},v_i^{(0)}+1\}.$$
(5)
2. Slow-to-accelerate effect

$$v_i^{(2)}=\mathrm{min}\{v_i^{(1)},x_{i+2}^{t1}x_i^{t1}2\}.$$
(6)
3. Deceleration due to other vehicles

$$v_i^{(3)}=\mathrm{min}\{v_i^{(2)},x_{i+2}^tx_i^t2\}.$$
(7)
4. Avoidance of collision

$`v_i^{(4)}`$ $`=`$ $`\mathrm{min}\{v_i^{(3)},x_{i+1}^tx_i^t1+v_{i+1}^{(3)}\}`$ (8)
5. Vehicle movement

$$x_i^{t+1}=x_i^t+v_i^{(4)}.$$
(9)

The velocity $`v_i^{(4)}`$ is used as $`v_i^{(0)}`$ in the next time step. (8) is the condition that the $`i`$-th car does not overtake its preceding $`(i+1)`$-th car, including anticipation. The accerelation (5) is the same as in the NS model, which is needed for a mild accerelating behaviour of cars. In the step 2, we call (6) as “slow-to-accelerate” instead of s2s. This is because this rule affects not only the behaviour of standing cars but also that of moving cars, which is considered to be a generalization of s2s rule.

It is not difficult to write down the new model in a single equation for general $`S`$. The result is

$`x_i^{t+1}`$ $`=`$ $`x_i^t+\mathrm{min}\{V_i^t,\underset{k=1,\mathrm{},S1}{\mathrm{min}}(x_{i+k}^tx_i^tk+V_{i+k}^t)\},`$ (10)
where the last term represents the collision-free condition explained in Fig. 2, and

$`V_i^t`$ $`=`$ $`\mathrm{min}\{V_{max},x_{i+S}^{t1}x_i^{t1}S,x_{i+S}^tx_i^tS,x_i^tx_i^{t1}+1\}.`$ (11)
The condition that there is no collision between the $`i`$-th and $`i+k`$-th cars ($`k=1,\mathrm{},S1`$) is given by

$$x_{i+k}^tx_i^tk+V_{i+k}^tV_i^t,$$
(12)
for $`S2`$ (if $`S=1`$ then we simply put $`k=1`$), which is identical to the last term in (10). In contrast to the NS model, the velocity of the preceeding car is taken into account in the calculation of the safe velocity in the step 4, i.e. our model also includes anticipation effects.

## 3 Metastable branches and their stability

Next, we investigate the fundamental diagram of this new hybrid model.

In Fig. 3, we observe a complex phase transition from a free to congested state near the critical density 0.2 $``$ 0.4. There are many metastable branches in the diagram, similar to our previous models in Euler form  or in other models with anticipation . We also point out that there is a wide scattering area near the critical density in the observed data (Fig. 1) which may be related to these metastable branches. As we will discuss later, these branches may account for some aspects of the scattering area observed empirically.

First, we discuss properties of the state in the metastable branches. In all cases it consists of pairs of vehicles that move coherently with vanishing headway (see Fig. 4). Cars are represented by black squares, and the direction of the road is horizontal right and time axis is vertical down. The corresponding velocity distributions are also given in Fig. 5. We see that there are stopping cars which velocity are zero only in the case of the lowest branch given in the state in Fig. 4 (e).

Next let us calculate the flow-density relation for each branch. In the metastable branches we find phase separation into a free-flow and a jamming region. In the former, pairs move with velocity $`v_f`$ and a headway of $`d_f`$ empty cells between consecutive pairs. In the jammed region, the velocity of the pairs is $`v_j`$ and the headway $`d_j`$. $`N_j`$ and $`N_f`$ are the numbers of cars in the jamming cluster and the free uniform flow, respectively. We assume $`N_f`$ and $`N_j`$ to be even so that there are $`N_f/2`$ and $`N_j/2`$ pairs, respectively. Then the total number of cars $`N`$ is given by $`N=N_j+N_f`$ and the total length of the system becomes $`l=(d_j+2)N_j/2+(d_f+2)N_f/2`$. Since the average velocity is $`\overline{v}=(N_fv_f+N_jv_j)/N`$ and density and flow of the system are given by $`\rho =N/l`$ and $`Q=\rho \overline{v}`$, we obtain the flow-density relation as

$$Q=2\frac{v_fv_j}{d_fd_j}+\left(v_j(d_j+2)\frac{v_fv_j}{d_fd_j}\right)\rho .$$
(13)
From the stationary states in Fig. 4 we have

$`(a):(v_f,v_j,d_f,d_j)=(5,4,6,4)`$
$`(b):(v_f,v_j,d_f,d_j)=(5,3,7,3)`$
$`(c):(v_f,v_j,d_f,d_j)=(5,2,8,2)`$
$`(d):(v_f,v_j,d_f,d_j)=(5,1,9,1)`$
$`(e):(v_f,v_j,d_f,d_j)=(5,0,10,0)`$ (14)
Therefore the resulting equations for each branch are

$$Q=1+c\rho $$
(15)
where $`c=1,1/2,0,1/2,1`$, which correspond to the branches $`A,B,C,D`$ and $`E`$ in Fig. 6, respectively. End points of the branches are found to be given by $`\rho _1=2/(d_f+2)`$ and $`\rho _2=2/(d_j+2)`$, where $`\rho _1`$ is the point that the metastable branches intersect the free flow branch, and $`\rho _2`$ is the maximal possible density in the metastable branches. Note that all $`\rho _2`$ lie on the line $`Q=2\rho +2`$, which is indicated as the broken line in Fig. 6.

Next let us now study the stability of each metastable branch. We mainly consider the density $`\rho =2/7`$ and, in particular, we will focus on the uniform flow represented by $`\mathrm{}11000001100000\mathrm{}`$, which shows the highest flow given in the point $`P_0`$ in Fig. 6.

Spatio-temporal patterns due to various kinds of perturbations are already seen in Fig. 4. Perturbation in this case means that some cars are shifted backwards at the initial configuration. The initial conditions for Fig. 4(a)-(e) are given as follows:

1. Very weak perturbation (one car is shifted one site backwards)$`\mathrm{}11000010100000\mathrm{}`$
2. Weak perturbation (one car is shifted two sites backwards)$`\mathrm{}11000100100000\mathrm{}`$
3. Moderate perturbation (one car is shifted three sites backwards)$`\mathrm{}11001000100000\mathrm{}`$
4. Strong perturbation (one car is shifted five sites backwards)$`\mathrm{}11100000100000\mathrm{}`$
5. Strongest perturbation (three cars are shifted backwards) $`\mathrm{}11111000000000\mathrm{}`$

The stationary state of (a), (b), (c), (d) and (e) are given by the points $`P_A`$, $`P_B`$, $`P_C`$, $`P_D`$ and $`P_E`$ in Fig. 6, respectively. That is, if the system in $`P_0`$ is perturbated, then the flow easily goes down to a lower branch in the course of time depending on the magnitude of the perturbation. Since the density does not change due to the perturbation, we obtain $`P_A:(2/7,9/7)`$, $`P_B:(2/7,8/7)`$, $`P_C:(2/7,1)`$, $`P_D:(2/7,6/7)`$ and $`P_E:(2/7,5/7)`$ by substituting $`\rho =2/7`$ into eq. (15).

We see a jamming cluster propagating backwards in the cases of (d) and (e) in Fig. 4. In other cases the jamming cluster propagates forward ((a) and (b)) or does not move (c). These facts are related to the gradient of the metastable branches which are given by $`c`$ according to eq. (15).

## 4 Stochastic generalization

Finally we will combine the above model with the NS model in order to take into account the randomness of drivers. The NS model is written in Lagrange form as

$`x_i^{t+1}=x_i^t+\mathrm{max}\{0,\mathrm{min}\{V,x_{i+1}^tx_i^t1,x_i^tx_i^{t1}+1\}\eta _i^t\}.`$ (16)
where $`\eta _i^t=1`$ with probability $`p`$ and $`\eta _i^t=0`$ with probability $`1p`$. The last term in the mininum in (16) represents the acceleration of cars. The randomness in this model is considered as a kind of random braking effect, which is known to be responsible for spontaneous jam formation often observed in real traffic . We also consider random accerelation in this model which is not taken into account in the NS model.

Thus a stochastic generalization of the hybrid model in the case of $`S=2`$ is similarly given by the following set of rules:

1. Random accerelation

$$v_i^{(1)}=\mathrm{min}\{V_{max},v_i^{(0)}+\eta _a\}.$$
(17)
where $`\eta _a=1`$ with the probability $`p_a`$ and $`\eta _a=0`$ with $`1p_a`$.
2. Slow-to-accelerate effect

$$v_i^{(2)}=\mathrm{min}\{v_i^{(1)},x_{i+S}^{t1}x_i^{t1}S\}.$$
(18)
3. Deceleration due to other vehicles

$$v_i^{(3)}=\mathrm{min}\{v_i^{(2)},x_{i+S}^tx_i^tS\}.$$
(19)
4. Random braking

$$v_i^{(4)}=\mathrm{max}\{v_i^{(3)}\eta _b,0\}$$
(20)
where $`\eta _b=1`$ with the probability $`p_b`$ and $`\eta _b=0`$ with $`1p_b`$.
5. Avoidance of collision

$`v_i^{(n+1)}`$ $`=`$ $`\mathrm{min}\{v_i^{(n)},x_{i+1}^tx_i^t1+v_{i+1}^{(n)}\}`$ (21)
with $`n4`$, which is an iterative equation that has to be applied until $`v`$ converges to $`v_i^{(n+1)}=v_i^{(n)}(v_i)`$.
6. Vehicle movement

$$x_i^{t+1}=x_i^t+v_i.$$
(22)

Again the velocity $`v_i`$ is used as $`v_i^{(0)}`$ in the next time step. Step 5 must be applied to each car iteratively until its velocity does not change any more, which ensures that this model is free from collisions. This is the difference between the deterministic and stochastic case. In the deterministic model it is sufficient to apply the avoidance of collision stage only once in each update, while in the stochastic case generically it has to be applied a few times in order to avoid collisions between successive cars.

The fundamental diagrams of this stochastic model for some values of $`p_a`$ and $`p_b`$ are given in Fig. 7. The randomization effect can be considered as a sort of perturbation to the deterministic model. Hence some unstable branches seen in the deterministic case disappear in the stochastic case, especially if we consider the random braking effect as seen in Fig. 7. Random accerelation itself does not significantly destroy the metastable branches. Moreover, from the spatio-temporal pattern it is found that spontaneous jam formation is observed only if we allow random braking. Random accerelation alone is not sufficient to produce spontaneous jamming. We also note that a wider scattering area appears if we introduce both random accerelation and braking.

## 5 Concluding discussions

In this paper we have proposed a new hybrid model of traffic flow of Lagrange type which is a combination of the BCA and the s2s model. Its stochastic extension is also proposed by further incorporating stochastic elements of the NS model and random accerelation. The model shows several metastable branches around the critical density in its fundamental diagram. The upper branches are unstable and will decrease its flow under perturbations. It is shown that the magnitude of a perturbation determines the final value of flow in the stationary state. Moreover, introduction of stochasticity in the model makes the metastable branches dilute and hence produces a wide scattering area in the fundamental diagram. We would like to point out that this metastable region around the phase transition density is similar to so-called synchronized flow proposed by . Our investigation shows that one possible origin of such a region is the occurance of many intermediate congested states near the critical density. If some of them are unstable due to perturbation or randomness, then a dense scattering area near the critical density is formed around the metastable branches. This is in some sense in between the two cases of a fundamental diagram based approach (with unique flow-density relation) and the so-called 3-phase model of  which exhibits a full two-dimensional region of allowed states even in the deterministic limit.

## Acknowledgment

This work is supported in part by a Grant-in-Aid from the Japan Ministry of Education, Science and Culture.

## References