text:
# AGAPEROS: Searching for variable stars in the LMC BarThis work is based on data collected by the EROS and DENIS collaborations.

## 1 Introduction

The late stages of stellar evolution are characterized by regular and irregular light variability, a well-known signature of the stellar pulsation of Asymptotic Giant Branch stars (hereafter AGB). These light changes allow to identify AGB stars over large distances and to derive the pulsation characteristics (periodicity, etc.), which are key parameters for understanding the fundamental properties of the highly extended atmospheres of these stars. The pulsational properties have a strong impact on the structure of AGB stars. Pulsation, as the driving mechanism for the stellar winds, plays a key role for the high mass loss rates reached during the AGB phase.

A new era in the study of variable red giants started, when microlensing surveys produced a large amount of light curves of these stars, especially for objects in the LMC. The pioneering work by Wood (Wood2000 (2000)), using data from the MACHO survey, showed that the red giant variables form four roughly parallel sequences in a period-magnitude diagram. Three of these sequences could be associated with fundamental, first and second overtone pulsation. The explanation of the fourth sequence is not clear yet (Wood Wood2000 (2000), Hinkle et al. Hinkle2002 (2002)). Recently, Cioni et al. (Cioni2001 (2001)) presented a survey of variable red giants in the LMC based on data from the EROS-2 microlensing survey (Lasserre et al. Lasserre2000 (2000)). The work of Cioni et al. focused on the logP$``$K relation confirming three of the relations found by Wood, and they discussed the behaviour of different groups of variables in near infrared colour-magnitude diagrams.

The present paper relies on the variability information contained in the AGAPEROS variable star catalogue (Melchior et al. Melchior2000 (2000)). Here, we extend and analyse the corresponding light curves, relying on the EROS-1 microlensing survey data set (Ansari et al. Ansari95 (1995), Aubourg et al. Aubourg95 (1995)). These data have been obtained between December 1991 and April 1994. The second half of the data set therefore overlaps with the MACHO survey. We combine the EROS data with IJ$`\mathrm{K}_\mathrm{S}`$ photometry of the DENIS survey (Epchtein et al. Epchtein97 (1997)), with an approach similar to the work of Cioni et al. Cioni2001 (2001). The intention of our work is to discuss the light change of red variables on a large and homogeneous sample and to compare the results for LMC red variables with the corresponding objects in the Galactic Disk and Bulge.

## 2 Data

We have studied the AGAPEROS catalogue of 584 variable stars detected over a 0.25 deg2 field in the LMC Bar (Melchior et al. 2000, hereafter referred to as Paper I). These stars have been selected on the basis of their variability on a 120-days window with a bias towards long-timescale variations ($`>`$ few days), which are not necessarily periodic. The original data set has been taken at ESO by the EROS-1 collaboration (Arnaud et al., 1994a,b, Renault et al. 1997) using a 40cm telescope, equipped with a wide field camera composed of 16 CCD chips, each of 400$`\times `$579 pixels of 1.21 arcsec (Arnaud et al. 1994b ). We use 9 chips, and study light curves in the red ($`\overline{\lambda }=670`$nm) filter. Table 1 gives a short description of the dataset. See Melchior et al. (Melchior98 (1998),Melchior99 (1999)) for a more detailed description of the data treatment.

The study of the position of these stars in the colour-magnitude diagram showed that this catalogue is dominated by a population of Long Timescale & Long Period variables, while a few “bluer” variables have also been detected. A cross-correlation with various existing catalogues showed that about 90$`\%`$ of those variable objects were undetected before. We extended the corresponding light curves to the whole EROS-1 database of the LMC (900-days). We improved the photometry of the corresponding light curves with image subtraction using the ISIS2.1 algorithm of Alard (Alard2000 (2000)). We follow the same definition of the magnitude system as in Paper I, but we rely here on image subtraction photometry.

The published DENIS catalogue (I, J and $`\mathrm{K}_\mathrm{S}`$) for the LMC (Cioni et al. Cioni2000 (2000)) has been used to make a cross-identification between the AGAPEROS variables and the DENIS magnitudes. A search radius of 3″was chosen to avoid misidentification. Out of the 584 variables 468 were detected by DENIS ($``$ 80%). Note that the positional accuracy of these variables is about 1″, as discussed in Paper I.

## 3 Classification of the variables

Classically, three types of variable red giants have been defined (General Catalogue of Variable Stars, GCVS, Kholopov et al. 1985- (88)): Mira-type variables show periodic large amplitude variations with time scales typically of the order of 200 to 500 days. Semiregular variables (SRVs) show a less regular behaviour and a smaller amplitude. A typical time scale of the variation can be found, but the light curve shows phases of irregularity as well. The GCVS has introduced a limiting amplitude of 2.5 mag to separate Miras and SRVs. While even within the GCVS this rule has not been strictly applied (see e.g. the SRV W Hya), several investigators used this simple criterion for classification (e.g. Alard et al. 2001, Cioni et al. Cioni2001 (2001)). The artificial nature of this division has been criticized already e.g. by Kerschbaum (Kerschbaum93 (1993)). The third group of variables are the irregular variables. It is still not clear if such stars really exist or if these objects are simply not observed well enough to detect the same amount of periodicity as in the SRVs (e.g. Lebzelter et al. LKH95 (1995)).

In this work, we used a different approach to classify the light curves of the red giants in our sample. The classic classification system depends on whether a more-or-less constant period can be found and also depends on an arbitrary amplitude limit. Here we adopt a new scheme which is based on how well the light curve can be described by one or two periods only and where amplitude plays no role. In this way we are able to separate the two effects amplitude and regularity.

We based the classification on the regularity of the light curve on a visual comparison of the light change with a combination of up to three sine curves. To derive the periods, a Fourier analysis of the light curves (based on the program Period98 by Sperl (Sperl98 (1998))) has been applied. Semiregular and irregular light changes result in a large number of peaks of similar strength in the Fourier spectrum (Lebzelter Lebzelter99 (1999), see below). Therefore the periods used for the fit have been selected from peaks in the periodogram by visual inspection. The amplitudes of the peaks were the starting point for the selection of the periods. Naturally, this selection is influenced by aliases. Fig.2 shows the typical spectral window of our data that has been used to identify spurious peaks. The selected periods were always cross checked by a visual comparison wiht the light curve. For unclear cases a second Fourier analysis was made with the primary period subtracted.

As a first approach to this large amount of light curve data we made no attempt to fit every detail of the light curves but identified the major period(s) to roughly resemble the overall light change. A more detailed fitting, as it was done by e.g. Kerschbaum et al. (Kerschbaum2001 (2001)) for a small number of SRVs in the solar neighborhood, is planned. Finally, we stress that the total available baseline of the data set did not allow to derive periodicities on time scales longer than 900 days. The classification is based on three years of observation and represents the behaviour of each object over the 900-days window. Stars classified as irregular may show some periodicity on a longer time scale or during a different time interval. The amplitudes were estimated visually from the lightcurve.

We classified the light curves on the regularity and type of their light change into four groups:

* Regular: A constant cycle length is observed for the available measurements. Some of these objects show some amplitude variations or bumps in their light curves. This group includes also stars which show a second period, if the two periods allow a very good fit of the light curve. While this group will include almost all objects that classically would have been classified as Miras, possible small amplitude regular variables will be found in this class as well. We therefore did not use the name ”Miras” for this group. The only regular pulsators we did not include here were cepheids and obvious binaries, which have been classified as Other.
* Semiregular: The cycle length is variable, but some kind of periodicity is visible. The stars show up to three strong peaks in the Fourier spectrum. In some cases the light change can be fitted rather well with three or four periods.
* Irregular: These stars do not show any significant periodicity in their light change. Their light change occurs on time scales typical for AGB stars (i.e. a few 10 to a few 100 days). Typically, the Fourier spectrum shows a large number of peaks with similar strength.
* Other: This group includes stars with a large fraction of bad data points, stars that turned out to be constant (misidentification due to some erroneous data points), and stars with a luminosity variation atypical for long period variables (including a number of binaries). A large fraction of these objects have also no DENIS data.

Fig. 1 shows a sample of regular, semiregular and irregular light curves. Note that our classification does not take into account the amplitude of the variation as in the GCVS classification. The examples were selected to represent the different expressions of variability found in the three groups. Among the regular variables, we included examples of amplitude variations (top left in Fig. 1), variables with two periods (top middle) and classical Mira variables (top right).

Naturally, this classification remains somewhat subjective. However, we made an attempt to check the homogeneity of our classification by using the Fourier spectra of the light curves. In Fig. 3, we plot the amplitude ratio of the strongest and the second strongest peak against the ratio of the second and the fifth strongest peak. The advantage of our sample is that all light curves have a similar sampling and therefore a similar spectral window. Examples for Fourier spectra and a spectral window are given in Fig. 4 and 2, respectively. A small number of stars has been excluded from this plot as their time coverage is not as good as for the majority of the sample.

As mentioned above, a semiregular or irregular light curve typically results in a number of peaks of similar strength in the Fourier spectrum. Stars classified as regular have only one or two strong peaks in their Fourier spectrum. They should therefore be found on the right-hand side and the top side of Fig. 3. Stars with a single period are on the right, stars with a second period in the upper left region of the plot. Note that there is no correction for aliases in this approach. From the spectral window (Fig. 2) one would expect to find stars with a single period at a ratio a1/a2 of about 2.

For the regular variables the fifth strongest peak is typically already at the noise level and was used as a reference point. Note that we did not use more than three periods for each object in the following analysis. On the other hand, irregular variables should be found in the lower left corner of the plot. Semiregular stars are expected in between. Fig. 3 shows this classification indicated by different symbols.

We observe that our classification criteria is coherent within our sample. However, for an individual object, Fig. 3 is not usable for classification as the borders between the three classes are not well defined.

For each star classified as regular, semiregular or irregular a typical amplitude of the light variation was determined. In the case of semiregular and irregular variables the light amplitude can change dramatically. In these cases, we used a mean value of the variation. As no standard Johnson filters have been used a direct comparison of the amplitude values found here and those given in the GCVS or the MACHO catalogue is not possible.

## 4 Comparison with results from the MACHO survey

We searched the MACHO Variable Star Catalogue (http://wwwmacho.mcmaster.ca/Data/MachoData.html) for variables within the fields covered by our sample and classified as LPV.WoodA, LPV.WoodB, LPV.WoodC and LPV.WoodD, respectively. The catalogue released on the web includes only a subsample of all variables found in the MACHO survey. It was therefore not surprising that we did not find all stars of our sample in the released MACHO catalogue, which is probably not complete in the area we are concerned with. In total, we found 36 MACHO LPVs that are within the fields we investigated. 25 of them had counterparts in our sample111In the course of this comparison, we found that the star 78.5861.10/80.6466.5194 has two entries in the web based MACHO catalogue. We assume that the remaining 11 stars are located on defects or borders of the CCD chips, but we did not investigate these objects further.

We applied the same analysis to the red MACHO light curves of these 25 stars. Results are listed in Table 2. The same classification was reached for 21 stars in our comparison, in two further cases either the MACHO or the AGAPEROS light curve was not of sufficient quality for the analysis. Interestingly, the two remaining objects with different classifications both show a higher degree of regularity in the MACHO data. In these cases our dataset was obviously not covering enough light cycles to reveal the regularity. We can estimate from this result that for less than 10% of the AGAPEROS light curves the regularity was not detected correctly.

22 of the LPVs in common with the MACHO catalogue agree in the main period within a few percent. The values we derived for the MACHO light curves are in good agreement with the results listed in the MACHO catalogue. However, the catalogue gives only one period for each object, so in cases we found multiple periods in our sample stars only one of them could be compared. Several of the secondary periods found in the AGAPEROS stars also agree quite well with values from the MACHO data. Very long periods could not be detected with the shorter time series of AGAPEROS data. The differences in some secondary periods illustrates the difficulty to derive a unique fit of these light curves with more than one period (see also Kerschbaum et al. Kerschbaum2001 (2001)).

The good agreement in the classification between the AGAPEROS and the MACHO data, computed over independent data sets (1991-1994 for AGAPEROS, 1992-2000 for MACHO), strengthens the validity of our approach.

## 5 General characteristics of the sample variables

Semiregular variables clearly dominate our sample of late type giants. 583 light curves have been analysed in total. 112 of them have been classified as other. Among the remaining 471 objects we classified 18 % as regular, 67 % as semiregular and 15 % as irregular. Due to different classification criteria a comparison of this result with other investigations is difficult. Cioni et al. (Cioni2001 (2001)) found that 65 % of the AGB stars are variable. As in our sample semiregular variables (in their case small amplitude red variables) are the dominant group of objects, only 12 % of their stars have been classified as Miras.

If more than one period is detected, we have chosen a ’primary’ period based on its amplitude. However, this decision is somewhat subjective when the amplitudes are similar. “Second periods” are present in at least 54 % of the stars classified as semiregular and about 18 % of the regular variables. For the semiregular variables the ratio between the short and long period is between 1 and 2 for about 36 % of the objects (Fig. 5), the other semiregulars have period ratios between 2 and 15. This result is similar to what has been found in previous investigations of SRVs in our Galaxy (e.g. Kiss & Szatmary Kiss2000 (2000)). Among the regular variables with two periods about one third of the stars show a period ratio below 2.

Comparing amplitude with period shows a maximum amplitude for periods between 300 and 400 days, given our 900-day window. No large amplitude variables with periods below 100 days have been found. Details are illustrated in Fig. 6 where the mean amplitude has been calculated for each period bin. Alard et al. (Alard2001 (2001)) found a similar increase of mean amplitude from small to long periods in a sample of AGB variables in the Galactic Bulge. However, their sample shows a maximum amplitude at periods around 250 days, as they excluded all miras from their analysis. Fig. 7 summarizes the amplitude distributions within each class of variables. All three classes are dominated by small amplitude objects. Large amplitude stars occur exclusively among regular variables, all irregular objects are small amplitude stars.

We will now mainly concentrate on regular and semiregular variables. Fig. 8 shows the period distribution of the AGAPEROS sample. Both groups of objects (regular and semiregular) show a maximum at the shortest periods: Therefore, the regular variables cannot be simply related to the class of Mira variables. No Galactic Mira with a period below 100 days is known. Furthermore, Miras typically show large amplitude variability while the short period regular stars in our sample all have small amplitudes. The ”classical” Miras probably form the second maximum in period distribution of the regular variables around 350 days, similar to the Galactic Miras. We conclude that the group classically known as semiregular variables may contain a substantial number of periodic variables. On the other hand it is known from long time series of Galactic semiregular variables that these stars can show phases of periodic behavior, so that these stars may exhibit semiregular behavior on longer time scales.

## 6 Near-Infared data

The near-infrared data from the DENIS survey allow to characterize the variables of our sample in more detail concerning their luminosity and chemical composition. Fig. 9 shows the $`\mathrm{K}_\mathrm{S}/(\mathrm{J}\mathrm{K}_\mathrm{S})`$ diagram for the AGAPEROS variables. One can clearly see that the majority of the sources are located above the tip of the Red Giant Branch (hereafter RGB-tip) which is for the LMC about 12.0 mag in $`\mathrm{K}_\mathrm{S}`$ (Cioni et al. 2000a ). We find regular and semiregular variables which are below the RGB-tip. These objects have rather short periods ($`<`$ 100 days) and could be AGB stars in the early evolutionary phase (early-AGB phase) or variable stars on the red giant branch.

Carbon-rich objects are characterized by their red (J–K) and (I–J) colour compared to the oxygen-rich sequence (see Cioni et al. Cioni99 (1999)). However, as noted by Loup et al. (Loup2002 (2002)) the colour-colour diagram is just a statistical tool to distinguish between oxygen-rich and carbon-rich objects. Fig. 10 shows the $`(\mathrm{I}\mathrm{J})_0`$ vs $`(\mathrm{J}\mathrm{K})_0`$ diagram. Obviously, the ratio of regular to semiregular variables is smaller for the oxygen-rich stars than for the carbon-rich objects. This suggests that the majority of the semiregular variables are less massive than Miras which prevents them from becoming carbon stars.

We do not find any significant difference in colours or luminosites between SRVs with one single period and SRVs with multiple periods.

### 6.1 Colour-Period diagrams

For the LMC bar, it is obvious from Fig. 11 that the AGAPEROS variables follow a tight logP vs I–J relation. It is important to emphasize that the I magnitudes of DENIS correspond to a single epoch measurement and thus the logP vs. I–J diagram is affected by the scatter due to the amplitude variation of each source.

Relying on MACHO data in the Galactic Bulge, Schultheis & Glass (Schultheis2001 (2001)) demonstrated that semiregular variables in the Galactic Bulge show a noticeable scatter in I–J (3–4 mag) along the logP vs I–J relation. The most significant difference between the Galactic Bulge and the LMC is the smaller range in I–J for the LMC ($``$ 2 mag) than for the Bulge ($``$ 4 mag). (see Fig. 11)

The I band for M stars is mostly affected by the strong TiO and VO molecular absorption (Turnshek et al. Turnshek85 (1985), Lancon & Wood Lancon2000 (2000)). Schultheis et al. (Schultheis99 (1999)) showed that lower metallicity is correlated with weaker TiO band intensities, corresponding to bluer I–J colours. The large scatter and the wide I–J range in the Galactic Bulge sample compared to the LMC might be explained by the wide spread in metallicity compared to the Magellanic Clouds. However, the difference in the I–J range between the Galactic Bulge ($`1<(\mathrm{I}\mathrm{J})_0<5`$) and the LMC ($`1<(\mathrm{I}\mathrm{J})_0<3`$) seems rather large. A more detailed quantitative analysis, using realistic model atmospheres of AGB stars (including metallic lines), is necessary to fully understand this systematic difference in the I–J colour between the Galactic Bulge and the LMC.

Fig. 12 displays the J–K colours of the AGAPEROS variables as a function of their period. The majority of the SRVs appear to follow a different period-colour relation with a slope flatter than the regular variables. For comparison, we indicated in Fig 12 the averaged colours of oxygen-rich Miras for the SgrI field (Glass et al. Glass95 (1995)). The majority of our long-period Miras ($`\mathrm{log}\mathrm{P}>250^\mathrm{d}`$) follow the location of the oxygen-rich Miras in SgrI. The carbon rich objects (J-K $`>`$ 1.6) seem to form a parallel sequence to the oxygen-rich Miras, while the long-period SRVs ($`\mathrm{P}>300^\mathrm{d}`$) do show clearly another period-colour relation. These stars are located on the sequence D in Wood’s diagram (see Wood et al. Wood99 (1999) and discussion below) and are SRVs with multiple periods. A few long-period Miras also follow this sequence. However, the scatter in this diagram increases for $`\mathrm{log}\mathrm{P}>2.3`$ due to the contribution of the circumstellar dust shell arising from mass loss. Schultheis et al. (Schultheis99 (1999)) and Schultheis & Glass (Schultheis2001 (2001)) obtain similar results for semiregular variables in the Galactic Bulge (see their Fig. 8).

### 6.2 $`\mathrm{K}_\mathrm{S}`$ vs logP diagram

In the Large Magellanic Cloud, the Miras and the SRVs seem to form distinct parallel sequences C,B,A which have been identified by Wood (Wood2000 (2000)) as pulsators in the fundamental, first and the next two higher overtones, respectively. Wood et al. (Wood99 (1999)) showed by comparison of observed periods, luminosities and period ratios with theoretical models, that Miras are radial fundamental mode pulsators, while semiregular variables can be pulsating in the 1st, 2nd or 3rd overtone, or even the fundamental mode. The pulsation mode derived by Whitelock & Feast (WF00 (2000)) from diameter measurements of Miras in the Milky Way suggests first overtone pulsation for Miras. However, observations of radial velocity variations of Miras (e.g. Hinkle et al. HHR (1982)) clearly favour fundamental mode pulsation (Bessell et al. BSW96 (1996)).

In Fig. 13, we distinguish between semiregular variables with one period and those having a second or even third pulsational period. The location of our regular variables is consistent with the PL-relation from Feast et al. (Feast82 (1989)) and Wood’s sequence C corresponding to fundamental mode pulsation. However, a few regular variables are also found to be located on sequence B and A (first and second overtones according to Wood (Wood2000 (2000))). The majority of the SRVs follow Wood’s sequence B although the scatter is rather large ($``$ 0.5 mag in $`\mathrm{K}_\mathrm{S}`$ at a given period). The SRVs situated on sequence A show very low amplitudes ($`<`$ 0.5 mag in $`\mathrm{R}_{\mathrm{EROS}}`$) and typically no secondary periods. While Cioni et al. (Cioni2001 (2001)) found no objects on sequence A, we could clearly confirm the existence of this PL-sequence. On sequence B and C, we find both single periodic and multiperiodic objects. The occurrence of single or multiple periodic behaviour does not depend on the luminosity.

Several data points also mark sequence D of Wood (Wood2000 (2000)). The large scatter in this part of the K-logP-diagram is due to the limited time window of our data set. We are therefore able to reproduce all four sequences found in the MACHO data. The PL-relation for SRVs found by Bedding & Zijlstra (Bedding98 (1998)) from local objects could not be confirmed with our data (see below).

## 7 Discussion

### 7.1 Variability on the AGB

According to Fig. 9 most of the variables in our sample are on the AGB. Therefore, we can use our results to discuss the variability during the AGB phase. Our classification system for the type of variability aims to measure the regularity of the light change. Even not taking into account variations in the amplitude of the light change, we show that most stars have light curves that cannot be fitted by the simple combination of one or two excited periods. Regular variations are found with a wide range in period, while semiregular variability typically occurs mainly on time scales below 150 days (see Fig. 8). In Fig. 14, we compare the period distribution of the semiregular variables in our sample with the Milky Way SRVs listed in the GCVS. While in both cases the maximum of the distribution is at short periods, the GCVS distribution shows a significantly larger fraction of stars with periods longer than 150 days. These long periods may have been missed by our rather short time window. One would also expect a bias of the GCVS sample towards large amplitude variables as most of the data used there are based on photographic measurements. Furthermore, the period distribution from the GCVS given in Fig. 8 includes only one (main) period per object, while for the AGAPEROS data we give also secondary periods found for these stars.

Due to the separation of amplitude and regularity in our classification system, we can explore the relation between these two quantities. We find that large amplitude variation occurs almost exclusively among the regular variables (see Fig. 7). However, there exist regular pulsators with small amplitudes. It is therefore not correct to classify all red variables below a certain amplitude limit as semiregular. A division into large and small amplitude variables seems to be more meaningful. Large and small amplitude variables are both found all along the AGB. This is illustrated in Fig. 15 where the $`\mathrm{R}_{\mathrm{EROS}}`$ light amplitude is plotted against the DENIS K band measurement. Towards the tip of the AGB the fraction of regular as well as large amplitude variables increases. Below the RGB-tip, amplitudes become on the average smaller. The occurrence of regular and semiregular as well as small and large amplitude variables on the AGB indicates that AGB stars have to be seen as a highly inhomogeneous group. One reason for this may be a difference in stellar mass as noted above.

Summarizing, large amplitudes are well correlated with regular pulsations, but we find no correlation between large amplitude and stellar luminosity nor between small amplitude variability and semiregularity of the light change. This result is in agreement with Wood et al. (Wood99 (1999)).

### 7.2 PL-relation

In the literature, the observed PL-relation of long-period variables is considered to be the same in different environments such as the LMC, the Galactic Bulge or globular clusters (see e. g. Glass et al. Glass95 (1995), Feast et al. Feast2002 (2002)). It is therefore independent of metallicity, contrary to the predictions of pulsation theory (see e. g. Wood & Sebo Wood96 (1996)). However, previous studies were restricted to Mira variables mainly due to limitation of sensitivity. Thanks to the microlensing surveys such as EROS, MACHO or OGLE, we can study systematically small amplitude variations over a (still rather small) time interval. Wood (Wood2000 (2000)) found for the SRVs in the LMC different PL-relations for different pulsational modes. However, these separations cannot be reproduced in the Galactic Bulge (Schultheis & Glass 2001). In addition, the PL-relation of the solar neighborhood (Bedding & Zijlstra Bedding98 (1998)) looks different. Why is the PL-relation the same for Miras in different galactic environments, but not for SRVs?

Fig. 13 shows the PL-relation for the AGAPEROS sample. On the one hand, the Mira variables, classically defined as long period and large amplitude stars, concentrate along Wood’s sequence C. Regular variables at shorter periods would not have been classified as Miras. On the other hand, the semiregular variables (both according to the classical and to our definition), are spread all over the K-log P-plane. Making one fit with all semiregular stars would not result in a K vs. log P relation. In the solar neighborhood, Bedding & Zijlstra (Bedding98 (1998)) note that the SRVs are actually found on two sequences: the first one corresponds to the LMC Mira PL-relation (Wood’s sequence C); the second one is located close to a PL-relation derived from Galactic globular cluster LPVs shifted 0.8 mag from the Whitelock globular cluster sequence (Whitelock Whitelock86 (1986)), as shown in Fig. 13. The Bedding & Zijlstra sequence, defined for SRVs, obviously mixes objects from Wood’s sequence B and C, as shown in Fig. 13. The increase towards longer periods is consistent with the larger fraction of long period SRVs in the GCVS (Fig. 14) assuming that the detection of long periodic small amplitude variations is biased towards bright objects. Therefore, three PL-sequences seem to be more appropriate for semiregular variables. Multiperiodic stars are found on all three sequences A, B and C (see Fig. 13). Sequence D is almost exclusively occupied by stars with two periods in agreement with the suggestion from Wood (Wood2000 (2000)) that these long periodic variations are either due to binarity or a pulsation mode resulting from an interaction of pulsation and convection. However, there are also a few regular pulsating variables on this sequence with only one period. These stars would be definitely worth further investigation.

Schultheis & Glass (Schultheis2001 (2001)) showed that the interpretation of the PL-relation of Bulge SRVs is rather complex due to the depth of the Bulge ($``$ $`\pm 0.35^{\mathrm{mag}}`$, see Glass et al. Glass95 (1995)) and the variable interstellar extinction. There is no clear separation of the four sequences. We also showed that the LMC variables are much more homogeneous in their metallicity than the Bulge AGB stars (Fig. 11). This would explain part of the scatter in the K vs. log P plot for the Bulge.

### 7.3 Number densities

The number of semiregular variables in comparison to the regular variables is about a factor of 3. If we use the selection criterion of Cioni et al. (Cioni2001 (2001)), i.e. all stars with $`\mathrm{R}_{\mathrm{EROS}}`$ amplitudes smaller than 0.9 mag are SRVs, we end up with a ratio of almost 37 between SRVs and Miras in our sample. This value is much higher than what was found by Cioni et al. ($``$5), so we assume that our sample is more complete at smaller amplitudes. In the Galactic Bulge, Alard et al. (Alard2000 (2000)) found that the proportion of SRVs with respect to Miras is about a factor of 20. Most recently, Derue et al. (Derue02 (2002)) found a similarly large ratio between semiregulars and miras in the Galactic spiral arms. However, this ratio is of course very sensitive to the classification of SRVs (see above). For the Galactic disk, Kerschbaum & Hron (KH92 (1992)) found equal number densities for Miras and semiregular variables. However, they note that their sample of semiregular variables is probably not complete due to the difficulties in detecting small amplitude variables.

Do we see in different environments the same ratio of SRVs to Mira variables or does it depend on metallicity? Vassiliadis & Wood (Vassiliadis93 (1993)) calculated lifetimes of the major evolutionary phases for different initial masses and different metallicities. They found that higher metallicity will increase the lifetime of the early-AGB but decrease the lifetime on the TP-AGB. Miras stars populate the TP-AGB, therefore in environments with higher metallicities, such as the Galactic Bulge the lifetime of the TP-AGB is shorter and thus the number densities should decrease. This might explain the correlation between the ratio of SRVs to Miras and metallicity. However, while a large fraction of our variables on the TP-AGB are regular variables222In this case regular variables and Miras can be assumed to be identical. also semiregular variables are found. Lebzelter & Hron (LH99 (1999)) have shown that for stars in the solar neighborhood stellar evolution goes from SRVs to Miras. The semiregular stars found at a similar luminosity as the Miras (see Fig. 13) are therefore probably not in the same evolutionary state or they have different masses. Comparison of the number densities with expected lifetime is therefore problematic.

A lower metallicity leads also to a shift of the AGB towards higher temperatures in the HR diagram. The visual light change of these cool variables is dominated by highly temperature sensitive molecules like TiO (e.g. Reid & Goldston RG02 (2002)). If the stellar temperature is higher, these molecules will play a minor role. Lower metallicity will also make the TiO bands weaker. Therefore one would expect that the visual amplitudes will in general be smaller for lower metallicity. This would favour small amplitude variability in metal poor environments and would explain the smaller fraction of large amplitude objects in the LMC compared to the Bulge. It would also be consistent with the complete lack of Miras in metal poor globular clusters (Frogel & Whitelock FW98 (1998)).

However, one has to be extremely careful concerning possible selection effects, in particular for small amplitude variables. A homogeneous survey of variable stars in different Galactic environments is therefore needed.

###### Acknowledgements.

ALM thanks the EROS collaboration and in particular Jean-Baptiste Marquette for his help with the light curves production with the image subtraction method. ALM is extremely grateful to Claude Lamy who performs the tremendous work of sorting the whole EROS-1 data set. TL has been supported by the Austrian Science Fund under project number P14365-PHY. MS is supported by the Fonds zur Förderung der wissenschaftlichen Forschung (FWF), Austria, under the project number J1971-PHY. We wish to thank Josef Hron for fruitful discussion. Finally, we wish to thank the referee for constructive comments. This 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 cooperative agreement AST-8809616, and by the Mount Stromlo and Siding Spring Observatory, part of the Australian National University.


format: text

text:
# Crossover of quantum Loschmidt echo from golden rule decay to perturbation-independent decay

## I Introduction

The quantum Loschmidt echo, $`M(t)`$, measures the overlap of the evolution of the same initial state under two slightly different Hamiltonians in the classical limit,

$$M(t)=|\mathrm{\Phi }_0|\mathrm{exp}(iHt)\mathrm{exp}(iH_0t)|\mathrm{\Phi }_0|^2.$$
(1)
Here $`H_0`$ is the Hamiltonian of a classically chaotic system and

$$H=H_0+\kappa V,$$
(2)
with $`\kappa `$ being a small quantity.

This quantity characterizes the stability of quantum motion under a small change of the Hamiltonian, e.g. by the interaction with the environment, it is thus called “fidelity” and of great interest in the fast developing field of quantum information.

Moreover, a relationship between this quantity and the Lyapunov exponent that characterizes the classical chaos has been established analytically by using the semiclassical theory by Jalabert and Pastawski , which has been confirmed numerically in several models . Therefore, this quantity has attracted a great attention from the community of quantum chaos ,

It has been shown that several regimes exist. In the golden rule regime of the perturbation parameter, above a perturbative border, the fidelity has been found to have a simple decaying behavior, $`M(t)\mathrm{exp}(\mathrm{\Gamma }_Lt)`$ , where $`\mathrm{\Gamma }_L`$ is the half width of the local density of states (LDOS) $`\rho _L`$, which has a Lorentzian form (Breit-Wigner form),

$$\rho _L(E)=\frac{\mathrm{\Gamma }_L/2\pi }{E^2+\mathrm{\Gamma }_L^2/4},$$
(3)
with $`\mathrm{\Gamma }_L=2\pi \rho U^2`$. Here $`\rho `$ is the density of states and $`U`$ is the typical transition matrix element of $`\kappa V`$ between eigenstates of $`H_0`$. As the perturbation parameter $`\kappa `$ is increased, the exponential decay of $`M(t)`$ deviates from $`\mathrm{exp}(\mathrm{\Gamma }_Lt)`$, when $`\mathrm{\Gamma }_L`$ becomes comparable to the band width of $`H_0`$ and $`\rho _L`$ deviates from the Lorentzian form notably. A decay with the Lyapunov exponent $`\lambda `$ of the underlying classical chaotic dynamics, $`M(t)\mathrm{exp}(\lambda t)`$, can appear in this regime of the perturbation parameter, when the (effective) Planck constant is small enough and the initial state is chosen suitably, e.g., a narrow wavepacket . An interesting feature of this decay of $`M(t)`$ is that it is perturbation-independent, in the sense that it is irrelevant to the strength of the perturbation $`\kappa V`$ and is determined by the classical behavior of the system $`H_0`$ (and also that of $`H`$ due to the smallness of the perturbation in the classical limit). Based on the semiclassical theory and the random matrix theory analysis, a transition of $`M(t)`$ from the $`\mathrm{\Gamma }_L`$-decay to the $`\lambda `$-decay has been conjectured to occur at $`\mathrm{\Gamma }_L=\lambda `$ . Numerical results in several models support the conjecture .

The validity of results of the semiclassical theory depends on the value of the (effective) Planck constant. When the (effective) Planck constant is not small enough, the validity should be checked by direct quantum mechanical calculations. In fact, the behavior of the fidelity in this case is still not quite clear. For example, whether there is a sharp transition from the $`\mathrm{\Gamma }_L`$-decay to the $`\lambda `$-decay, or the transition occur in a finite regime of the perturbation parameter. It is even unclear whether $`M(t)`$ could decay with the Lyapunov exponent. Indeed, a perturbation-independent but slower than $`\mathrm{exp}(\lambda t)`$ decay has been observed in the kicked top model , the mechanism of which is still not clear. Meanwhile, the random matrix theory seems unsuitable for such problems, since perturbation-independent decay of $`M(t)`$ is usually initial-state-dependent, the feature of which is hard to be captured by the random matrix theory treatment.

The quantity used in this paper in analyzing properties of the fidelity $`M(t)`$ is its Fourier transform, denoted by $`F_M(E)`$ in what follows. As will be shown in Sect. II that $`F_M(E)`$ is determined directly by properties of the initial state and of the (quasi)energy eigenstates of the two classically chaotic systems concerned. The role played by $`F_M(E)`$ in understanding the behavior of $`M(t)`$ is similar to that by the LDOS for the survival probability $`P(t)=|\alpha |\mathrm{exp}(iHt)|\alpha |^2`$.

In this paper, properties of $`F_M(E)`$ and their relations to those of $`M(t)`$ will be investigated numerically in the kicked top model . The fidelity in this model has been studied in Refs.  and here we will concentrate on issues unaddressed in the previous work. Specifically, in addition to properties of $`M(t)`$ related to the function $`F_M(E)`$, we focus ourself on the following problems:

(i) The mechanism of the perturbation-independent decay slower than $`\mathrm{exp}(\lambda t)`$.

(ii) Due to the small difference between $`H_0`$ and $`H`$ in the classical limit, the classical counterpart of the fidelity $`M(t)`$ should change quite slowly when $`t`$ is small enough (see Ref.  for numerical results in another model). In the quantum mechanical case, a similar phenomenon should be of practical interest in quantum computing. However, it was not clear whether or not the phenomenon exists.

The paper is organized as the follows. The function $`F_M`$ is introduced in Sect. II. In Sect.III, we shall present our numerical investigation of the problems mentioned above, in particular, for properties of the amplitude and the phase of the function $`F_M(E)`$ with initial coherent states. The perturbation-independent decay slower than $`\mathrm{exp}(\lambda t)`$ will be shown to be due to large quantum fluctuations in individual fidelity, which are caused by the coherence possessed the initial coherent states, when the (effective) Planck constant is not small enough. It will be shown that the fidelity $`M(t)`$ does not decay obviously within a short initial time, the length of which is determined by some properties of the function $`F_M(E)`$. The relationship between the decaying rate of $`M(t)`$ after the short initial time and some properties of the function $`F_M(E)`$ will also be shown numerically. Conclusions and discussions will be given in Sect. IV.

## II Function $`F_M(E)`$ as the Fourier transform of $`M(t)`$

In this paper, the eigenstates of the two Hamiltonians $`H_0`$ and $`H`$ are denoted by $`|\alpha `$ and $`|\beta `$, respectively, with eigenenergies $`E_\alpha `$ and $`E_\beta `$,

$$H_0|\alpha =E_\alpha |\alpha ,H|\beta =E_\beta |\beta .$$
(4)
The expanding coefficient of $`|\beta `$ in $`|\alpha `$ is indicated by $`C_{\beta \alpha }`$, $`C_{\beta \alpha }=\alpha |\beta `$. The LDOS $`\rho _L(E)`$ of an eigenstate $`|\alpha `$ of $`H_0`$ is determined directly by the eigensolutions the two Hamiltonians concerned,

$$\rho _L(E)=\underset{\beta }{}|C_{\beta \alpha }|^2\delta \left(E(E_\beta E_\alpha )\right).$$
(5)
As is known, when the initial state $`|\mathrm{\Phi }_0`$ is an eigenstate of $`H_0`$, the fidelity amplitude $`m(t)`$,

$$m(t)=\mathrm{\Phi }_0|\mathrm{exp}(iHt)\mathrm{exp}(iH_0t)|\mathrm{\Phi }_0,$$
(6)
is just the Fourier transform of the LDOS $`\rho _L(E)`$, and the form of the LDOS is important in the study of $`M(t)`$, which is just the survival probability $`P(t)`$ in this case. However, as shown in Ref. , the LDOS can not explain the perturbation-independent decay of $`M(t)`$. In fact, the perturbation-independent decay of $`M(t)`$ appears only for some special set of initial states, that is, it is initial-state-dependent. Although a rigorous condition is still lacking, a sufficient condition for such initial states is expressed in a rough way that an initial state should be a narrow wavepacket. When the decay of $`M(t)`$ is initial-state-dependent, one should use a quantity that is more general than the LDOS in analyzing properties of $`M(t)`$.

The function $`F_M`$, the Fourier transform of the fidelity $`M(t)`$, is a suitable candidate, which is introduced as the follows. Let us first express the fidelity amplitude $`m(t)`$ in the form

$$m(t)=f_m(ϵ)e^{iϵt}𝑑ϵ,$$
(7)
where

$`f_m(ϵ)={\displaystyle \underset{\alpha ,\beta }{}}A_{\alpha \beta }\delta \left[ϵ(E_\alpha E_\beta )\right],`$ (8)
$`A_{\alpha \beta }=\mathrm{\Phi }_0|\beta \alpha |\mathrm{\Phi }_0C_{\beta \alpha }^{}.`$ (9)
Note that the quantity $`A_{\alpha \beta }`$ defined in Eq. (9) is invariant under the transformation

$$|\alpha e^{i\phi _\alpha }|\alpha ,|\beta e^{i\phi _\beta }|\beta ,|\mathrm{\Phi }_0e^{i\phi _0}|\mathrm{\Phi }_0,$$
(10)
where $`\phi _\alpha ,\phi _\beta `$, and $`\phi _0`$ are arbitrary phases. That is, it is irrelevant to the relative phases among the states $`|\alpha `$, $`|\beta `$, and $`|\mathrm{\Phi }_0`$. The fidelity $`M(t)`$ can be expressed as

$$M(t)=|m(t)|^2=F_M(E)e^{iEt}𝑑E,$$
(11)
where

$$F_M(E)=f_m(ϵ)f_m^{}(ϵE)𝑑ϵ.$$
(12)
The quantity $`F_M(E)`$ is in fact the correlation function of $`f_m(ϵ)`$. Note that $`F_M(E)`$ can be calculated directly from properties of the states $`|\alpha `$, $`|\beta `$, and the initial state $`|\mathrm{\Phi }_0`$. In the following section, we will investigate numerically the relationship between properties of $`M(t)`$ and those of $`F_M(E)`$. To this end, it is convenient to write $`F_M(E)`$ as

$$F_M(E)=\rho _F(E)e^{i\theta _F(E)}.$$
(13)
Some arguments can be given to the reason why the random matrix theory can be used to predict the decaying behavior of the fidelity $`M(t)`$ in the golden rule regime, but not in the regime where the perturbation-independent decay of $`M(t)`$ appears. The initial states, whose fidelity is found to have perturbation-independent exponential decay controlled by the Lyapunov exponent, such as narrow wave packets in the configuration space , coherent states , momentum eigenstates , etc., usually have wide spreading in the chaotic states $`|\alpha `$ ($`|\beta `$). Such initial states, which are regular in the classical limit, possess certain kind of coherence. It must be pointed out that the semiclassical theory does not predict the perturbation-independent decay of the fidelity for an arbitrary initial state . As shown in the previous paragraph, the function $`F_M(E)`$, equivalently, $`M(t)`$, is basically determined by the quantity $`A_{\alpha \beta }`$ in Eq. (9). In the golden rule regime, the LDOS $`\rho _L(E)`$ is narrow, i.e., $`C_{\beta \alpha }`$ has many small components; as a result, for the quantity $`A_{\alpha \beta }`$, components of the initial states in the chaotic states, namely, $`\mathrm{\Phi }_0|\beta \alpha |\mathrm{\Phi }_0`$, are not equally effective in the whole (quasi)energy region of $`H_0`$ and $`H`$. This should suppress the coherence of the initial states, which makes it possible to treat the components of the initial state taking part in the evolution of $`M(t)`$ effectively as random numbers, like in the random matrix theory treatment in Refs. . On the other hand, when the perturbation is so strong that the width of the LDOS is comparable to the (quasi)energy band width of $`H_0`$, the coherence possessed by the initial states should play a non-negligible role in the evolution of the fidelity. The random matrix theory can not be applied in this case, since it cannot describe the coherence in the initial states.

## III  Numerical Study of Fidelity and Its Fourier Transform

### A  The model

The unperturbed Hamiltonian $`H_0`$ of the kicked top model used in this paper is

$$H_0=\frac{\pi }{2\tau }S_y+\frac{K}{2S}S_z^2\underset{n}{}\delta (tn\tau ),$$
(14)
where $`S`$ is the total angular momentum, $`\tau `$ is the period, and $`K`$ is a parameter adjusting the strength of the kicks. Without loss of generality, the period $`\tau `$ is set to be unit, $`\tau =1`$. The model describes a vector spin that undergoes a free precession around the $`y`$-axis and is periodically perturbed by kicks around the $`z`$-axis. The time evolution of an initial state $`|\mathrm{\Phi }_0`$ at $`t=0`$ is governed by the Floquet operator

$$F_0=\mathrm{exp}[i\frac{K}{2S}S_z^2]\mathrm{exp}[i\frac{\pi }{2}S_y],$$
(15)
where $`\mathrm{}`$ has been set to be unit. The classical limit of the system, which is obtained by letting $`S\mathrm{}`$ with $`1/S`$ serving as the effective Planck constant, is fully chaotic for $`K9`$.

The perturbation $`\kappa V`$ is chosen in the same way as in Ref. , i.e., a slightly delayed periodic rotation of constant angle around the $`x`$-axis,

$$V=\frac{\pi }{2\tau }S_x\underset{n}{}\delta (tn\tau ϵ).$$
(16)
The Floquet operator of the system $`H`$ is

$$F=\mathrm{exp}[i\kappa \frac{\pi }{2}S_x]F_0,$$
(17)
which gives the time evolution of an initial state,

$$|\mathrm{\Phi }(t)=F^t|\mathrm{\Phi }_0,$$
(18)
where $`t=0,1,2,\mathrm{}`$. For this model, discussions in Sect. II are still valid, with the eigenstates of $`H_0`$ and $`H`$ changed to those of the Floquet operators $`F_0`$ and $`F`$, respectively, with the corresponding eigenenergies changed to quasi-energies.

In numerical calculations, for the sake of convenience, the perturbation parameter $`\kappa `$ is written as

$$\kappa =j\times 10^3.$$
(19)
Unless addressed explicitly, the value of $`S`$ is 500 for the numerical calculations discussed in this paper. Since the time $`t`$ in the evolution equation (18) takes integer values only, according to Eq. (11), the domain of the variable $`E`$ in the function $`F_M(E)`$, which is $`[4\pi ,4\pi ]`$ in the general case, can be reduced to $`[\pi ,\pi ]`$, with the definition of $`F_M(E)`$ changed accordingly, namely,

$$F_M^{\mathrm{fold}}(E)=\underset{n}{}F_M(E+2n\pi ),$$
(20)
where $`n`$ takes the possible values in $`(0,\pm 1,\pm 2)`$. The function $`F_M^{\mathrm{fold}}(E)`$ is more closely related to properties of $`M(t)`$ than the basic one in Eq. (12). It is this function that will be used in what follows in the numerical investigation of the kicked top model and, for brevity, it will be denoted by $`F_M(E)`$.

Initial states studied in this paper are in most cases coherent states of the SU(2) group ,

$`|z=Ae^{z^{}S_+}|M=S`$ (21)
$`z=tan{\displaystyle \frac{\theta }{2}}e^{i\varphi },`$ (22)
where $`|M`$ is the eigenvector of $`S_z`$, with the eigenvalue $`M`$ ($`M=S,S+1,\mathrm{},S`$), and $`A`$ is the normalization coefficient. The coherence in a coherent state $`|z`$ can be easily seen when it is expanded in the states $`|M`$. However, in the study of the fidelity, $`|z`$ is required to be expanded in the eigenstates of the Floquet operators $`F_0`$ and $`F`$. Although the expanding coefficients in such chaotic states must contain information on the coherence possessed by the coherent state, manifestation of the information is not so easy. In Sect. III D, we will show numerically that the function $`F_M(E)`$ introduced in Sect. II can supply some of the information.

### B  Perturbation-independent decay of fidelity and quantum fluctuations

Basic properties of the fidelity $`M(t)`$ in the kicked top model have been studied in Ref. . Here, as mentioned in the introduction, we are interested in a perturbation-independent decay of the averaged $`M(t)`$ with initial coherent states, which decays more slowly than the exponential decay predicted by the semiclassical theory. In order to understand this phenomenon, we study the individual $`M(t)`$. In the regime of the perturbation parameter where the perturbation-independent decay appears, individual $`M(t)`$ is found to have large fluctuations. See Fig. 1 for some examples of $`K=13.1`$ and $`j=6`$, where $`M(t)`$ are plotted in the logarithm scale. On the other hand, in the golden rule regime, e.g., $`K=13.1`$ and $`j=1`$, no obvious fluctuation in individual $`M(t)`$ is observed. The large fluctuations are in fact quantum effects caused by the coherence in the coherent states. According to the arguments given in Sect. II, they should appear when the expanding coefficients of the states $`|\beta `$ in $`|\alpha `$ spread in the whole quasi-energy band of $`F_0`$ effectively, that is, when the half width of the LDOS becomes comparable to the bandwidth of the quasi-energy, which is found numerically for $`j3`$. In the classical limit $`S\mathrm{}`$, the fluctuations should disappear.

As shown in Fig. 1, each $`M(t)`$ has a ‘shoulder’ in a short initial time, denoted by $`t_s`$ in what follows, within which it changes slowly. The exponential-type decay of $`M(t)`$ appears only after $`t_s`$. In order to study the average decaying behavior of $`M(t)`$, their ‘shoulders’ should be subtracted. To this end, we shift the time variable $`t`$ to $`t_d=tt_s`$ and take average over $`M(t_d)`$. Here we meet a problem, since individual $`M(t)`$ has large fluctuations. The averaged $`M(t_d)`$ may be mainly determined by a small fraction of the $`M(t_d)`$ taken for averaging, if they have extraordinarily large values. Therefore, in addition to the standard averaged $`M(t_d)`$, denoted by $`M_1(t_d)`$,

$$M_1(t_d)=\frac{1}{N}\underset{i=1}{\overset{N}{}}M(t_d,\mathrm{\Phi }_i),$$
(23)
we also calculate another averaged $`M(t)`$, denoted by $`M_2(t_d)`$, by taking the logarithm of $`M(t_d)`$ before performing the summation,

$$\mathrm{ln}M_2(t_d)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\mathrm{ln}M(t_d,\mathrm{\Phi }_i).$$
(24)
Here $`M(t,\mathrm{\Phi }_i)`$ is the fidelity of an initial coherent state $`|\mathrm{\Phi }_i`$ defined in Eq. (1), with the dependence on $`|\mathrm{\Phi }_i`$ written explicitly.

Numerical results for $`M_1(t_d)`$ (upper) and $`M_2(t_d)`$ (bottom) are shown in Fig. 2, with increasing $`\kappa `$ ($`j`$ from 1 to 8), when $`K=13.1`$. In calculating the two quantities, $`t_s`$ for each $`M(t)`$ of $`j=1`$ and 1.5 was taken to be the first $`t`$ at which $`M(t+1)<0.9`$. For $`j2`$, $`t_s`$ was determined by the first $`t`$ at which $`M(t)/M(t+1)>1.5`$. Only states satisfying $`M(t_s)>0.85`$ were used in averaging. The reason of taking two methods in calculating $`t_s`$ is that $`M(t)`$ decays slowly at $`j=1`$ and 1.5. The number $`N`$ of coherent states taken for averaging is 1000. Results in Fig. 2 for $`M_1(t_d)`$ are in consistence with those in Ref. , namely, the decay of $`M_1(t_d)`$ saturates at an exponential decay which is slower than the one with the Lyapunov exponent. However, when the large fluctuations in individual $`M(t)`$ are suppressed by the second averaging procedure, Fig. 2 shows that the decay of the averaged fidelity $`M_2(t_d)`$ saturates at an exponential decay that is close to the Lyapunov decay. The deviation of the saturation decay of $`M_1(t_d)`$ from the one predicted by the semiclassical theory is caused by the large fluctuations in individual $`M(t)`$. In fact, this can also be seen roughly in Fig. 1.

To show the difference between the decay of $`M_1(t_d)`$ and that of $`M_2(t_d)`$ in Fig. 2 in a quantitative way, we calculate their decaying exponents $`\gamma _1`$ and $`\gamma _2`$,

$$M_1(t_d)\mathrm{exp}(\gamma _1t_d)M_2(t_d)\mathrm{exp}(\gamma _2t_d),$$
(25)
the values of which are presented in Fig. 3. Numerically, $`\gamma _i`$ ($`i=1,2`$) are calculated by the best linear fitting to $`\mathrm{ln}M_i(t_d)`$, with the fitting lines fixed to the values of $`\mathrm{ln}M_i(0)`$ at $`t_d=0`$. The points $`t_d`$ used in fitting are those satisfying $`1M_2(t_d)0.008`$. In fig. 3, we see that, in the golden rule regime ($`j=1`$ and 1.5), $`\gamma _1`$ is close to $`\gamma _2`$ because of the small fluctuation in individual $`M(t)`$, and both of them are close to the half width of LDOS, $`\mathrm{\Gamma }_L`$, as predicted by the random matrix theory. The difference of the three quantities becomes obvious, when $`j2`$. In the parameter regime $`j6`$, $`\gamma _2`$ fluctuates around the Lyapunov exponent $`\lambda =16.5`$, while $`\gamma _1`$ is obviously smaller than $`\lambda `$. The value of $`\gamma _2`$ at $`j=14`$ is obviously larger than the Lyapunov exponent. This may be due to the reason that the perturbation cannot be regarded as small at this perturbation value. In fact, when the perturbation is so strong that the classical perturbation theory breaks down, $`M(t)`$ may decay faster than the Lyapunov decay . As discussed above, when $`S`$ is increased and we are more and more close to the classical limit, the fluctuations in $`M(t)`$ shown in Fig. 1 should become small and small, which would make $`\gamma _1`$ approach to $`\gamma _2`$. Indeed, we calculate the case of $`S=1000`$ ($`\kappa `$ reduced to half to those of $`S=500`$), where $`\gamma _2`$ are found to be fluctuating around $`\lambda `$ as well and $`\gamma _1`$ a little larger than the values at $`S=500`$. However, $`\gamma _1`$ of $`S=1000`$ are still obviously smaller than $`\lambda `$, showing that $`S=1000`$ is still not in the deep semiclassical regime.

### C Dependence of fidelity on initial states

In different regimes of the perturbation parameter, the initial state can play different roles in influencing the behavior of the fidelity. In the Golden rule regime, the LDOS is narrow and, as discussed in Sect. II, the random matrix theory can be used to predict the behavior of the fidelity. In this regime, there should be no difference between the decaying behavior of the fidelity of an initial coherent state and that of an arbitrary initial state, that is, the decay of the fidelity should be initial-state-independent. Indeed, in a study of the influence of sub-Planck scale structures  on the fidelity in Ref. , the averaged fidelity of the initial states

$$|\mathrm{\Phi }_0(T)=\mathrm{exp}(iH_0T)|\mathrm{\Phi }_0^c,$$
(26)
where $`|\mathrm{\Phi }_0^c`$ are coherent states, has been found numerically to be independent on the parameter $`T`$ and almost the same as that of initial coherent states. We compare the averaged fidelity of initial coherent states with that of initial random states of the form

$$|\mathrm{\Phi }_0^{\mathrm{random}}=\underset{\alpha }{}c_\alpha |\alpha ,$$
(27)
where $`c_\alpha `$ are random complex numbers satisfying the normalization condition. It is found that the former ones have ‘shoulders’, while the latter ones do not. With the ‘shoulders’ of the former ones subtracted, the decay of the corresponding averaged fidelity is found to be close to each other in the golden rule regime. Some of the results are plotted in Fig. 4, where 1000 initial states were used in averaging for each case.

In the perturbation-independent regime of the fidelity of initial coherent states, $`M(t)`$ of $`|\mathrm{\Phi }_0(T)`$ is found to be decaying in part of the evolution time in a way similar to that of initial coherent states . This implies that, although the system $`H_0`$ is classically chaotic, the time evolution operator $`\mathrm{exp}(iH_0T)`$ does not destroy the coherence possessed by a coherent state. Indeed, we find that random superpositions of coherent states defined by

$$|\mathrm{\Phi }_0^{\mathrm{rs}}(n_c)=\underset{n=1}{\overset{n_c}{}}e^{i\phi _n}|z_n,$$
(28)
where $`\phi _n`$ are random phases and $`|z_n`$ are coherent states in Eq. (21), still have some coherence properties of the coherent states. The averaged fidelity with $`|\mathrm{\Phi }_0^{\mathrm{rs}}(n_c)`$ being initial states at the perturbation parameter $`j=6`$ are presented in Fig. 5. We see that the averaged $`M(t)`$ with $`n_c>1`$ has a rapid decrease between $`t=0`$ and 1, which is faster than the Lyapunov decay for $`n_c8`$. This indicates that part of the coherence in the coherent components $`|z`$ of the initial states $`|\mathrm{\Phi }_0^{\mathrm{rs}}(n_c)`$ is destroyed by the random phases in the initial states. In fact, when $`\alpha |\mathrm{\Phi }_0`$, the components of initial states in the eigenstates $`|\alpha `$, can be treated as random numbers, the random matrix theory analysis predicts that the decay of the averaged $`M(t)`$ is controlled by the averaged LDOS , the half width of which is larger than the Lyapunov exponent at $`j=6`$ (Fig. 3).

After the initial transient rapid decrease, the remanent coherence in the initial states $`|\mathrm{\Phi }_0^{\mathrm{rs}}(n_c)`$ plays a role similar to that of initial coherent states. Indeed, the averaged $`M(t)`$ changes slowly in some short time intervals after $`t=1`$, which are analogues of the ‘shoulders’ of the fidelity of initial coherent states, then, decreases exponentially in a way similar to that of initial coherent states, until they become close to the saturation value, which is about $`1/(2S)`$. The larger the value of $`n_c`$ is, the smaller the remanent coherence will be. The fidelity of random superpositions of $`|\alpha `$ in Eq. (27) decays faster than the Lyapunov decay shown by the dotted straight line in Fig. 5, as predicted by the random matrix theory.

### D Properties of $`F_M(E)`$ and relation to properties of $`M(t)`$

As shown in Sect. II, the Fourier transform of $`M(t)`$, the function $`F_M(E)`$, is directly determined by the states $`|\alpha `$, $`|\beta `$ and $`|\mathrm{\Phi }_0`$. In this subsection, we study properties of the function $`F_M(E)`$ in the kicked top model numerically, in particular, its absolute value part $`\rho _F(E)`$ and its phase part $`\theta _F(E)`$. Unless addressed explicitly, initial states $`|\mathrm{\Phi }_0`$ studied in this subsection are coherent states. It will be shown that properties of the central part of $`F_M(E)`$ are closely related to the decaying rate and the width of the ‘shoulder’ of the averaged fidelity $`M(t)`$.

Typical shapes of $`\rho _F(E)`$ and $`\theta _F(E)`$ in cases of different values of the perturbation parameter are shown in Fig. 6, where the same coherent state is used. For $`\rho _F(E)`$ of $`j=1`$ and 2, the values of $`\rho _F`$ are plotted in the logarithm scale. Each $`\rho _F(E)`$, symmetric with respect to $`E=0`$, is composed of three or more parts, with a main central part. A manifestation of the coherence in the coherent state $`|\mathrm{\Phi }_0`$ is that the corresponding phase functions $`\theta _F(E)`$ have related divisions, as can be seen in Fig. 6. (If $`|\mathrm{\Phi }_0`$ is taken as a random state $`|\mathrm{\Phi }_0^{\mathrm{random}}`$ in Eq. (27), values of the function $`\theta _F(E)`$ will be close zero.)

The central part of each $`\rho _F`$ is found to be fitted well by the Lorentzian form, with a width denoted by $`\mathrm{\Gamma }_F`$, to be distinguished from $`\mathrm{\Gamma }_L`$ for the width of the LDOS. See Fig. 7 for an example of the Lorentzian fit to the central part of the $`\rho _F`$ in Fig. 6 of j=4. The Lorentzian form used in fitting is $`\rho _L(E)+a`$, where $`\rho _L(E)`$ is the function on the right hand side of Eq. (3) and $`a`$ is a fitting constant. Therefore, the closeness of the central part of $`\rho _F(E)`$ to the fitting Lorentzian form does not predict directly an approximate exponential decay of $`M(t)`$ with an exponent $`\mathrm{\Gamma }_F`$. In fact, the contribution of the central part of $`F_M(E)`$ to $`M(t)`$ in Eq. (11) gives both positive and negative results. The final positive results for $`M(t)`$ come out, only when the contributions of the other parts of $`F_M`$ are taken into account. However, the average values of $`\mathrm{\Gamma }_F`$, which are shown by solid squares in Fig. 3, are found numerically to be close to the values of $`\gamma _1`$ for $`M_1(t_d)`$, in the whole parameter regime of $`j`$ from 1 to 14. When $`S`$ is taken equal to 1000, the change of the corresponding results of $`\mathrm{\Gamma }_F`$ are found to be smaller than that of $`\gamma _1`$. Since $`\gamma _1`$ increases only a little, when $`S`$ is changed from 500 to 1000, it is not clear at the present stage whether or not $`\mathrm{\Gamma }_F`$ can approach to $`\gamma _1`$ in the classical limit.

Numerical results for the phase function $`\theta _F(E)`$ show that it is approximately linear in the region of $`E`$ corresponding to the central part of $`\rho _F`$ (see Fig. 6 for examples),

$$\theta _F(E)k_\theta E,$$
(29)
where $`k_\theta `$ is the slope of the approximate linear behavior. Distributions of $`k_\theta `$ are shown in Fig. 8 for some perturbation parameters. The positions of the peaks of the distributions are close to $`k_\theta =1`$, while the distribution of $`j=4`$ has another small peak at $`k_\theta =1.5`$. To illustrate the influence of $`\theta _F(E)`$ on behaviors of $`M(t)`$ clearly, let us write $`M(t)`$ in the form

$$M(t)=\rho _F(E)e^{i\theta _t(E)}𝑑E,$$
(30)
where $`\theta _t(E)=\theta _F(E)Et`$. Some examples of the behavior of $`\theta _t(E)`$ are shown in Fig. 9. With increasing $`t`$, the approximate slope of $`\theta _t(E)`$ in the central region of $`F_M(E)`$ changes from positive to negative. For $`t2`$, the larger $`t`$ is, the steeper $`\theta _t(E)`$ will be. As shown below, this feature of $`\theta _t(E)`$ can be used to estimate the width of the ‘shoulder’ of $`M(t)`$.

In order to give an estimation of the ‘shoulder’ width of $`M(t)`$, we note that the exponential decay of $`M(t)`$ should appear when $`\theta _t(E)`$ is steep enough, specifically, when $`\theta _t(W_F/2)`$ is comparable to $`\pi `$, where $`W_F`$ is the width of the region of $`E`$ occupied by the main central part of $`\rho _F`$. Quantitatively, we use $`t_{\mathrm{s}\theta }`$ to denote the average of the largest $`t`$ satisfying

$$|k_\theta t|\frac{W_F}{2}R_c\pi ,$$
(31)
where $`R_c`$ is a quantity used to show the closeness to $`\pi `$. The average width of the ‘shoulders’ of $`M(t)`$, denoted by $`t_{\mathrm{sM}}`$, are calculated from the averaged $`M(t)`$ directly. Since $`t`$ has integer values only in the fidelity $`M(t)`$ of the kicked top model, we use the linear interpolation in calculating $`t_{\mathrm{sM}}`$ by $`M(t_{\mathrm{sM}})=M_c`$, where $`M_c`$ is a quantity measuring the lower border of the ‘shoulder’ of the averaged $`M(t)`$. In studying the relation between $`t_{\mathrm{sM}}`$ and $`t_{\mathrm{s}\theta }`$, for a given value of $`M_c`$, the value of $`R_c`$ is determined by the smallest value of

$$\underset{j=1}{\overset{8}{}}(t_{\mathrm{s}\theta }t_{sM})^2,$$
(32)
in the parameter regime $`1j8`$. Figure 10 shows two examples of $`M_c=0.9`$ with $`R_c0.75`$ and $`M_c=0.8`$ with $`R_c0.9`$, where we see that $`t_{\mathrm{s}\theta }t_{\mathrm{sM}}`$ for $`j`$ between 1 and 8. The difference between $`t_{\mathrm{s}\theta }`$ and $`t_{\mathrm{sM}}`$ becomes increasingly large, when $`j`$ is larger than 8, which is related to the fact that the perturbation $`\kappa V`$ is not quite small for $`j>8`$. Variation of $`R_c`$ with respect to $`M_c`$ is presented in Fig. 11, where $`R_c`$ is seen to be almost linear with $`M_c`$ for $`M_c`$ between 0.7 and 0.95.

## IV Conclusions and Discussions

In this paper, we study the crossover of the fidelity from the golden rule regime to the perturbation-independent decay regime. We show that the deviation of the perturbation-independent decay of the averaged fidelity from the prediction of the semiclassical theory is due to large quantum fluctuations in individual fidelity caused by the coherence in the initial coherent states. When the quantum fluctuations are suppressed by an appropriate averaging procedure, the perturbation-independent decay of the fidelity is found to be close to the semiclassical prediction.

We also show that when the effective Planck constant is not small enough to guarantee the validity of the semiclassical theory, the crossover does not occur at a point of the perturbation parameter, but takes place in a region of the parameter (Fig. 3). The crossover in this case can be analyzed by neither the semiclassical theory, because the effective Planck constant is not small enough, nor the random matrix theory, since the perturbation-independent decay of the fidelity is initial-state-dependent.

Furthermore, a function $`F_M(E)`$ that is in fact the Fourier transform of the fidelity, is introduced to study the fidelity in the crossover region. This function is determined directly by the initial state and the eigenstates of the two Floquet operators concerned. Numerically, for initial coherent states, the absolute value part and the phase part of the function $`F_M(E)`$ are found to have related division into several sub-parts, which is a manifestation of the coherence possessed by the initial coherent states. In the whole crossover region, properties of the central part of $`F_M(E)`$ are found to be closely related to both the exponent of the exponential decay of the averaged fidelity and the length of the short initial time period, within which the fidelity almost does not change due to the coherence possessed by the initial coherent state and the smallness of the perturbation in the classical limit.

The function $`F_M(E)`$ is studied mainly numerically in this paper. A detailed analytical analysis of the function, e.g., by making use of its expression as the correlation function of $`f_m(ϵ)`$ given in Sect. II, would supply further understanding of the behaviors of the fidelity, in analogue to the relation between the LDOS and the survival probability. Information of the coherence in a coherent state may be hidden in its expanding coefficients in the eigenstates of classically chaotic systems. Numerical results in this paper show that the function $`F_M`$ supplies a useful method of extracting such information. We would like to mention that some other quantities discussed in this paper may serve this purpose as well in other situations, e.g., the quantity $`A_{\alpha \beta }`$ introduced in Sect. II, which is irrelevant to the arbitrary phases that can be given to the initial coherent state and the quasi-energy eigenstates of the two classically chaotic systems.

## ACKNOWLEDGMENTS

The authors are grateful to G. Casati for valuable discussions. The work was supported in part by the Academic Research Fund of the National University of Singapore and the DSTA of Singapore. We also thank Dr. Wang Jiao for providing the programme for diagonalization of unitary matrix.


text:
# 1 Introduction

## 1 Introduction

Instabilities during growth of surfaces are of crucial importance for fabrication of devices . Of particular interest in recent time is the instability of the two-dimensional (2D) layer-by-layer growth against the formation of coherently strained (dislocation-free) three-dimensional (3D) islands of nanometer scale in highly mismatched epitaxy. The latter is known as a “coherent Stranski-Krastanov” (SK) growth , and is a subject of intense research owing to possible optoelectronic applications as lasers and light emiting diodes \[2-4\]. That is why much effort has been made in the last decade to determine the equilibrium shape of the crystallites as a function of the volume \[5-8\], the change of shape during growth , the strain distribution within the coherent islands and their energy \[10-20\], the kinetics of growth of arrays of quantum dots and the physical reason of the narrow size distribution \[21-23\], which is often experimentally observed.

The physical reason of occurrence of the Stranski-Krastanov growth mode is generally inderstood. Too much strain energy accumulates into the film during the initial planar growth, and the strong adhesion exerted by the substrate (which is the reason for the planar growth) disappears beyond several atomic diameters. A wetting layer (WL) composed of an integer number of equally strained monolayers is thus formed. The growth continues further by the formation of 3D crystallites, in which the additional surface energy is overcompensated by the strain relaxation. In other words, the higher energy phase representing a homogeneously strained planar film is replaced beyond some critical thickness by a lower energy phase of (completely or partially) relaxed 3D crystallites.

Although the essential physics seems clear, too many questions of fundamental character remain to be answered. As the atoms on top of the surface of the wetting layer do not “feel” energetically the presence of the substrate and both the wetting layer and the 3D islands consist of one and the same material, we can consider as a first approximation the formation of coherent 3D clusters in SK growth as homoepitaxial growth on an uniformly strained crystal surface. If so, it is not clear what is the thermodynamic driving force for 3D islanding if the islands are coherently strained to the same degree as the underlying wetting layer. This question is closely connected with the structure and energy of the boundary between the 3D islands and the wetting layer. The energy of this boundary is often taken equal to zero . This means a complete wetting of the 3D islands by the substrate (the WL) which rules out the 3D islanding from a thermodynamic point of view. It is also not clear why coherent 3D islands are observed in compressed rather than in expanded overlayers, and at values of the misfit $`\epsilon _0=\mathrm{\Delta }a/a`$ that are huge for materials with directional and brittle covalent bonds (InAs/InP (3.2%) , Ge/Si (4.2%) , InAs/GaAs (7.2%) , CdSe/ZnSe (7.6%) ). The only exception, to the authors’ knowledge, of expanded overlayer, is the system PbSe/PbTe (-5.5%) . Other question is whether the misfit should be greater than some critical value in order for the coherent 3D islanding to take place. Are two-dimensional monolayer height islands necessary precursors for the formation of 3D islands as suggested by some authors ? If yes, is there a critical volume size (or a size of the 2D island) for the 2D-3D transformation to occur? What is the pathway of the latter, does it pass through a series of intermediate states with increasing thickness, and are these states stable or metastable? In this paper we make an attempt to answer at least qualitatively some of the questions posed above.

The thermodynamic driving force for occurrence of one or another mode of growth should be given by the difference $`\mathrm{\Delta }\mu =\mu (n)\mu _{3D}^0`$ of the chemical potentials $`\mu (n)`$ of the film, and $`\mu _{3D}^0`$ of the bulk crystal of the same material. The film chemical potential depends on the thickness measured in number $`n`$ of monolayers owing to the thickness distribution of the misfit strain and the attenuation of the energetic influence of the substrate \[31-33\]. If we deposit a crystal $`A`$ on the surface of a crystal $`B`$ the thermodynamic driving force can be written in terms of interatomic energies $`\mathrm{\Delta }\mu =E_{AA}\Phi `$ where $`\Phi =1E_{AB}/E_{AA}`$ is the so-called adhesion parameter which accounts for the wetting of the substrate by the overgrowth . $`E_{AA}`$ and $`E_{AB}`$ are the energies per atom to disjoin a half-crystal $`A`$ from a like half-crystal $`A`$ and from an unlike half-crystal $`B`$, respectively. $`E_{AB}`$ is in fact the adhesion energy which includes in itself the thickness distribution of the strain energy due to the lattice misfit, and the attenuation of the bonding with the substrate . The adhesion parameter $`\Phi `$ is the same which accounts for the influence of the substrate on the work of formation of 3D nuclei of different material on top of it in the classical nucleation theory . Replacing the bonding energies $`E_{AA}`$ and $`E_{AB}`$ by the corresponding surface energies gives the famous 3-$`\sigma `$ criterion of Bauer for the mode of growth $`\mathrm{\Delta }\mu =a^2[\sigma _A+\sigma _{AB}(n)\sigma _B]`$ , where $`a^2`$ is the area occupied by an atom at the interface.

The thickness dependence of the film chemical potential is schematically illustrated in Fig. 1. In the two limiting cases of Volmer-Weber (VW) (incomplete wetting, $`0<\Phi <1`$) and Frank-van der Merwe (FM) growth (complete wetting, $`\Phi 0,\epsilon _00`$) $`\mu (n)\mu _{3D}^0`$ goes asymptotically to zero from above and from below, respectively, but changes sign in the case of SK growth ($`\Phi 0,\epsilon _00`$. In the latter case, beyond the maximum, we consider the 3D islands as the overlayer material $`A`$, and the wetting layer as the substrate crystal $`B`$. Thus the strained wetting layer and the relaxed 3D islands represent necessarilly different phases in the sense of Gibbs. The wetting layer can be in equilibrium only with an undersaturated vapor phase, whereas the 3D islands are in equilibrium with a supersaturated vapor. The dividing line is $`\mathrm{\Delta }\mu =0`$ at which the wetting layer cannot grow thicker and the 3D islands cannot nucleate and grow. Thus the adhesion parameter $`\Phi =\mathrm{\Delta }\mu /E_{AA}`$ relative to the cohesion energy $`E_{AA}`$ is in fact equal to the thermodynamic driving force for the occurence of one or another mode of growth. In other words, we can treat the SK mode as a FM mode driven by complete wetting ($`\mathrm{\Delta }\mu <0`$), followed by VW mode driven by incomplete wetting ($`\mathrm{\Delta }\mu >0`$). (The more rigorous definition is $`d\mu /dn<0`$ or $`d\mu /dn>0`$ ). The question is how the lattice misfit can lead to incomplete wetting ($`\Phi >0`$) on the surface of the wetting layer if the energetic influence exerted by the substrate is already lost, i.e. $`E_{AB}E_{AA}`$. In the classical SK mode the incomplete wetting is due to the introduction of misfit dislocations (MDs) . Once we know the answer of this question in the case of a coherent SK growth we could easily find the answers of the others.

## 2 Model

We consider an atomistic model in $`1+1`$ dimensions (substrate + height) which we treat as a cross section of the real $`2+1`$ case. An implicit assumption is that in the real $`2+1`$ model the monolayer islands have a compact rather than a fractal shape and the lattice misfit is one and the same in both orthogonal directions. Furthermore, we exclude from our considerations the possible interdiffusion and the subsequent gradient of strain as found recently by Kegel et al  in the case of InAs/GaAs quantum dots. The 3D islands are represented by linear chains of atoms stacked one upon the other as in the model proposed by Stoop and van der Merwe  and later by Ratsch and Zangwill , each upper chain being shorter than the lower one by one atom. In this sense the lateral size, and particularly the height of the islands are discrete parameters, whereas in most of the theoretical considerations they are taken as continuous variables .

In a previous paper , we used the method of computation proposed by Ratsch and Zangwill , which is based on the well-known model of Frenkel and Kontorova . The latter treats the overlayer as a linear chain of atoms subject to an external periodic potential exerted by a rigid substrate . Ratsch and Zangwill accepted that each layer (chain) presents a rigid sinusoidal potential to the chain of atoms on top of it. The potential trough separation of the lower chain is taken constant and equal to the average of all trough separations. As the strains of the bonds that are closer to the free ends are smaller, the average bond strain of each upper chain is closer to zero. In other words, the lattice misfit decreases from $`\epsilon _0`$ at the island’s base to zero at the apex. This method is, however, inadequate to describe properly a thickening overlayer because of one basic assumption, namely, the rigidity of each monolayer upon formation of the next one on top of it. This assumption rules out the relaxation and redistribution of the strains in the lower layers when upper layers are added. In particular, this method does not allow to compute the structure and energy of the interfacial boundary between the wetting layer and the 3D islands upon thickening of the latter.

For the above-mentioned reasons, in the present work we make use of a simple minimization procedure. The atoms interact through a Morse potential that can be easily generalized to vary its anharmonicity by adjusting two constants $`\mu `$ and $`\nu `$ ($`\mu >\nu `$) that govern separately the repulsive and the attractive branches, respectively \[44-46\],

$`V(x)=V_o\left[{\displaystyle \frac{\nu }{\mu \nu }}e^{\mu (xb)}{\displaystyle \frac{\mu }{\mu \nu }}e^{\nu (xb)}\right],`$ (1)
where $`b`$ is the equilibrium atom separation. For $`\mu =2\nu `$ the potential (1) turns into the familiar Morse potential, which has been used in the present work for the case $`\nu =6`$.

The pair potential designed by Tersoff for description of the properties of materials with directional covalent bonds like Si contains an additional parameter which accounts for the local atomic environments around the neighboring atoms . He showed that most of the properties of Si could be computed with an error smaller than 1%, compared with experimental data and ab initio calculations, by accounting only for the first neighbor interactions. For this reason, we occasionally consider only interactions in the first coordination sphere in order to mimic the directional bonds that are characteristic for the most semiconductor materials.

Our programs calculate the interaction energy of all the atoms as well as its gradient with respect to the atomic coordinates, i.e. the forces. Relaxation of the system is performed by allowing the atoms to displace in the direction of the gradient in an iterative procedure until the forces fall below some negligible cutoff value. The calculations were performed under the assumption that the substrate (the wetting layer) is rigid. This assumption is strictly valid in the beginning of the 2D-3D transformation when the 3D islands are still very thin .

Yu and Madhukar computed recently, by making use of the Stillinger-Weber interatomic potential  in a molecular dynamics study, the distribution of the strains and stresses in and around a 3D Ge island having a shape of a full pyramid with a length of the base edge 326 Å and a heigth of 23 monolayers . They found that the atoms in the middle of the first atomic plane of a coherent 3D Ge island are displaced upwards by 0.6Å that is approximately half of the interplanar spacing of Ge(001) (1.4Å), whereas the atoms at the island’s edges are displaced slightly downwards. The same holds for the vertical displacements of the atoms belonging to the uppermost Si plane. As the vertical displacements strongly influence the adhesion of the islands to the wetting layer we also performed preliminary calculations in which the uppermost three monolayers were allowed to relax. The results of these calculations demonstrated qualitatively the same behavior as in the case of a rigid substrate. For this reason, we present here only the results obtained under the assumption of the rigid substrate. Detailed systematic studies of the effect of the substrate relaxation will be published elsewhere.

## 3 Results

Fig. 2(a) shows the horizontal displacements of the atoms of the base chain of a coherently strained island, for a value of the misfit of 7%. The displacements are referenced to the sites the atoms would occupy if they belonged to the next complete monolayer, which would then be a part of the wetting layer. It can be seen that the end atoms are strongly displaced as in the model of Frenkel and Kontorova  and of Frank and van der Merwe . Increasing the island height leads to greater displacements of the end atoms. The reason is the effective increase of the strength of the lateral interatomic bonding in the overlayer with greater thickness as predicted by van der Merwe et al. . According to these authors an island with a bilayer height could be approximately simulated by a monolayer height island but with twice stronger lateral bonds. Fig. 2(b) shows the vertical displacements of the base atoms relative to the interplanar spacing between the monolayers belonging to the wetting layer. It is obvious that the vertical displacements are due to the climbing of the atoms on top of underlying atoms as a result of the horizontal displacements. The thicker the islands the greater are the horizontal displacements (for reasons discussed above) and in turn the vertical displacements. The results shown in Fig. 2 clearly demonstrate that the bonds that are close to the island’s edges are much less strained compared with these in the middle, in agreement with the results obtained by Ashu and Matthai  and Orr et al.  but contrary to the finding of Yu and Madhukar .

The interconnection between the vertical and the horizontal displacements is beautifully demonstrated in Fig. 3 where they are shown in an island containing two MDs. The horizontal displacements in this case are greatest in the cores of the MDs and so are the vertical displacements. This figure shows in fact the physical reason for the incomplete wetting in the classical SK mode. The adhesion is weaker owing to the introduction of MDs.

In order to illustrate the effect of the atom displacements on the adhesion of the separate atoms belonging to the island’s base chain, we plot their energy of interaction with the underlying wetting layer (Fig. 4) for coherently strained islands. As seen the atoms that are near to the chain ends (island’s edges) adhere much more weakly with the substrate. The influence of the potential anharmonicity is clearly demonstrated. Only one or two end atoms in the expanded chain adhere more weakly to the substrate whereas more than half of the atoms at both ends in the compressed chain are weakly bound. The figure demonstrates in fact the physical reason for the coherent SK mode which is often overlooked in theoretical models. Moreover, it is a clear evidence of why compressed rather than expanded overlayers exhibit greater tendency to coherent SK growth.

It follows from Fig. 2 that increasing the island thickness leads to weaker adhesion of the 3D islands to the wetting layer and, in turn, to the stabilization of the coherent 3D islands. This is clearly demonstrated in Fig. 5 which shows the dependence of the mean adhesion parameter $`\Phi `$ on the islands’ height for positive and negative values of the misfit. It is calculated as the average of the interaction energy between the base chain atoms and those of the wetting layer and is referenced to the corresponding value for a non-misfitting monolayer of the same size. It can be seen that $`\Phi `$ saturates beyond a thickness of about 5 monolayers as expected. Note that in this case the incomplete wetting ($`\Phi >`$ 0) is due solely to the misfit, the bonding in both phases $`A`$ (the wetting layer) and $`B`$ (the 3D islands) being nearly one and the same. What is more important is that the adhesion parameter in compressed islands is visibly larger than that in expanded islands which is due to the anharmonicity of the interatomic potential. This behavior clearly shows the greater tendency of the compressed overlayers to form coherent 3D islands. Another important feature that characterizes the mean adhesion parameter is its large absolute value. It is comparable with the values that lead to 3D islanding in VW mode of growth on chemically unlike surfaces .

Fig. 6 shows the mean adhesion parameter $`\Phi `$ as a function of the misfit both negative and positive for coherent islands as well as for islands containing one and two MDs. As discussed in the Introduction this is in fact the thermodynamic driving force for 3D islanding. Several interesting properties are observed. First, the mean adhesion parameter of compressed coherent islands is greater than that of expanded islands. This means that the thermodynamic driving force for coherent 3D islanding is greater in compressed rather than in expanded overlayers. In the absence of MDs the incomplete wetting is due to the displacements of the end atoms (see Fig. 2). In expanded overlayers the end atoms interact with their neighbors by the weaker attractive branches and vice versa. As a result $`\Phi _\mathit{0}^+>\Phi _\mathit{0}^{}`$.

On the other hand, the opposite is observed for dislocated islands. This is very easy to understand bearing in mind that in the classical (dislocated) SK mode $`E_{AB}E_{AA}E_{MD}`$, and $`\Phi E_{\mathrm{𝑀𝐷}}/E_{\mathrm{𝐴𝐴}}`$ , $`E_{MD}`$ being the energy per atom of the MDs. MDs have higher energy in expanded overlayers as they represent regions with higher density of atoms which repulse each other with the stronger repulsive branches of the potential. It is exactly the opposite in compressed overlayers, so that $`E_{MD}^+<E_{MD}^{}`$ and $`\Phi _{\mathrm{𝑀𝐷}}^+<\Phi _{\mathrm{𝑀𝐷}}^{}`$. This means that in the classical SK growth the thermodynamic driving force of formation of dislocated 3D islands is greater in expanded rather than in compressed overlayers.

Another property is that the adhesion parameter of islands containing two MDs appears as a continuation of that of the dislocation-free islands. This is also easily understandable having in mind the similarity of the model with that of Frank and van der Merwe . Dislocation-free solutions exist until the misfit reaches the so-called metastability limit at which the end atoms reach the crests between the next potential troughs and two dislocations (because of the symmetry of the model) are simultaneously introduced at both free ends. The energetic barrier for this process is equal to zero (for a review see Ref. ().

In order to answer the questions posed in the Introduction we compare the energies per atom of mono- and multilayer islands (frustums of pyramids) with different thickness varied by one monolayer. The pyramids are bounded with the steepest (60) sidewalls as they have the lowest energy in models in $`1+1`$ dimensions . As calculated, the energy represents a sum of the strain energy and the energy of the surfaces relative to the energy of the same number of atoms in the bulk crystal . Fig. 7(a) demonstrates the energies per atom vs the total number of atoms of monolayer and bilayer height islands at $`\epsilon _0=0.03`$. As seen the monolayer height islands are always stable against the bilayer islands. The latter means that the thermodynamics do not favor coherent 3D islanding. Monolayer height islands will grow and coalesce until they cover the whole surface. MDs will be then introduced to relieve the strain. Fig. 7(b) demonstrates the same dependence (including also thicker islands) but at larger value of the misfit $`\epsilon _0=0.07`$. This time the behavior is completely different. The monolayer islands are stable against the bilayer islands only upto a critical volume $`N_{12}`$, the bilayer islands are stable in turn against the trilayer islands upto a second critical volume $`N_{23}`$, etc. This behavior is precisely the same as in the case of VW growth where the interatomic forces (the wetting) predominate and the lattice misfit plays an additional role . The same result (not shown) has been obtained in the case of expanded overlayers ($`\epsilon _0<0`$) with the only exception that monolayer height islands are stable against multilayer islands upto much larger absolute values of the misfit.

The mono-bilayer transformation is the first step of the complete 2D-3D transformation. Studying the critical size $`N_{12}`$ as a function of the misfit (see Fig. 8) shows the existence of critical misfits beyond which the formation of multilayer islands can only take place. Below the critical misfit the monolayer height islands are stable irrespective of their size and the growth will continue in a layer-by-layer mode until MDs are introduced to relax the strain. The nearly twice larger absolute value of the negative critical misfit is obviously due to the anharmonicity of the atomic interactions. The weaker attractive interatomic forces lead to smaller displacements both lateral and vertical of the end atoms and in turn to stronger adhesion. The latter requires larger misfit in order for the 3D islanding to take place.

## 4 Discussion

The existence of critical misfit clearly shows that the origin of the 3D islanding in the coherent SK growth is the incomplete wetting which is due to the atomic displacements near the islands edges. As seen in Fig. 6 the mean adhesion parameter $`\Phi `$, or which is the same, the thermodynamic driving force $`\mathrm{\Delta }\mu `$ for coherent 3D islanding has practically the same values as that in the case of the classical (dislocated) SK mode at sufficiently large values of the misfit. Moreover, the comparison of Fig. 2 and Fig. 3 shows that at a given misfit and a given thickness, there is a critical lateral size (or a critical volume) beyond which MDs are spontaneously introduced to relieve the strain. It in fact determines the transition from the coherent to the classical (dislocated) SK growth which in the real case should be accompanied with the change of the shape. All the above leads to the conclusion that the physical reason for both the classical and coherent SK mode is one and the same.

The average adhesion depends strongly on the anharmonicity of the interatomic forces. Expanded islands adhere more strongly to the wetting layer and the critical misfit beyond which coherent 3D islanding is possible is much greater in absolute value compared with that in compressed overlayers. As a result coherent SK growth in expanded films could be expected at very (unrealistically) large absolute values of the negative misfit. The latter, however, depends on the materials parameters (degree of anharmonicity, strength of the chemical bonds, etc.) of the particular system and cannot be completely ruled out. Xie et al  studied the deposition of Si<sub>0.5</sub>Ge<sub>0.5</sub> films in the whole range of 2% tensile to 2% compressive misfit on relaxed buffer layers of Si<sub>x</sub>Ge<sub>1-x</sub> starting from $`x=0`$ (pure Ge) to $`x=1`$ (pure Si). They found that 3D islands are formed only under compressive misfit larger than 1.4%. Films under tensile misfit were thus stable against 3D islanding in excellent agreement with the predictions of our model.

The existence of a critical misfit for 2D-3D transformation to occur both in compressed and expanded overlayers has been noticed in practically all systems studied so far. Pinczolits et al.  have found that deposition of PbSe<sub>1-x</sub>Te<sub>x</sub> on PbTe(111) remains purely two dimensional when the misfit is less than 1.6$`\%`$ in absolute value (Se content $`<30\%`$). Leonard et al. have successfully grown quantum dots of In<sub>x</sub>Ga<sub>1-x</sub>As on GaAs(001) with $`x=0.5`$ ($`\epsilon _03.6\%`$) but 60Å thick 2D quantum wells at $`x=0.17`$ ($`\epsilon _01.2\%`$). Walther et al.  found that the critical In content is approximately $`x=0.25`$, or $`\epsilon _01.8\%`$. As commented before, a critical misfit of 1.4% has been found by Xie et al upon deposition of Si<sub>0.5</sub>Ge<sub>0.5</sub> films on relaxed buffer layers of Si<sub>x</sub>Ge<sub>1-x</sub> with varying composition .

A rearrangement of monolayer height (2D) islands into multilayer (3D) islands has been reported by Moison et al.  who established that the 3D islands of InAs begin to form on GaAs at a coverage of about 1.75 ML but then the latter suddenly decreases to 1.2 ML. This decrease of the coverage in the second monolayer could be interpreted as a rearrangement of an amount of nearly half a monolayer into 3D islands. The same phenomenon has been noticed by Shklyaev, Shibata and Ichikawa in the case of Ge/Si(111) . Voigtländer and Zinner noted that Ge 3D islands in Ge/Si(111) epitaxy have been observed at the same locations where 2D islands locally exceeded the critical wetting layer thickness of 2 bilayers . Bhatti et al.  and Polimeni et al.  also reported the coexistence of large pyramids and small flat islands. These observations show that the 2D islands really appear as precursors for the 3D islands.

The question of the existence and particularly the stability of the intermediate states is more difficult to answer. Rudra et al. measured photoluminescence (PL) spectra of InAs layers deposited on InP(001) at two different temperatures (490 and 525C) and buried in the same material . When the layers were grown at 490C and the capping layer was deposited immediately after the deposition of the InAs the spectrum consisted of a single line. If the InAs layer was annealed for 10 s before capping with InP the spectrum consisted of 8 lines. At 525C 3 lines were observed already in absence of annealing. The above observations could be explained by formation and coexistence of islands with different thickness varying by one monolayer. Colocci et al.  performed PL studies of InAs deposits on GaAs(001) with thickness slightly varying around the critical thickness of 1.6 monolayers for the onset of the 3D islanding. They observed an increasing number of luminescence lines with increasing film thickness. These lines were attributed to families of 3D islands with similar shape but with heights differing by one monolayer. Flat platelets, 2 - 6 monolayers high, have been observed during the growth of GaN/AlN heterostructures .

Although the above results seem to be in an excellent qualitative agreement with the theoretical predictions of the model, the thermodynamic stability of islands with quantized height of one monolayer, and the existence of a critical misfit is still debated . The reason of the discrepancy of our results with those of Duport et al.  most probably stems from the implicit assumption, made by the above authors, that the widths of the lower, $`R`$, and the upper, $`R^{}`$, bases, and particularly the height $`h`$, of the crystal having a shape of a frustum of a pyramid, they consider, are continuous variables. This is correct if the crystals are sufficiently large. However, the continuum approximation is not acceptable in the beginning of the 2D-3D transformation when the islands are still very small (and thin). It is also not applicable in the limit $`hR`$ for the same reasons. The question of existence of a critical misfit follows logically if we accept that the intermediate states with heights differing by one monolayers exist and are thermodynamically stable in consecutive intervals of the volume.

In conclusion, accounting for the anharmonicity of the real interatomic potentials in a model in $`1+1`$ dimensions, we have shown that coherent 3D islands can be formed on the wetting layer in the SK mode predominantly in compressed overlayers. Coherent 3D islanding in expanded overlayers could be expected as an exception rather than as a rule. The thermodynamic driving force for 3D islanding on the wetting layer of the same material is identified as the weaker adhesion of the atoms near the islands edges. This should also facilitate the 2D-3D transformation. Overcoming a critical lattice misfit appears as a necessary condition for the formation of coherent 3D islands. Monolayer height islands of a critical size appear as necessary precursors of the 3D islands. The 2D-3D transformation from monolayer islands to 3D pyramids takes place through a series of intermediate states with heights increasing by one monolayer. The intermediate states are thermodynamically stable in consecutive intervals of the volume.

###### Acknowledgements.

The authors are indebted to the Instituto Universitario de Ciencia de Materiales “Nicolás Cabrera” for granting research visits which enabled scientific collaboration. This work was supported by the Spanish CICyT through project Nr. MAT98-0965-C04-02.