text:
# Study of Variable Stars in the MOA Database11footnote 1: Long-Period Red Variables in the Large Magellanic Cloud ⋆ Partly based on DENIS data obtained at the European Southern Observatory

## 1 Introduction

There are many types of variable stars whose mechanism for variability (mainly pulsation, rotation, eruptive processes, and geometric effects) is well studied and well understood observationally and theoretically. Among them, the pulsating variable stars have been investigated vigorously because they yield information on their stellar mass and their interior physical properties. They are a powerful tool for measuring distances in the Galaxy and extra-galactic systems (1) because they are luminous and (2) because they have been shown empirically to obey a period-luminosity (PL) relation.

The PL relation for classical Cepheids was discovered by Leavitt (1908) and used to measure the distance of galaxies. However, there is a discrepancy between the distance to the LMC derived from the PL relations for Cepheids and that derived for RR Lyraes (Walker 1992). This discrepancy may be due to uncertainties in the absolute calibration and a sensitivity to metallicity effects (Caputo 1997). It is important to establish PL relations of any new objects which are independent of these variable stars.

A PL relation for long period red variables (Miras and semi-regular variables (hereafter SRs)) will be a useful distance indicator because these stars are numerous and have high luminosities. However, they are the least studied star type because of their complex pulsation mechanism. Difficulties arise from uncertainties in our understanding of the convective energy transport and the opacity resulting from the many molecular absorption lines.

We present in this paper the results of our study of LMC Miras and SRs using the MOA database. This database was generated from the MOA project microlensing survey.

A PL relation for Miras of the LMC was presented by Glass & Lloyd-Evans (1981) and Glass & Feast (1982). They observed 11 Miras in the infrared bands $`J,H,K`$. However, there was a large amount of scatter in the plot ($`\sigma 0.25`$ mag) because their observation time was limited and it did not cover the whole phase of variations.

To reduce the scattering in the plot due to poor phase coverage, extensive infrared observations have been made (Feast et al. 1989; Hughes & Wood 1990), and a narrow PL relation was obtained ($`\sigma 0.13`$ mag). Feast (1996) showed that the Miras in the LMC follow a narrow PL relation in both absolute bolometric magnitude, $`M_{bol}`$, and absolute $`K`$-magnitude, $`M_K`$, and he concluded that all Miras pulsate in the same mode.

Wood & Sebo (1996) found two sequences in a $`K`$ band PL diagram for the small amplitude SRs in the LMC. One was the well-known Mira sequence, and the other (eight stars) was located at the position which is parallel to the Mira sequence and their periods were approximately half the periods of Miras on the PL diagram. This discovery encouraged a hypothesis that the pulsation mode of these stars in the secondary sequence must be in a higher mode. Using the Hipparcos database of nearby stars, Bedding & Zijlstra (1998) discovered the SRs which are located in the secondary sequence of Wood & Sebo (1996). Bedding and Zijlstra interpreted the sequence of SRs as evolutionary (see Whitelock (1986)).

Feast (1999) suggested that both Miras and SRs pulsate in the same mode, possibly both in the first overtone. The arguments were based on the angular diameter measurements of nearby Miras and the PL relationship of the Miras and SRs in globular clusters. He agreed with Bedding and Zijlstra (1998) that the sequence which is located on the shorter period side of the Mira sequence in the PL plane is an evolutionary track ending at the Mira sequence. Feast also suggested that the Mira sequence would be a collection of variable stars with different chemical abundances and/or with different stellar mass. In this case, it was not necessary that Miras and SRs pulsate in different modes.

As described above, there are two theories on the short-period stars in the PL diagram for long-period red variable stars. Considering that the pulsation features of the Miras and SRs depend on several factors, other physical properties must also be considered. The study of stellar colour will be interesting because this gives us the temperature of stellar surface, stellar radius and mass which are the important factors of the pulsation period.

Recently, the continuous observation of a great number of stars using CCD cameras which cover a large image area became possible. The recent report of Wood (2000) based on the $`K`$ band photometry and the period derived from the MACHO database indicated plausible sequences on the PL diagram. This diagram shows at least five sequences falling parallel to the Mira sequence. According to Wood, the three sequences containing the Mira sequence were explained as being different pulsation modes. However, more information is required to settle this issue. The features of the other two sequences are unknown, and it is expected that further independent observations and detailed investigation will reveal the nature of these stars.

We worked on the problem using the MOA database which contained data for a large number of stars in the LMC. The period, magnitude colour, amplitude and periodicity for each star was determined.

We found an independent group of short period stars in a sample of 146 stars. We investigated the properties of these stars through their colour, amplitude and periodicity, and discuss the possibility of that they are higher mode pulsators.

Following an introduction to the MOA project, we describe our observation and photometry systems in Sections 2 and 3. To obtain the light curves of variable stars, a careful calibration of each image frame is required. This procedure is described in Section 4. In Sections 5 and 6, the data reduction and period analysis methods are described. In Section 7, the colour-magnitude diagram (CMD) of the LMC stars observed using our system is presented and the period-magnitude relation will be given in Section 8. The period-colour relation is presented and discussed in Section 9. The result of the conversion of MOA magnitudes to $`K`$ band values will be given in Section 10, along with a $`K,\mathrm{log}P`$ diagram. In Section 11 we discuss the possibility of the existence of multiple sequences in the ($`\mathrm{log}P,K`$) plane with the amplitude distribution. In Section 12 we summarise our results.

## 2 The MOA Project Observations

### 2.1 Observational Instruments

The MOA project commenced observations in May 1996 with the aim of observing microlensing events towards the LMC, SMC and Galactic Bulge. The observational program was undertaken at the University of Canterbury’s Mt. John Observatory ($`44^{}S,1030\mathrm{m}`$ above sea level) in the centre of the South Island of New Zealand.

The MOA project uses a 61cm Ritchey-Chrétien Cassegrain telescope with modified optics to give a wide FOV at f/6.25.

For the first set of observations in 1996 (Series 1), we used the original f/13.5 telescope optics and monitored $`3\times 10^5`$ stars in three LMC fields with the mosaic CCD camera, MOA-cam1. MOA-cam1 was constructed with nine TI TC215 CCDs with $`1000\times 1018`$ pixels whose pixel size was $`12\mu `$m $`\times 12\mu `$m. From January 1997 to July 1998, in Series 2, the optics were altered to f/6.25. MOA-cam1 covered 3 fields in the LMC bar ($`1\times 10^6`$ stars) and 2 fields in the SMC ($`4\times 10^5`$ stars). The current set of observations, Series 3, started in August 1998. For this series, a new camera, MOA-cam2 is in use. MOA-cam2 consists of three abutted SITe $`2047\times 4095`$ pixel thinned, back-illuminated CCDs with $`15\mu `$m $`\times 15\mu `$m pixels. The pixel size is $`0.81^{\prime \prime }`$, giving a field of view of $`0.5^{}\times 0.9^{}`$ for each chip. The details of MOA-cam2 are in Yanagisawa et al. (2000). In this series we are observing 16 LMC fields ($`4.4\times 10^6`$ stars), 8 SMC fields ($`9.3\times 10^5`$ stars) and 14 Galactic Bulge fields.

We used two broad-band filters to maximise photon collection. The MOA blue filter ($`MOA_B`$) covers $`395620n`$m ($`\lambda _{\mathrm{eff}}500`$ nm) and MOA red filter ($`MOA_R`$) covers $`6201050`$ nm ($`\lambda _{\mathrm{eff}}700`$ nm). The spectral transmission curves for each filter are shown in Fig. 1.

The MOA instrumental photometry has been transformed to the standard Johnson $`V`$ and Kron-Cousins $`R`$ ($`R_{KC}`$) by comparisons with about one hundred stars imaged by Hubble Space Telescope (HST). We describe these details in Section 4.

The variable star study based on the MOA database was reported in Hearnshaw, Bond, Rattenbury, et al. (2000), Noda, Takeuti, Bond, et al. (2000), and Takeuti, Noda, Bond, et al. (2000).

### 2.2 OBSERVATION DETAILS

Regular 300 second exposures were made towards the Large Magellanic Cloud and the Small Magellanic Cloud. 180 second exposures were made for the Galactic Bulge (GB) fields. Under the current system, it takes 70 seconds to read out the 3 CCD images for each exposure. Dark frames for each exposure were prepared every photometric night and dome flats or dark-sky flats were also frequently taken.

Observational data were archived on DLT tapes and all the data were analysed in both Japan and New Zealand. Series 1, 2 and 3 images were analysed using DoPHOT photometry software (Schechter et al. 1993). The size of the Series 2 database created by the subsequent photometry is $``$ 40 Gbytes.

A new online analysis system based on image subtraction has been developed in New Zealand. In the year 2000 this system was able to find several microlensing events in real time and issue alerts to the astronomical community via the WWW222http://www.phys.canterbury.ac.nz/$``$physib/alert/alert.html. We are re-analyzing some parts of our data with this method.

## 3 DATA ANALYSIS

Before undertaking data analysis using the DoPHOT program, all image files were prepared with dark and the flat field frames. DoPHOT was used in a fixed position warm start mode, using the best photometric images as templates. It takes about 50 minutes to analyze one LMC frame with our computer. There were more than 400 measurements in the Series 2 (1.6 year) and Series 3 (1.4 year) database. We obtain only relative magnitudes from the DoPHOT analysis. We also use non-standard filters, thus a transformation of our instrumental magnitude to the standard $`V`$ and $`R_{KC}`$ is required. This transformation is described in the next section.

## 4 CALIBRATION FOR MOA PHOTOMETRY

The transformation of the MOA photometry to the standard system was carried out using Hubble Space Telescope (HST) image data. More details are given by Kato (2000). We defined the transform relation for MOA-cam1 and MOA-cam2 separately because their CCD quantum efficiencies were quite different (the quantum efficiency of MOA-cam2 was twice as high as that for MOA-cam1). Although the filter system is exactly the same, this will cause serious differences in the relative magnitudes returned by DoPHOT. The DoPHOT magnitudes were calculated using photon quantities only, so a difference in CCD quantum efficiency will affect the DoPHOT output directly. To get standard magnitudes for MOA-cam2, we performed the procedure described below.

Photons coming from a star show a complicated scattering pattern around a central position caused by atmospheric turbulence. A point-spread function (PSF) models the scatter. The DoPHOT code estimated each magnitude using an elliptical Gaussian with seven free parameters as a model for a stellar PSF. A common PSF was adopted for all stars (bright stars and faint stars) in a frame, however it was clear that there were significant differences between the profiles for bright stars and that for faint stars. The former were weighted by the photon statistics of the star itself, and the latter were weighted by the photon statistics of the sky, so the latter should show larger scattering. The detailed study of the PSF indicated that the function often differs from the elliptical Gaussian (Alard & Lupton 1998). The error due to this discrepancy is approximately $`\pm `$0.05 mag for the bright stars and $`\pm `$0.09 mag for the faint stars. It was not a critical problem if we define a mean magnitude with more than 100 measurement points to study stellar intensity, because the scatter due to PSF model inadequacy usually occurs at random.

As we were using large CCD chips, we have to determine the position dependent error for each CCD chip. In order to investigate this, we observed the globular cluster NGC3201 at various positions on each chip, and compared stellar magnitudes with published magnitudes obtained using $`B,V,R,I`$ filters by Alcaino et al. (1989). We found that the error was smaller ($`\pm 0.015`$ mag) than the DoPHOT PSF photometry error ($`\pm 0.09`$ mag). This was a systematic error however, and did not disappear even if we deal with a mean magnitude.

In the third step, also using the NGC3201 data, we computed the chip to chip differences for the magnitude zero point:

$`MOA_{B(chip3)}=MOA_{B(chip1)}0.23`$ (1)
$`MOA_{R(chip3)}=MOA_{R(chip1)}0.31`$ (2)
$`MOA_{B(chip3)}=MOA_{B(chip2)}0.34`$ (3)
$`MOA_{R(chip3)}=MOA_{R(chip2)}0.41`$ (4)
We used the chip3 scale as our standard as it gave the most linear response.

For the final step, we converted our $`MOA_B`$ and $`MOA_R`$ magnitudes to standard Johnson $`V`$ and Kron-Cousins $`R`$ ($`R_{KC}`$) respectively. The brightness of stars was acquired from an ADU (Analog Digital Unit, 1 ADU corresponded to 1.3 electrons for MOA-cam2) value which corresponded to the number of photons measured with the assumed PSF via the DoPHOT photometry. We define the zero-point of magnitude for MOA as 1 ADU. Therefore all of the obtained magnitudes are negative, and we have to correct the zero-point of our magnitude in order to obtain actual magnitudes. We analyzed the HST data of the LMC obtained with the Wide Field Planetary Camera2 (WFPC2) to correct for atmospheric seeing. Since the filters for WFPC2 had been transformed to the Johnson $`V`$ and $`R_{KC}`$ by Holtzman et al. (1995), we were able to compare the converted HST magnitudes directly with those of our one hundred stars.

Finally, we found a relation to convert the $`V_m`$ and $`R_m`$ for MOA-cam2;

$$V_m=MOA_B0.16(\pm 0.01)(MOA_BMOA_R)+27.41(\pm 0.01)$$
(5)
$$R_m=MOA_R+0.29(\pm 0.01)(MOA_BMOA_R)+27.24(\pm 0.01)$$
(6)
$$Colour(V_mR_m)=0.53Colour(MOA_BMOA_R)+0.18$$
(7)
For MOA-cam1;

$$V_m=MOA_B0.10(\pm 0.01)(MOA_BMOA_R)+25.55(\pm 0.01)$$
(8)
$$R_m=MOA_B+0.31(\pm 0.02)(MOA_BMOA_R)+25.45(\pm 0.1)$$
(9)
$$Colour(V_mR_m)=0.59Colour(MOA_BMOA_R)+0.10$$
(10)
Hereafter, we deal with $`V_m`$ and $`R_m`$ as the MOA standard magnitudes. In Fig. 2, we show the Colour Magnitude Diagram (CMD) obtained by MOA-cam2 for our observation fields in the LMC bar. The limiting magnitude was around 20.5 mag for MOA-cam2, and 19 mag for MOA-cam1. Since our passbands were broad and covered frequencies from violet to near infra-red bands, the colour is not sensitive for blue stars.

To minimize the effects of the non-uniformity of the sky and system on the raw magnitudes (i.e. DoPHOT output magnitudes) we:

1. Divide the CCD chips into small regions (200 $`\times `$ 200 pixels)
2. Select eighty bright stars which had more than 100 measurement points in each small region
3. Regard stars with a small scatter as ‘non-variable’ stars, i.e., where the number of observed points more than 4 sigma from the average magnitude was less than 15
4. Calculate magnitudes for all stars by comparison with the ‘non-variable’ stars. That is, the MOA standard magnitudes were determined relative to the magnitudes of the near ‘non-variable’ stars in each small region.

Through the revision of our calibration procedure, we obtained progressive and interesting diagrams compared with our previous reports (Hearnshaw, Bond, Rattenbury, et al. 2000; Noda, Takeuti, Bond, et al. 2000; Takeuti, Noda, Bond, et al. 2000). Details are given in later sections.

## 5 Data Reduction

We selected long-period red variables from the Series 2 database using the following criteria:

1. Each star should be ‘photometric’ and have a ‘single star shape’ as judged by DoPHOT.
2. Light curves should be sufficiently long ($``$ 100 points).
3. Light curves of each star must have large variability, i.e., the reduced $`\chi ^2`$ for a fit to mean brightness should be larger than $`\chi _{Vm}^2\chi _0^2`$ and $`\chi _{Rm}^2\chi _0^2`$. Here $`\chi _0^2`$ corresponded to the upper $`10\%`$ limit of the reduced $`\chi ^2`$ distribution. For the MOA database $`\chi _0^2=4.3`$.

To obtain large amplitude variables, we applied the additional criteria:
4. Low amplitude stars at an early stage of the study are rejected by removing those stars which had a magnitude difference between the fourth highest and fourth lowest magnitude values of less than 1.3 magnitudes. This difference was larger than the amplitude $`\delta R_m`$ determined in Section 6. This selection was done for the results of MOA-cam1 system without using the colour dependence in equations (8) - (10). This sample of large amplitude stars might be affected when we use the correct colour relation.
5. We select redder stars for large amplitude; $`(V_mR_m)0.4`$ for the MOA-cam1 system.

We obtained 1,067 light curves for the candidates of variable stars with large amplitude after applying criteria (iv) and (v) on the Series 2 database. Light curves from Series 2 and Series 3 were cross referenced to stars from their respective templates. We rejected unreliable light curves, and finally obtained 908 light curves. In Fig. 3 we show an example of a light curve.

We also found 485 small amplitude variable star candidates by applying a further criterion after numbers (i) to (iii) above. We analyse the large amplitude variables in this work. We will present the small amplitude variable star results in a future paper.

## 6 Period Analysis (PDMM)

As the first stage in the period analysis, we used the Phase Dispersion Minimization method for MOA (PDMM) which is a PDM method (Stellingwerf 1978) modified to operate on large photometry datasets. The PDMM algorithm finds a likely light curve period together with a reliability parameter, $`\theta `$, for a large number of light curves automatically. $`\theta =0`$ indicates a complete cyclic change and $`\theta =1`$ indicates no periodicity. The amplitude, the difference between the magnitudes of the brightest and the faintest bins, was also calculated. We converted the magnitudes of each bin to intensities and calculated the mean intensity. We denote the mean magnitude values as $`V_m`$ and $`R_m`$. We executed this code using both colours ($`V_m,R_m`$) for the 908 selected light curves.

Since our dataset included scattered data points, the PDMM algorithm sometimes gave integral multiples or sub-multiples of the period, or in the worst case an incorrect result. To avoid simple errors in the periods obtained by the automatic analysis, we checked all code output carefully, and ran the PDM code manually to determine the most likely periods. We accepted a period which had good agreement between both colours. In Fig. 4, we present an example of PDMM figures. The upper diagram shows the reliability of the testing period. A small value of $`\theta `$ indicates a likely value for the period. The bottom diagram shows the light curve folded with the most likely period.

The total dataset (Series 2 to Series 3) extends over about 3 years ($``$1,060 days), so periods less than $``$400 days should be well determined because these have more than 2.5 cycles of variation within the observation period. We classified 908 candidate light curves by their features of periodicity. Among 908 stars, 329 stars show clear variability and of these, 259 stars had periods which were relatively well-determined. We rejected light curves like eclipsing binaries, then chose light curves whose periods were longer than $`logP1.5`$ (30 days). We also required $`\mathrm{log}P2.6`$ (400 days) reflecting the dataset time span. Finally, 146 light curves were left. The reliability parameter $`\theta `$ for these curves is less than 0.67 (the mean $``$ 0.38) with 15 folding bins. We show the distribution of $`\theta `$ in Fig. 5.

Other candidates, the stars with poor determined period, short period or small amplitudes, will be the subject of future work. Here we give observational values for 146 stars in Table 1 (ID, $`\alpha `$, $`\delta `$, $`R_m`$, $`V_m`$, $`\delta R_m`$, $`\delta V_m`$, Period, $`\theta `$). ID is the MOA database star number, ($`\alpha ,\delta `$) is the right ascension and declination respectively. Comparing the coordinates of 629 stars in Tables 2 $``$ 6 of Hughes & Wood (1990) with our sample, we identified at least 34 stars as a control. The ID numbers of their tables are indicated in the Remarks column. Spectral types given in their tables are also shown. The 34 stars consisted of 21 oxygen-rich stars and 13 carbon-rich stars. In Fig. 6, we show the relation between the period determined via the PDM method and the period tabulated in Hughes & Wood (1990) for 30 stars whose periods are well-determined by us. It is clear that the periods are coincident and the differences between MOA and Hughes & Wood (1990) is less than 2 %. We used these 30 stars for the conversion to the $`K`$ magnitude system as described in the following section. Here, we note that Hughes & Wood (1990) tabulated two periods, 170 days and 319 days, for the star nlmc3-288533. We adopted 319 days as their period because we obtained 307.6 days for the period of this star via our analysis using the PDM method.

## 7 Colour Magnitude Diagram (CMD)

Fig. 2 shows the distribution of our sample of 146 variable stars on a colour-magnitude diagram (CMD) with 9,050 randomly chosen stars for reference. In this figure, the 146 stars were located at the tip of the Giant Branch. According to the stellar evolution theory, stars with initial mass of about $`18M_{}`$ evolve into Asymptotic Giant Branch (AGB) stars with carbon and oxygen cores. Although we only required a large variation in each light curve and restricted the colour condition and period in our selection criteria for variable stars (Section 6), these stars were located in one particular region in the CMD. Consequently we presumed that these stars form a group which were in a common evolutionary stage and pulsate by a common mechanism. However there was an extremely bright star (indicated by a $``$ mark) in this CMD, nlmc3-214266, which we have mentioned below. Also we will describe the seven stars with + marks in Section 11. nlmc3-214266 and these seven stars are denoted by (1) or (a) $``$ (g) respectively in Table 1.

## 8 Period-Luminosity Relation

We show a plot of $`V_m`$ magnitude versus period in Fig. 7. Open circles indicate stars whose $`R_m`$ amplitudes ($`\delta R_m`$) are smaller than 2.2 mag. Filled circles indicate the stars whose amplitude are larger than 2.2 mag. The line indicates the regression fit for the stars with $`\mathrm{log}P2.2`$;

$$V_m=2.01(\pm 0.43)\mathrm{log}P+12.10(\pm 1.05)$$
(11)
The trend to be fainter in $`V_m`$ has been known in previous studies. The star, which was extremely bright in the CMD, nlmc3-214266 ($`V_m=14.69`$, $`\mathrm{log}P=2.05`$, $``$ mark), was also much brighter than the other stars with similar periods in this diagram. We also show the ($`\mathrm{log}P,R_m`$) plot in Fig. 8. The slope of this distribution for the stars with $`\mathrm{log}P2.2`$ was less than that in ($`V_m,\mathrm{log}P`$) plane but still had a downward tendency. The best fit line for these stars is;

$$R_m=1.13(\pm 0.33)\mathrm{log}P+12.84(\pm 0.80)$$
(12)
The star, nlmc3-214266, was again quite bright ($`R_m=13.80`$) and was separate from the other samples. The scatter of $`V_m`$ magnitude ($`\sigma 0.52`$) was larger than that of $`R_m`$ ($`\sigma 0.39`$), which was possibly due to the difference between oxygen-rich stars and carbon-rich stars. We will discuss this in Section 10.

## 9 Period-Colour Relation

We show the period-colour ($`\mathrm{log}P,V_mR_m`$) plot in Fig. 9. The upward line indicates the regression line for the stars with $`\mathrm{log}P2.2`$;

$$V_mR_m=0.88(\pm 0.16)\mathrm{log}P0.74(\pm 0.38)$$
(13)
In this figure, the colour of most stars with longer periods ($`\mathrm{log}P2.2`$) were much redder than that of the short period stars, and the relationship was more constrained ($`\delta 0.19`$) compared with that shown in Fig. 7 and Fig. 8.

### 9.1 nlmc3-214266

In Fig. 9, star nlmc3-214266 is positioned in the main group. Comparing with the other stars of similar period (open and filled circles), the effective temperature of this star should be similar since the colour ($`V_mR_m`$) was almost equivalent. The brightness of nlmc3-214266 was high in both $`V_m`$ and $`R_m`$, therefore its radius must have been large compared with the stars in the main group. The star was almost certainly HV 1011 which was a suspected Cepheid (Kurochkin 1992). Kurochkin was unable to obtain a period for this star and our results do not conclusively suggest that this star is a classical Cepheid. When the pulsation period is the same, a larger star must have a higher mass considering the relation between the period and the mean density of stars. Therefore we assumed that nlmc3-214266 was a massive star evolved from the main sequence, and not an AGB star. The light curve is presented in Fig. 10. The star nlmc3-214266 was a luminous semi-regular variable. In the Galaxy, V810 Centauri is one of the variable stars which is more luminous and has a higher temperature than the AGB in the CMD. Kienzle, Burki, Burnet, and Meynet (1998) estimated that the spectral type of this star was F8 Ia, and the period was approximately 156 days, $`M_{bol}=8.5`$, and the initial mass was $``$ 25 $`M_{}`$. Such a heavy star will become a supernova rather than a AGB star (e.g. Becker, 1998).

### 9.2 Short Period Red Stars

In Fig. 9, several stars whose periods were $`\mathrm{log}(P/d)<2.15`$ were distributed towards the red region and were separated from the main group of stars. If the colours are different between stars with a similar period and luminosity, their mass or pulsation property ($`Q=P\sqrt{\overline{\rho }}`$) should be different. It is quite important to investigate their luminosity because the relation between luminosity and colour gives information about their radius, i.e., mass and pulsation parameter ($`Q`$) of each star. Although our $`V_m`$ and $`R_m`$ are trustworthy they cover too small a wavelength range to enable the energy distribution of such red stars to be discussed directly. Therefore we study the relationship between our magnitudes and the near-IR luminosity in Section 10.

### 9.3 nlmc1-239264

We note that there were no data points for the minima of the light curve for nlmc1-239264. Therefore the mean magnitudes, $`V_m`$ and $`R_m`$, for this star were calculated by using only the bright phases. This gives just the brighter limit of the mean magnitude for this star. Regarding the period, we assumed that the determined period for the star ($`P=129.5`$ d, or $`\mathrm{log}P=2.11`$) was accurate because there were three maxima in the light curve. On Figs. 7 and 8, the star was located near the extension of the least squares fitted line for stars with $`\mathrm{log}P2.2`$. When we used a dataset covering complete cycles, the stars were located at fainter positions. We used $`V_mR_m=1.44`$ as the blue or lower limit of the colour average in Fig. 9. The star was located about half a magnitude redder than the least squares fitted line in Fig. 9. Considering the incompleteness of our dataset, we may assume that the star nlmc1-239264 should be located at a fainter position on Figs. 7 and 8, and at a redder or higher position on Fig. 9.

### 9.4 nlmc3-306907

The star, nlmc3-306907, is the remarkable star because of its unique colour ($`V_mR_m`$) variation, i.e., the colour was bluer when $`V_m`$ and $`R_m`$ were fainter. The light variation is shown in Fig. 11. Among 146 samples, nlmc2-151839 and nlmc3-191203 also showed the same property, though to a lesser extent. The colour variations of all the other stars were fainter when the $`V_m`$ and $`R_m`$ were brighter, which is a common feature of radial pulsators. Such an extraordinary property might indicate the development of a particular molecular band caused by a different chemical composition in the stellar surface.

## 10 PERIOD TO $`K`$-MAGNITUDE RELATION

It is well known that Miras fall on a narrow PL sequence on the ($`\mathrm{log}P,K`$) plane (Feast et al. 1989; Hughes & Wood 1990). Wood & Sebo (1996) discovered a second sequence located on the short period side of Miras, while Bedding & Zijlstra (1998) found strong evidence that Semi-Regulars (SRs) near the Sun’s position in the ($`\mathrm{log}P,K`$) plane lie in a region which coincides with the second sequence found by Wood & Sebo (1996). Also the existence of at least two (or four) sequences other than Miras was indicated from the analysis of MACHO database (Wood 2000).

The $`K`$-band ($`2.2\mu `$m) corresponds closely to the maximum radiation flux of long-period red variables, so that the discrepancy between $`K`$-magnitudes and bolometric magnitudes should be relatively small. Because previous studies have been done with $`K`$ magnitudes, we have to convert our observational magnitude $`V_m,R_m`$ into $`K`$-magnitudes in order to compare our results with the established PL relation.

First, to define our conversion scale, we identified 34 stars which were the same as those published by Hughes & Wood (1990) from our 146 sampled long-period red variable stars (see Section 6). The 34 identified stars consisted of 21 oxygen-rich (M-type) stars and 13 carbon-rich (C-type) stars. Furthermore, we selected 19 oxygen-rich stars whose periods were well-determined in our analysis, and used them to calibrate our $`V_m,R_m`$ magnitudes to $`K`$ values;

$$R_mK_{shv}=2.10(V_mR_m)+1.22$$
(14)
$`K_{shv}`$ was the $`K`$-magnitude tabulated in Hughes & Wood (1990) ($`shv`$ is the prefix of the star number). Hereafter we write $`K_m`$ in place of $`K_{shv}`$ as our custom $`K`$ magnitude. We discuss here a few reasons for the scatter ($`\sigma 0.15`$ mag) in this relation. Firstly, and possibly most importantly, the $`K`$-magnitudes obtained by Hughes & Wood (1990) were determined with only a single observation, whilst we were dealing with the average magnitudes of at least one hundred measured points. Although amplitudes of long-period red variable stars are large ($`0.6`$ mag) in the $`K`$-band, a single epoch observation provides the magnitudes on the various phase of light variations, so it should have caused extra scattering in our $`K_m`$ conversion. Secondly, our transformed magnitudes ($`V_m`$ and $`R_m`$) which were obtained with our broad band filters do not correspond to a standard V or $`R_{KC}`$ (also $`K`$) exactly. Finally, as mentioned in Section 4, there was an uncertainty of approximately $`\pm 0.015`$ mag dependent on the image position on each CCD chip. The result from comparing the $`K`$ magnitudes of Hughes & Wood (1990) and the average MOA magnitudes had a similar scatter to that due to the intrinsic variability of the stars. We also calculated the probable $`K`$ magnitudes by extrapolating the observed light curves backwards to the time of observation of Hughes & Wood (1990). The scatter found in this case was larger and most probably due to the imperfect periodicity of Mira variables. Since the resulting relation was nearly the same as equation (14) that expression was used in the following study.

Now we estimated the coincidence of our $`K_m`$ conversion with formula (14). In Fig. 12, we show a ($`V_mR_m,R_mK_{shv}`$) plot for 19 oxygen-rich stars (open circles) and 11 carbon-rich stars (filled circles) which were common in Hughes & Wood (1990) and our catalogue, and show the well-determined periods by our analysis. 19 oxygen-rich stars were used to derive formula (14) which is indicated by the solid line. It is obvious that the conversion formula (14) fits the oxygen-rich Miras (O-Miras) ($`\sigma 0.15`$ mag) well but does not fit the carbon-rich stars ($`\sigma 0.34`$ mag). The dashed line indicates another conversion formula using 11 carbon-rich stars:

$$R_mK_{shv}=3.16(V_mR_m)+0.22$$
(15)
We also show the ($`\mathrm{log}P,K_m`$) diagram for these 19 oxygen-rich stars and 11 carbon-rich stars in Fig. 13. The two lines in this figure were the sequences of O-Miras from previous studies (Feast et al. 1989; Hughes & Wood 1990). Our converted magnitudes $`K_m`$ for 19 oxygen-rich stars agreed well with these classical relations ($`\sigma 0.18`$ mag). In contrast, those for 11 carbon-rich stars appear below the lines ($`\sigma 0.25`$ mag). According to the results from previous observations, the PL relation for O-Miras (oxygen-rich Miras) and C-Miras (carbon-rich Miras) are similar in the $`K`$-band (Feast et al. 1989). In Fig. 12 and 13, we conclude that our conversion formula (14) makes a good enough fit to the oxygen-rich stars, but does not fit to the carbon-rich stars.

We consider here why such a disagreement arises for the carbon-rich stars. As is well known, there are roughly two kinds of spectral types, M and C, for long-period variable stars depending on the relative abundance of oxygen to carbon in their atmospheres. This causes a remarkable difference in colour tendency (Smak 1964). The B-band ($`V_m`$-magnitude in this paper) flux of oxygen-rich (M-type) stars was smaller than that for carbon-rich (C-type) stars. This was due to the existence of some strong TiO bands resulting from their rich abundance of oxygen. Therefore carbon-rich stars are bluer than oxygen-rich stars at the same effective temperature, which means that the value of ($`R_mK_{shv}`$) for carbon-rich stars is greater than that for oxygen-rich stars even if the colour index ($`V_mR_m`$) is the same (see Fig. 12). This explains why the $`K`$-magnitudes for carbon-rich stars were less luminous in Fig. 13, when we derived the conversion scale, $`K_m`$, using only oxygen-rich stars.

The above results were derived from the data by Hughes & Wood (1990). Within these stars (21 O-Miras and 13 C-Miras), the $`K`$-magnitudes of nlmc3-306907 and nlmc3-288533 were reported in Table 4 of Trams et al. (1999). The former is an oxygen-rich star and the latter is a carbon-rich star. Their results were shown in Table 2 together with $`K_{shv}`$.

Peak-to-peak modulation was found in the light curve of nlmc3-288533. The magnitudes tabulated in Table 2 were brighter than $`K_m`$, 11.35 and 11.20, of each star. The differences of magnitudes at the light maxima and $`K_m`$ were 1.00 and 0.55 mag respectively. Such differences can be accepted as the light variations of this type of variable stars.

When we adopt the $`K`$-magnitude by Trams et al.(1999), $`K_T`$, in Fig. 12, the position of nlmc3-306907 moves upward and seems to belong to the carbon-rich stars. However, this star should not be regarded as a typical oxygen-rich star because of the unusual colour variation it showed as mentioned in Section 9.4. The position of nlmc3-288533 in Fig. 12 will be slightly lower than the current position.

Figure 14 shows a similar diagram using the $`K`$-magnitude ($`K_S`$ : $`2.15\mu `$m) of the DENIS catalogue (http://cdsweb.u-strasbg.fr/denis.html). Among 146 samples, 139 stars have found, and $`K`$-magnitudes were given for 119 stars in the DENIS catalogue. The correlation between $`K_S`$ and $`K_m`$ is given by the equation,

$$K_S=0.9(\pm 0.09)K_m+0.73(\pm 1.07)$$
(16)
with a standard deviation of 0.47. Such a scatter may be explained as the results of intrinsic light variability. In Fig. 14, oxygen-rich and carbon-rich stars were located in the lower and the upper regions as seen in Fig. 12. The star nlmc3-306907 was again located in the carbon-rich star region. No $`K_S`$ data for nlmc3-288533 has been found. If nlmc3-306907 is a typical oxygen-rich star, this indicates that some oxygen-rich stars may exist at the upper side of the oxygen-rich relation in Fig. 12. Therefore brighter oxygen-rich stars than given by equation (14) may exist.

We show a ($`\mathrm{log}P,K_m`$) diagram for our sample of 146 stars in Fig. 15. In this figure, the solid diagonal line is the LMC PL relation for O-Miras of Hughes & Wood (1990) and the dotted line shows the track for O-Miras of Feast et al. (1989). Three groups, designated $`A`$, $`B`$ and $`C`$ are evident in Fig. 15. They are discussed in detail in Section 11. The stars in the most populous group lie close to these lines. nlmc3-306907 is located in the same region as carbon-rich stars. Another star, nlmc2-64521, is also located below the established relation. Therefore it is clearly hard to decide the all star below these lines are carbon-rich stars.

One of the 13 carbon-rich stars, nlmc2-64540, lies on these lines as an exception, and the remaining 12 carbon-rich stars fall below these lines. The slope of $`K_{shv}`$-magnitude versus MOA Colour (see Fig. 12 and formula (15)) for carbon-rich stars is steeper than that for oxygen-rich stars, so these lines intersect. This was why this exceptional carbon-rich star was located near the line of formula (14) and lies on the sequences for O-Miras in Fig. 15 (also in Fig. 13).

We note it was possible that at least 50 $``$ 60 carbon-rich stars were included in the sample of 146 stars because the ratio of oxygen-rich (the sample of 21 stars) and carbon-rich (the sample of 13 stars) stars was approximately $`62\%`$ and $`38\%`$ respectively. These could generate the extra scatter in the ($`\mathrm{log}P,K_m`$) plot.

When we converted the MOA data into $`K`$ magnitudes by using equation (15), most of stars including 11 carbon-rich stars (around group $`C`$ in Fig. 15) were found near the empirical LMC PL relation for C-Miras, but some of them remained in the fainter domain of the diagram. The isolation of the short period variable stars (seven + mark stars and a $``$ star, defined in the next section) became clearer in this case.

## 11 DISCUSSION

### 11.1 $`\mathrm{log}PK_m`$ diagram

It is difficult to estimate $`K`$-magnitudes from the large photometry database as it was originally designed to detect microlensing events. However we have now succeeded in estimating $`K`$-magnitudes using the properties of oxygen-rich and carbon-rich stars as described in the previous Section. Based on these results, we can investigate the ($`\mathrm{log}P,K_m`$) diagram more carefully.

As described in Section 10, our conversion of $`V_m`$ and $`R_m`$ into $`K_m`$-magnitude could not be applied for carbon-rich stars and they are located below the main sequence of Figures 13 and 15. Therefore the scatter in the ($`\mathrm{log}P,K_m`$) plane (Fig. 15) must be associated with these stars. However, even taking the above fact into account, the scatter in the plot is too large to be described by one sequence. To demonstrate this point, in Fig. 16 we show a histogram of vertical distance from the line of O-Mira sequence by Hughes & Wood (1990):

$$K=3.86(\mathrm{log}P2.4)+11.5$$
(17)
The error bars are equal to the square root of the number of stars in each of the three groups indentified in Figure 16.

A vertical separation of zero means it was on the Hughes & Wood (1990) track in Fig. 15. The most populous group ,$`A`$, lies at this position. This group contained the stars with largest amplitude. There was a small peak (group $`B`$, eight stars) around $`+1.5`$ mag (corresponding to a leftward shift in period in Fig. 15). The group $`B`$ stars are relatively small amplitude variables. Another peak which is relatively broad can be identified at around $`0.6`$ mag (group $`C`$, corresponding to a rightward shift in period in Fig. 15). We define each group in Fig. 16 using the following limits; the stars distributed between 0 and 0.6 mag as group $`A`$, the stars distributed greater than 0.9 mag are as group $`B`$, and stars distributed between $`0.9`$ and $`0.3`$ mag as group $`C`$.

In Fig. 17 we show the amplitude distribution, $`\delta R_m`$ for the stars in each group. In Fig. 18 we show the ($`\mathrm{log}P,\delta R_m`$) diagram with the same symbols as Fig. 15. In Fig. 17, it is clear that most of largest amplitude ($`\delta R_m2.2`$ mag) variable stars are group $`A`$ stars, which provides strong proof that these stars are Miras.

In contrast, the amplitudes of stars in group $`B`$ were smaller ($`\delta R_m1.4`$ mag) than most of the others as shown in Fig. 17 and Fig. 18. This group includes the seven short period stars and nlmc3-214266. This group is located on the very short period side ($`\mathrm{log}P2.1`$) of the ($`K_m,\mathrm{log}P`$) plane (Fig. 15). The group of seven short period stars is located in similar positions in the period-colour or period-magnitude diagrams (Fig. 7, 8, 9). All of these stars showed relatively inferior periodicity: the mean value of $`\theta `$ was $`0.43`$, compared with the mean of $`\theta 0.38`$ for the stars in group $`A`$. We obtained periodic variable stars via our selection criteria (Section 5), but we did not strictly eliminate stars whose periodicities were relatively poor. So the stars in group $`B`$ were somewhat semi-periodic compared to the stars in group $`A`$ (see Fig. 5).

### 11.2 The Group $`B`$ Stars

To examine $`K`$ magnitudes of the short period variable stars, we compared our estimates with $`JHK`$ observations made by the 1.9 m telescope at the South African Astronomical Observatory. Star identification, magnitudes and date of observations are listed in the first five columns of Table 3. Mean $`K`$ magnitudes derived by equation (14) are tabulated in the 6th column. The first three stars of the table are stars (a), (b), and (g) of group $`B`$, and the last one was the probably non-AGB star ($``$ in Fig. 15) of group $`B`$ (see Table 1). The differences in $`K`$ mag obtained at the SAAO and derived averaged $`K`$-magnitude, $`K_m`$, were 0.33, greater than 1.02, 0.06, and 0.04, respectively. Except for the second star, nlmc1-239264, the differences do not seem remarkable. We predicted successive light maxima of these four stars using the observed light curves, and have found that the SAAO observations were done near a predicted light maximum for the star nlmc1-324832, and near their predicted light minima for the other three stars. If their light variations were sufficiently regular, the SAAO $`K`$-magnitudes of nlmc1-324832 had to be brighter than 10.47. This contradicts the observed brightness, so that the regularity of the light variation should be suspected. For the other three stars, $`K`$-magnitudes were fainter than their $`K_m`$. The observed magnitude for nlmc1-239264 was fainter than our adopted mean magnitude, and the difference was 1.02 mag. Such a difference looked too large when compared with the other two stars. However, the true difference between $`K_{SAAO}`$ and the mean $`K`$ would be less than 1.02 mag because the adopted mean magnitude should be brighter than true mean magnitude. Therefore the difference may be acceptable with the true amplitude in $`R_m`$ greater than 1.38 mag. The question remains as to whether these differences originate from the uncertainty in $`K_m`$ or the semi-regular nature of these variables. As indicated in Fig. 14, the $`K_S`$ magnitudes for seven stars in group $`B`$ are were located along the line of the O-Miras. This tends to support the use of equation (14) in estimating the $`K`$-magnitudes. Further analysis of our database will resolve this question.

Amongst the eight stars of group $`B`$, nlmc3-214266 should be a massive star, not a star evolved from the red giant branch or the first red giant branch as pointed out in Section 9.1. We examine the properties of the other seven stars. The differences between the brightness of the group $`B`$ stars and of the group$`A`$ stars is evident in Fig. 15. The faintest star, nlmc1-118904, was equally bright in $`K_S`$ with the other six stars of group $`B`$. The amplitude of nlmc1-239264 was larger than the other six stars, but the position on the period-luminosity diagram was not peculiar. We have to study the nature of this star by comparison with more abundant samples to determine whether or not this star is different from the others.

In Fig. 15, the stars in group $`B`$ were located parallel to the sequence of group $`A`$ (corresponding to the Mira sequence of Hughes & Wood (1990)) and separated by $`\delta \mathrm{log}P0.4`$. The seven stars in group $`B`$ were located at a similar position to the short-period red (SP-red) stars which were reported by Feast & Whitelock (2000) analyzing Galactic Miras. While the SP-red stars had large amplitude and Mira-like properties, our group $`B`$ stars had small amplitudes except nlmc1-239264 and showed less regular variations, different from the Galactic Miras. The definition of the Galactic Miras is that their amplitude in standard $`V`$-light should be larger than 2.5 mag. It is a natural classification because there are few stars with such amplitudes. However, such a gap was not distinct in $`I`$-band observations of long period variables in the LMC (Whitelock, 1997). If we try to define Miras in the LMC by Fig. 17, the boundary might be at around 2.2 mag, though the small sample size makes that inconclusive. Therefore we discuss the long period variables of the LMC with no sharp distinction between Miras and semi-regulars.

The theoretical models of long-period variables predict the ratio of fundamental to first or second overtone period should be $`\delta \mathrm{log}P0.30.4`$. The stars in group $`B`$ were very similar to those compared with the higher mode pulsators by Wood and Sebo (1996). Higher mode pulsators are also expected to have smaller amplitudes than fundamental mode pulsators. The smallness of the amplitudes of these stars in Fig. 17 and Fig. 18 supported this point. Therefore if the stars in group $`A`$ were fundamental pulsators, the seven stars in group $`B`$ could be first or second overtone pulsators. On the other hand, as described above, these stars might be the counterparts of short-period red (SP-red) stars reported by Whitelock, Marang, and Feast (2000). These authors point out that the kinematics of such stars were different from the blue stars with similar periods. Such differences indicate a considerably different stellar age and imply a difference of evolutionary process due to a distinct chemical composition. A further development in the theory of stellar structure is needed before it will be possible to judge whether these seven stars are counterparts of SP-red stars (and are higher mode pulsators) or not.

To study these stars further we plotted, in Fig. 19, the $`V_mR_m`$ colour against $`K_m`$ magnitude with the symbols being the same as in Fig. 15. The seven stars of group $`B`$ are located at the reddest and highest luminosity region where longer-period variables are also located. This suggests that the radius of these stars is very similar to that of longer-period stars. The difference in the period of stars of similar radius is explained by the difference in the mass $`M`$ and the pulsation property $`Q`$. The estimate of the pulsation property for these stars has not been completed because of the uncertainty in the pulsation theory, but the effect of the chemical composition and the stellar mass on $`Q`$ looks negligible for a given mode (Barthes and Tuchman 1994). The difference in the period-temperature relation found in the globular clusters (Whitelock 1986) has been resolved using different chemical compositions, but the difference in the period of the stars of which the colour is the same cannot be explained this way. The difference in the stellar mass and the mode of pulsation should be studied to explain the difference in the period. When the distribution of periods with stellar physical properties is continuous, it is reasonable to suppose that this is due to changes in stellar mass or other stellar parameters such as chemical composition. Since the periods of these seven stars are separated from their long-period counterparts, the enhancement of a different mode is suggested. Further progress in data reduction is required to examine the existence of the gap in the periods of this luminous red domain of the colour-$`K`$ diagram.

When we consider the excitation of different modes, we have to study the reason for such a mode selection. In the classical Cepheids and RR Lyrae stars, the higher mode pulsators have been found on the blue side, or small stellar radii side, of the colour-magnitude diagram. Therefore it seems difficult to compare directly the group $`B`$ stars with higher mode pulsators of the Cepheid Instability Strip. On the other hand, it may be interesting to compare the surface gravity between them. The surface gravity of the higher mode pulsators is high for stars in the Cepheid Instability Strip. This implies that the surface gravity of the group $`B`$ stars might be also high. The critical frequency of the stellar atmosphere (i.e., the boundary frequency of reflection or transmission at the stellar surface) must be high for a high surface gravity. The enhancement of high frequency, or short period pulsation in group $`B`$ stars suggests high stellar mass when we suppose the same situation for the transmission and reflection for all pulsating stars.

### 11.3 The Group $`C`$ Stars

There are 46 stars in group $`C`$. Many carbon-rich stars might be included in this group. 12 group $`C`$ stars are the same as the carbon-rich stars tabulated in Hughes & Wood (1990). There was a possibility that about 50 $``$ 60 stars among the sample of 146 are carbon-rich stars as mentioned in Section 10. These carbon-rich stars were estimated to be too faint in our $`K_m`$ conversion. Also, in the amplitude distribution (Fig. 17), there was no obvious difference between group $`A`$ and group $`C`$. We conclude most group $`C`$ stars were carbon-rich and should be included in group $`A`$. However, we need to distinguish oxygen-rich stars from carbon-rich stars in order to investigate whether a long-period sequence exists, and if the strip structure is real or not. Future study of the relationship using a $`B`$-band filter should help resolve this issue.

### 11.4 Strip Structure

Wood (2000) presented the sequential structure of the PL diagram using a large sample. Because our sample size is small, it was hard to confirm this structure, except for the existence of a separate group of red short-period stars. However, it should be noted that these red short-period stars were of small amplitude ($`\delta R_m`$) and less regular (small $`\theta `$). These will be useful characteristics in the study of the nature of these stars. In the period-$`K`$ diagram, the stars, nlmc3-214266 and nlmc1-239264, were also contained in group $`B`$ implying that it is necessary to investigate the properties of individual star even in large scale photometry in order to understand the period-luminosity relation. The individual star discovered by Wood (2000) should also be checked to determine whether it is to be included in group $`B`$ and the in the set of Galactic red SP stars.

## 12 SUMMARY

We selected long-period red LMC variables from the three year MOA database, and executed a careful period analysis using the PDMM and PDM codes. We selected 146 variable stars using criteria which focus on long-period red variables (Miras and SRs), and investigated their properties. Their observational data were tabulated.

In the ($`\mathrm{log}P,V_m`$) and ($`\mathrm{log}P,R_m`$) diagrams, the luminosity of stars in the major cluster tended to decrease as the period increases. In the ($`\mathrm{log}P,V_mR_m`$) diagram, the colours were redder as the periods were longer.

We transformed $`V_m,R_m`$ magnitudes to $`K_m`$ values using 19 oxygen-rich stars, and presented the period-luminosity (PL) relation in the ($`\mathrm{log}P,K_m`$) plane. Many of the stars in our sample, including the large amplitude stars, fell near the position corresponding to the O-Mira sequence by Feast et al. (1989) and Hughes & Wood (1990). A small number of stars (8) formed a separate group in the period-magnitude or period-colour diagrams. This group (group $`B`$) was found in the bright and short period region of the ($`\mathrm{log}P,K_m`$) plane compared with the O-Mira sequence. These stars showed less periodic variation, and six of them had small amplitudes.

The separation in $`\mathrm{log}P`$ between them and the main group was $``$ 0.4. They were on the red side of the period-colour diagram and were located at the reddest and brightest tip of the colour-$`K_m`$ diagram similar to the long-period members. If our conversion scale is adequate, it indicates that some of the stars pulsate in a short period compared with the stars which may have a similar radius.

These facts were possibly explained by the assumption proposed by Wood & Sebo (1996): that those six stars with short periods and small amplitudes were higher mode pulsators than the stars on the O-Mira sequence. Two particular stars, nlmc3-214266 and nlmc1-239264, were included in group $`B`$. These stars showed different features from the other six variable stars. Their pulsation property possibly differs from the other six stars.

In the ($`\mathrm{log}P,K_m`$) diagram, below the O-Mira sequence, many stars were located in a wide region. Since almost all identified C-Miras were found in this domain, we supposed that a large number of such variable stars were C-Miras which were more luminous than the magnitude indicated here. In order to reduce the systematic scattering in the ($`\mathrm{log}P,K_m`$) diagram, it was necessary to distinguish carbon-rich stars from oxygen-rich stars in advance.

The difficulty in distinguishing oxygen-rich and carbon-rich stars comes from the fact that our blue filter covers both the standard $`B`$ and $`V`$ bands. It is expected that in the future these stars will be observed in the $`B`$ or $`V`$ bands which will distinguish between carbon-rich stars and oxygen-rich stars. Alternatively, using any passband which avoids the effect of TiO bands to estimate the bolometric magnitude of the stars could be used. Improvement in the magnitude calibration and confidence in the interpretation of the plural sequences is expected to follow from the analysis of further large observational databases.

## 13 Acknowledgments

We are grateful to A. Tanaka for his collaboration in preparing the PDM code. We wish to thank P.A.Whitelock for her useful comments and help. We also thank H. Saio for helpful discussions. This work is supported by a grant-in-aid for scientific research (A) of the Japan Ministry of Education, Science, Sports and Culture, and also by the Marsden Fund of the Royal Society of New Zealand.


format: text

text:
# Surface critical exponents at a uniaxial Lifshitz point

## I Introduction

Critical phenomena at surfaces have been extensively studied theoretically during the last three decades.  The surface phase diagram of the three-dimensional semi-infinite Ising model with only nearest neighbor couplings is well established. Introducing two different ferromagnetic interactions depending on whether the neighboring spins are both located at the surface, $`J_s0`$, or not, $`J_b>0`$, two typical scenarios are encountered. If the ratio of the surface coupling $`J_s`$ to the bulk coupling $`J_b`$, $`r=J_s/J_b`$, is smaller than a critical value, $`r_{sp}1.50`$ for the simple cubic lattice,  the system undergoes at the bulk critical temperature $`T_c`$ an ordinary transition, with the bulk and surface ordering occurring at the same temperature. Beyond this critical value one first observes the surface transition where the surface alone orders at a temperature $`T_s>T_c`$, followed by the extraordinary transition of the bulk at $`T_c`$. At the critical ratio $`r_{sp}`$, the special transition point is encountered, displaying critical properties distinct from those of the ordinary or the surface transition. The critical exponents of the different surface universality classes are known with rather high precision. However, far less is known on the surface critical behaviour of models with competing interactions which is the subject of the present paper.

The axial next-nearest-neighbor Ising (ANNNI) model  is the best known of these models, its bulk properties being intensively investigated since many years. Here, competition between ferromagnetic nearest neighbor and antiferromagnetic next-nearest neighbor couplings takes place in one direction. The Hamiltonian of the three-dimensional version, defined on a cubic lattice, may then be written in the form

$``$ $`=`$ $`J_b{\displaystyle \underset{xyz}{}}s_{xyz}\left(s_{(x+1)yz}+s_{x(y+1)z}+s_{xy(z+1)}\right)`$ (2)
$`+\kappa _bJ_b{\displaystyle \underset{xyz}{}}s_{xyz}s_{xy(z+2)}`$
where $`J_b>0`$ and $`\kappa _b>0`$ are coupling constants. The direction of competing interactions (here, the $`z`$-direction) is also called axial direction.

Spatially modulated phases due to competing interactions are observed, among others, in magnetic systems, alloys or ferroelectrics.  Some compounds (as, for example, BCCD  or NaV2O5 ) display very rich phase diagrams with a multitude of commensurately or incommensurately modulated phases. A further remarkable feature of these systems is the possible existence of a special multicritical point called Lifshitz point.  At a Lifshitz point, a disordered, a uniformly ordered and a periodically ordered phase become indistinguishable. Various systems have been shown to possess a Lifshitz point.

Recently, there has been renewed interest in the uniaxial Lifshitz point encountered in the phase diagram of the three-dimensional ANNNI model. This strongly anisotropic equilibrium critical point, located at $`\kappa _b^L=0.270\pm 0.004`$ and $`k_BT_c^L/J_b=3.7475\pm 0.0005`$,  is characterized by the anisotropy exponent $`\theta =\nu _{}^L/\nu _{}^L`$ where $`\nu _{}^L`$ and $`\nu _{}^L`$ are the exponents of the bulk correlation lengths parallel and perpendicular to the $`z`$-axis. The value $`\theta 0.49`$ has been obtained in a second-order $`ϵ`$-expansion. Recent thorough determinations of Lifshitz point critical exponents using different techniques (field-theoretical calculations  and Monte Carlo methods ) have yielded excellent agreement. In addition, the scaling of the spin-spin correlator at the Lifshitz point was determined  and its form shown to agree with the prediction resulting from a generalization of conformal invariance to the strongly anisotropic scaling at the Lifshitz point.

Whereas the properties of the bulk ANNNI model are well established, far less is known about the behaviour of the ANNNI model in the presence of surfaces. Early Monte Carlo studies  investigated ANNNI samples with free surfaces but did not pay special attention to surface properties as they were exclusively interested in bulk behaviour. Thin ANNNI films with free surfaces were studied recently and they were shown to present a distinct phase diagram for every film thickness. An early attempt to study surface critical behaviour near the Lifshitz point was undertaken by Gumbs  in the framework of mean-field theory. A much more elaborated mean-field treatment of the semi-infinite ANNNI model was presented recently by Binder, Frisch, and Kimball.  They considered two different surface orientations with respect to the axis of competing interactions: surfaces oriented perpendicular  or parallel  to this axis. For both cases the mean-field surface critical exponents at the Lifshitz points were determined and two different sets of critical exponents were obtained. This dependency of the critical behaviour on the surface orientation may be explained by the anisotropic scaling at the Lifshitz point. These authors also predicted the existence of a surface transition at the Lifshitz point where only the surface orders, in variance with Gumbs  who claimed that at the Lifshitz point the surface could not order before the bulk. It is worth noting that no results beyond mean-field theory are available on the surface critical behaviour of the ANNNI model close to the Lifshitz point.

In this paper I present the first Monte Carlo study of the surface critical behaviour at the uniaxial Lifshitz point encountered in the ANNNI model. This study uses a cluster flip algorithm, especially designed for simulating models with competing interactions, which has been successfully employed in a recent investigation of the bulk ANNNI model.

The paper is organized in the following way. The next Section is devoted to the presentation of semi-infinite ANNNI models with two different surface orientations. Some details of the numerical method are also presented. Section III deals with surfaces perpendicular to the axis of competing interactions, paying special attention to the surface phase diagram and to the determination of surface critical exponents. Surfaces oriented parallel to this special direction are treated in Section IV. A short summary and outlook conclude the paper.

## II Models and method

Due to the anisotropy of the ANNNI model, surfaces with different orientations are not equivalent. The two surface orientations with respect to the axial direction considered in this paper are the followings (see Figure 1): surfaces perpendicular to the axis of competing interactions (case A) and surfaces parallel to this axis (case B). The same orientations were treated in recent mean-field studies.

For case A modified surface couplings with strength $`J_s`$ connecting neighboring surface spins are introduced in addition to the usual ANNNI interactions, see Eq. (1) and Figure 1a. Three different scenarios have to be distinguished, depending on the value of the bulk competition parameter $`\kappa _b`$. When $`\kappa _b`$ is smaller than the Lifshitz point value $`\kappa _b^L`$, the bulk undergoes a second order phase transition between the disordered high temperature phase and the ordered, ferromagnetic, low temperature phase at the critical temperature $`T_c(\kappa _b)`$. This phase transition belongs to the universality class of the 3D Ising model. Consequently, the surface phase diagram will resemble that of the 3D semi-infinite Ising model, with a possible shift of the location of the special transition point, $`r_{sp}(\kappa _b)`$, as function of $`\kappa _b`$. At the Lifshitz point, $`\kappa _b=\kappa _b^L`$, the recent mean-field treatment yields a surface phase diagram for the semi-infinite ANNNI model similar to the Ising model, but with critical exponents which differ from those of the Ising model. The conjectured existence  of a surface transition at the Lifshitz point, and therefore of a special transition point, is in agreement with the Monte Carlo results presented in the next Section. Finally, for strong axial next-nearest neighbor bulk couplings, $`\kappa _b>\kappa _b^L`$, the bulk phase transition belongs to the universality class of the 3D $`XY`$ model. Therefore, for weak surface couplings, the ordinary transition critical behaviour should be identical to that of the three-dimensional semi-infinite $`XY`$ model, whereas for strong surface couplings the surface effectively decouples from the bulk and a two-dimensional bulk Ising critical behaviour is expected at the surface transition. These two critical lines meet at a multicritical point which should have interesting properties. However, I will in this study not consider the latter case as I am mainly interested in the surface critical behaviour at the Lifshitz point.

Case B with surfaces oriented parallel to the axial direction (see Figure 1b) is the most interesting but also the most demanding case. Introducing modified nearest neighbor, $`J_s`$, and axial next-nearest neighbor couplings, $`\kappa _s0`$, in the surface layer leads to intriguing and very complex situations. For example, a multicritical point shows up where an ordinary transition with a modulated bulk meets a surface transition to a floating incommensurate phase in the surface layer. In this study I will not discuss these pecularities but concentrate exclusively on the case $`\kappa _s\kappa _b^L`$. For these values of the surface competition parameter $`\kappa _s`$ the corresponding two-dimensional ANNNI model presents a transition, belonging to the 2D Ising universality class, from the disordered phase to the low temperature ferromagnetic phase.

Surface quantities are determined by simulating ANNNI films consisting of $`L_x\times L_y\times L_z`$ spins, with free boundary conditions perpendicular to the surfaces and periodic boundary conditions otherwise. Hence, for case A the system consists of $`L_z`$ layers with $`L_x\times L_y`$ spins per layer, whereas for case B the samples are formed by $`L_y`$ layers containing $`L_x\times L_z`$ spins. The semi-infinite models are then obtained in the limit $`L_x`$, $`L_y`$, $`L_z\mathrm{}`$. Some care is needed when choosing the shape of the finite ANNNI films in the vicinity of the Lifshitz point. On the one hand, the special finite-size effects coming from the anisotropic scaling at this point  are best taken into account  by choosing anisotropic samples with an increased number of sites perpendicular to the axial direction. On the other hand, computation of surface quantities is usually done by simulating samples where the linear dimension of the surfaces exceeds the film thickness. This does not pose any problem for case A, as here the axial direction coincides with the direction perpendicular to the surfaces. Therefore, for this case, systems of anisotropic shape with $`L_z`$ layers and $`L_x^2`$ spins per layer are simulated, $`L_z`$ ranging from 10 to 120 and $`L_x`$ from 10 to 60. For case B, however, surfaces are parallel to the axial direction. In order to balance the competing finite-size effects, I have chosen to simulate cubes with $`L_x^3`$ spins, $`L_x`$ ranging from 10 to 80. For comparison, a few simulations have also been done for samples with anisotropic shapes.

Critical phenomena are best studied numerically with non-local Monte Carlo methods. Recently, we proposed a cluster-flip algorithm, based on the one-cluster flip algorithm of Wolff,  especially designed for simulating spin systems with competing interactions.  Below, this algorithm, which combines the cluster algorithm approaches for systems with long-range ferromagnetic couplings  with that for systems with nearest neighbor random couplings, is reformulated in order to take modified surface couplings into account.

Let $`i`$ be a lattice site characterized by the spin $`s_i`$ and already belonging to the cluster we iteratively build up. A nearest neighbor site $`j`$ with spin $`s_j`$ is added to the cluster with probability

$$\frac{1}{2}\left(1+\text{sign}\left(s_is_j\right)\right)\left(1\mathrm{exp}\left[2J_s/(k_BT)\right]\right)$$
(3)
if sites $`i`$ and $`j`$ are both located in the surface layer, and with probability

$$\frac{1}{2}\left(1+\text{sign}\left(s_is_j\right)\right)\left(1\mathrm{exp}\left[2J_b/(k_BT)\right]\right)$$
(4)
otherwise. A lattice site $`k`$ with spin $`s_k`$ axial next-nearest neighbor to $`i`$ is included with probability

$$\frac{1}{2}\left(1\text{sign}\left(s_is_k\right)\right)\left(1\mathrm{exp}\left[2J_s\kappa _s/(k_BT)\right]\right)$$
(5)
for case B if both spins are surface spins, and with probability

$$\frac{1}{2}\left(1\text{sign}\left(s_is_k\right)\right)\left(1\mathrm{exp}\left[2J_b\kappa _b/(k_BT)\right]\right)$$
(6)
otherwise. The final cluster, which is flipped as a whole, presents two pecularities which are worth mentioning. First, it contains spins of both signs, in variance with the Wolff cluster method, where all spins belonging to a cluster to be flipped have the same sign. Second, the cluster spins are not always connected by nearest neighbor bonds, which is again different from the traditional cluster flip method. This algorithm works well in the vicinity of the Lifshitz point, as demonstrated in the recent numerical study of the bulk Lifshitz point critical behaviour.

During the simulation, layer-dependent quantities are computed. I discuss these quantities here only for case A, the corresponding quantities for case B are obtained by replacing $`z`$ by $`y`$ and $`L_z`$ by $`L_y`$. Of great interest when studying surface properties are the magnetization per layer

$$m(z)=\frac{1}{L_xL_y}\left|\underset{xy}{}s_{xyz}\right|$$
(7)
and the susceptibility per layer

$$\chi (z)=\frac{L_xL_y}{k_BT}\left[\left(\frac{1}{L_xL_y}\underset{xy}{}s_{xyz}\right)^2\left(m(z)\right)^2\right].$$
(8)
The surface magnetization is then $`m_1=m(z=1)=m(z=L_z)`$, whereas the response of the surface magnetization to a surface field is $`\chi _{11}=\chi (z=1)=\chi (z=L_z)`$. From the profiles $`m(z)`$ and $`\chi (z)`$ one also obtains the surface excess quantities

$$m_s=\underset{z=1}{}\left(m_bm(z)\right)$$
(9)
and

$$\chi _s=\underset{z=1}{}\left(\chi _b\chi (z)\right)$$
(10)
where $`m_b`$ and $`\chi _b`$ are the bulk magnetization and the bulk susceptibility. The response of the surface to a bulk field, $`\chi _1`$, the energy, the specific heat, and the Binder cumulant have also been computed.

Thermal averages are extracted by generating, after equilibration, $`5\times 10^5`$ clusters, and error bars are obtained by averaging over at least ten different realizations using different random numbers.

## III Surfaces perpendicular to the axial direction

Semi-infinite ANNNI models with surfaces oriented perpendicular to the direction of competing interactions exhibit a crossover from Ising surface critical behaviour for $`\kappa _b=0`$ to Lifshitz point surface criticality for $`\kappa _b=\kappa _b^L`$. For the purpose of studying this crossover, simulations were not only done for the value of $`\kappa _b`$ at the Lifshitz point, $`\kappa _b^L=0.27`$, but also at $`\kappa _b=0.15`$ and $`\kappa _b=0.24`$. The corresponding critical temperatures are listed in Table I. Note that the same values of $`\kappa _b`$ have been used in a study of ANNNI bulk critical behaviour.

Various values of the surface coupling strength, $`J_s`$, have been considered, with $`J_s`$ ranging from 0 to $`3J_b`$. The resulting surface phase diagrams are displayed in Figure 2, together with the corresponding phase diagram of the semi-infinte Ising model, $`\kappa _b=0`$. The horizontal lines indicate the bulk critical temperatures which decrease with increasing $`\kappa _b`$. One also observes a shift of the location of the special transition point to lower values of $`r=J_s/J_b`$ for larger values of the competition parameter $`\kappa _b`$, in accordance with the recent mean-field results. This shift is mainly due to the decrease of the bulk critical temperature. The different surface transition lines in Figure 2 merge for $`J_sJ_b`$ and form one line independent of the value of $`\kappa _b`$. This follows from the fact that for very strong surface couplings the surface effectively decouples from the bulk. At the Lifshitz point value $`\kappa _b=0.27`$, the phase diagram is also similar to the Ising case, with an ordinary transition, a surface transition, an extraordinary transition, and a special transition point. This is in agreement with the mean-field treatment of Ref. , but disagrees with the results obtained in Ref. . Based on our data, the special transition point of the semi-infinite ANNNI model at the Lifshitz point is located at $`r_{sp}^L=1.15\pm 0.05`$.

Crossover phenomena showing up in the surface magnetization may be analysed by plotting the effective exponent

$$\beta _{1,eff}=\text{d}\mathrm{ln}m_1/\text{d}\mathrm{ln}t$$
(11)
where $`t=\left(T_cT\right)/T_c`$ is the reduced temperature. On approaching $`T_c`$, $`t0`$, $`\beta _{1,eff}`$ becomes the critical exponent $`\beta _1`$ of the surface magnetization. In the following, data unaffected by finite-size effects are usually displayed, finite-size dependences being circumvented by adjusting the size of the sample.

Effective surface magnetization exponents obtained for $`\kappa _b=0`$, 0.15, and 0.24, with $`J_s=J_b`$ are shown in Figure 3. Clearly, for a fixed value of $`t`$, $`\beta _{1,eff}`$ decreases with increasing $`\kappa _b`$. In the limit $`t0`$, however, the effective exponents all tend asymptotically to the critical surface magnetization exponent $`\beta _10.80`$ of the semi-infinite Ising model at the ordinary transition, in accordance with the Ising character of the ANNNI bulk transition for $`\kappa _b<\kappa _b^L`$. The observed increase of the corrections to scaling with $`\kappa _b`$ reflects not only the closer proximity of the bulk Lifshitz point but also the reduced distance of $`r=J_s/J_b=1`$ to the special transition point $`r_{sp}(\kappa _b)`$. This is illustrated in Figure 3 by the data obtained for $`\kappa _b=0.24`$ with $`J_s=0.5J_b`$ where corrections to scaling are greatly reduced compared to the case $`J_s=J_b`$

Some of the results obtained for $`\kappa _b=\kappa _b^L`$, with $`r<r_{sp}^L`$, are displayed in Figures 4 and 5. For this choice of the interactions, the surface undergoes an ordinary transition at the Lifshitz point critical temperature $`T_c^L`$. Both the magnetization per layer, shown in Figure 4, and the susceptibility per layer (not shown) vary non-monotonously close to $`T_c^L`$. Whereas a non-monotonic behaviour of the layer susceptibility is also observed close to the ordinary transition in thick Ising films, the layer magnetization increases monotonically in Ising films from its surface value to the bulk value. The non-monotonic behaviour of the layer magnetization close to the Lifshitz point ordinary transition may be explained in the following way. Inside the bulk the value of the magnetization mainly results from the balancing of competing influences in axial direction. Close to the surface, however, some negative contributions are missing, which then yields a tendency to stronger local ordering, resulting in a maximum of the layer magnetization near the surface. Similar non-monotonic profiles also show up in the semi-infinite critical Ising model in the presence of a weak surface field. It should be noted that in our case the appearance of non-monotonic profiles depends strongly on the value of the competition parameter $`\kappa _b`$. This behaviour is not observed for the studied values of $`\kappa _b`$ smaller than $`\kappa _b^L`$. For $`\kappa _b=\kappa _b^L`$ the maximum in the magnetization profile shows up for all values of the surface coupling $`J_s`$ leading to the ordinary transition. When the temperature is increased, the maxima both of the layer magnetization and of the layer susceptibilty are shifted towards the center of the system.

The surface magnetizations and the corresponding effective exponents obtained for two different values of the coupling ratio $`J_s/J_b`$ are compared in Figure 5. As expected, the amplitudes (see Figure 5a) and the corrections to scaling (see Figure 5b) differ. It is only close to $`T_c`$, for $`t0.07`$, that the effective exponents (and therefore the corrections to scaling) become identical. Furthermore, the effective exponents then vary almost linearly with temperature. A linear extrapolation yields the value $`\beta _1^L=0.618\pm 0.005`$ for the surface magnetization critical exponent at the Lifshitz point ordinary transition. This value is clearly smaller than the value $`\beta _1=0.80\pm 0.01`$ obtained at the ordinary transition of the semi-infinite Ising model. It is worth noting that mean-field approximation yields the common value $`\beta _1^{MF}=1`$ for both cases.

Estimates of Lifshitz point surface critical exponents at the ordinary transition are gathered in Table II. Hereby, the critical exponent $`\gamma _1^L=0.84\pm 0.05`$ of the susceptibilty $`\chi _1`$, i.e. of the response of the surface to a bulk field, can be obtained by analysing the corresponding effective exponent in the same manner as discussed previously. The response of the surface to a surface field, $`\chi _{11}`$, however, is finite at $`T_c^L`$, with a cusp-like singularity when approaching the critical point from high temperatures. Instead of trying to extract the corresponding critical exponent $`\gamma _{11}^L`$ directly from $`\chi _{11}`$, it is easier to study the derivative $`\chi _{11}^{}=\text{d}\chi _{11}/\text{d}T`$ which in our case diverges on approach to the critical point with the power law $`\chi _{11}^{}(t)^{\gamma _{11}^L1}`$. This then yields the value $`0.06\pm 0.02`$ for $`\gamma _{11}^L`$.

Further Lifshitz point surface critical exponents are obtained by analysing the power law behaviour of the excess quantities $`m_s`$ and $`\chi _s`$, see Eq. (8) and (9), close to the critical point: $`m_st^{\beta _s^L}`$ and $`\chi _st^{\gamma _s^L}`$. Based on our data, the excess quantities critical exponents are estimated in the present case to be $`\beta _s^L=0.14\pm 0.04`$ and $`\gamma _s^L=1.69\pm 0.07`$.

A closer inspection of Table II reveals that various scaling relations are fulfilled, thus demonstrating the reliability of our estimates. Indeed, inserting the Lifshitz point bulk critical exponents $`\gamma _b^L=1.36\pm 0.03`$ and $`\beta _b^L=0.238\pm 0.005`$  as well as $`\nu _{}^L=0.348`$, the following scaling relations are readily verified:

$$\beta _s^L=\beta _b^L\nu _{}^L$$
(12)
$$\gamma _s^L=\gamma _b^L+\nu _{}^L$$
(13)
$$2\gamma _1^L\gamma _{11}^L=\gamma _s^L$$
(14)
Lifshitz point surface quantities have also been studied near the special transition point with $`J_s=1.15J_b`$, see above. At temperatures below $`T_c^L`$ magnetization profiles are extremely flat: only very close to the surface does one observe a slow increase of the layer magnetization which reaches a maximum in the layer immediately below the surface layer. The surface magnetization itself is slightly smaller than the bulk magnetization. Estimated values of Lifshitz point surface critical exponents at the special transition point are listed in Table II. As both the bulk and the surface become critical at the special transition point, the responses of the surface to a bulk field as well as to a surface field diverge. The corresponding critical exponents are estimated to be $`\gamma _{1,sp}^L=1.28\pm 0.08`$ and $`\gamma _{11,sp}^L=0.76\pm 0.05`$ where effective exponents have again been analysed. For the surface magnetization critical exponent the value $`\beta _{1,sp}=0.22\pm 0.02`$ is obtained. Error bars given for the special transition point critical exponents include the uncertainty in the localisation of this multicritical point.

It is worth noting that the Lifshitz point surface critical exponent at the special transition point given in Table II agree within the error bars with the estimates  of the corresponding exponents for the semi-infinite Ising model. The change from Ising to Lifshitz point bulk critical behaviour when $`\kappa _b\kappa _b^L`$ seems to have only a minor impact on the computed surface critical behaviour at the special transition point. However, this may not hold for the crossover exponent, which describes the merging of the surface transition line with this point, as different values are obtained in mean-field approximation for the Ising and for the Lifshitz point case. In the present work, no attempt was made to determine the crossover exponent.

## IV Surfaces parallel to the axial direction

Two different approaches have been chosen for studying the critical behaviour of surfaces oriented parallel to the axial direction. In the first approach the value of the surface competition parameter $`\kappa _s`$ is set to zero so that only the strength of the surface nearest neighbor coupling $`J_s`$ is varied. In the second approach the bulk and surface competition parameters have the same value, $`\kappa _s=\kappa _b`$, with $`J_s`$ ranging from 0 to $`6J_b`$. Simulations done for the latter case may be viewed as a preparatory work for the study of the more complex multicritical points mentioned in Section II. In this work I mainly discuss Lifshitz point surface criticality, i.e. $`\kappa _b=\kappa _b^L`$, at the ordinary transition and at the special transition point.

First one has to note that the critical value of the coupling ratio $`r_{sp}=J_s/J_b`$, needed for the surface to get critical by itself, now depends both on the bulk and on the surface competition parameters. For example, at the bulk Lifshitz point $`r_{sp}^L=1.75\pm 0.05`$ for $`\kappa _s=0`$, whereas for $`\kappa _s=\kappa _b^L`$ the special transition point is tentatively located at $`r_{sp}^L4.3`$. This remarkable shift of the critical coupling ratio to larger values in the latter case may be explained by a stronger weakening of the effectice surface couplings, due to the competing surface interactions, than of the effective bulk couplings.

At the bulk Lifshitz point, the magnetization per layer increases monotonically from its surface value to its bulk value, as shown in Figure 6. This is in marked contrast to the behaviour encountered when the surface is perpendicular to the axial direction, as discussed in the previous Section, see Figure 4. For the surface orientation considered here, the layers are parallel to the direction of competing interactions, so that no axial contributions are missing close to the surface and no tendency to stronger local ordering is expected.

Lifshitz point surface critical exponents at the ordinary transition, which have been computed in a similar way as for case A, are gathered in Table III. These values should be compared with the corresponding values of case A surfaces listed in Table II. Figure 7 shows effective exponents derived from the surface magnetizations obtained by simulating samples with different surface orientations. Clearly, both types of surfaces yield different asymptotic values $`\beta _1^L`$ for the effective exponent, thus demonstrating that the values of surface critical exponents depend on the surface orientation. Different values are also obtained for the critical exponents $`\gamma _{11}^L`$, $`\beta _s^L`$, and $`\gamma _s^L`$. For $`\gamma _1^L`$ the situation is not so clear as rather similar values follow from the analysis of effective exponents. One may note that mean-field theory yields for $`\gamma _1^L`$ the same value for both surface orientations.

Whereas the scaling relation (14) is fulfilled independently of the orientation of the surface, the scaling relations (12) and (13) have to be modified for surfaces parallel to the axial direction. Indeed, the behaviour of excess quantities is governed close to a bulk critical point by the bulk correlation length along the direction perpendicular to the surface. The Lifshitz point being an anisotropic critical point characterized by two correlation lengths diverging with different critical exponents, $`\nu _{}^L`$ respectively $`\nu _{}^L`$ has to be used for case B respectively A surfaces. This then yields for surfaces parallel to the axial direction the scaling relations

$$\beta _s^L=\beta _b^L\nu _{}^L$$
(15)
and

$$\gamma _s^L=\gamma _b^L+\nu _{}^L.$$
(16)
Using $`\nu _{}^L=0.709`$, a nice agreement with the numerical estimates for $`\beta _s^L`$ and $`\gamma _s^L`$ is observed.

I also studied the Lifshitz point surface critical behaviour of surfaces oriented parallel to the axial direction close to the special transition point, choosing $`\kappa _s=0`$ and $`J_s=1.75J_b`$, with $`\kappa _b=\kappa _b^L`$. The estimates for the critical exponents $`\beta _{1,sp}^L`$, $`\gamma _{1,sp}^L`$, and $`\gamma _{11,sp}^L`$ have been included in Table III. Again, the quoted error bars take into account the sample averaging as well as the uncertainty in the location of the special transition point. Interestingly, the values of the critical exponents are very similar to the values obtained for the other type of surfaces considered in Section III and quoted in Table II. This is a strong indication that at the bulk Lifshitz point the surface critical behaviour at the special transition point may not depend on the surface orientation.

## V Conclusion and outlook

Whereas surface criticality of models with only ferromagnetic couplings has been the subject of numerous investigations, studies of surface critical behaviour in systems with competing interactions are very scarce. The aim of the present work is to fill some gaps by presenting the first Monte Carlo study of surface critical behaviour in semi-infinite three-dimensional ANNNI models, paying special attention to the ordinary transition as well as to the special transition point, both at the bulk Lifshitz point. Prior works  exclusively focused on the mean-field treatment of this problem.

Overall, the numerical results, which have been obtained by studying different surface orientations with respect to the axis of competing interactions, are in qualitative agreement with results obtained in the recent mean-field treatments. The location of the special transition point is a function of the bulk competition parameter as well as of the surface competition parameter, when present. At the bulk Lifshitz point, the surface diagram is similar to that of the semi-infinite $`3D`$ Ising model, with a surface transition line ending in a special transition point. This is in accordance with recent mean-field results, but disagrees with the early treatment of Gumbs. Finally, at the bulk Lifshitz point, surface critical exponents at the ordinary transition depend on the orientation of the surface, as predicted by mean-field theory.

The main results of our study are the Monte Carlo estimates of surface critical exponents at the bulk Lifshitz point gathered in Tables II and III. This is the first time that values for these exponents have been determined numerically. Whereas at the ordinary transition surface critical exponents depend on the surface orientation, the computed values at the special transition point do not show a clear dependency on the orientation of the surface with respect to the axial direction. Furthermore, these values agree rather well with the values obtained in the Ising model at the special transition point. In order to clarify the situation at this multicritical point renormalization group calculations are desirable. However, with regard to the immense technical difficulties already encountered in the field-theoretical computation of bulk critical exponents at the Lifshitz point,  this seems to be a formidable task not easily achievable in the near future.

A host of open and interesting problems remains for future work. In the present study, surface criticality has only been investigated for values of the bulk competition parameter $`\kappa _b\kappa _b^L`$, thus neglecting the interesting case when the bulk undergoes a phase transition to a modulated structure. Similarly, for surfaces parallel to the axial direction, the surface competition parameter $`\kappa _s`$ has been chosen in such a way that at the surface transition a uniformly ordered phase appears in the surface layer. However, for large values of $`\kappa _s`$ one will observe a surface transition to a floating incommensurate phase in the surface layer. Various interesting scenarios follow, one of the most intersting being a multicritical point where a surface transition line with a floating phase in the surface layer meets an ordinary transition line with a modulated bulk.

A further interesting question concerns the scaling form at the Lifshitz point of spin-spin correlation functions along and perpendicular to the surface. A recent generalization of conformal invariance to the strong anisotropic scaling at the Lifshitz point has yielded predictions for the scaling form of bulk two-point correlators which were found to agree perfectly with Monte Carlo results. Similarly, predictions for two-point correlation functions in semi-infinite systems at the Lifshitz point can be obtained in the framework of this theory and may then be compared to numerical data. Work along the sketched lines is planned for the future.

###### Acknowledgements.

I thank Malte Henkel for useful discussions, and him and Alfred Hüller for a critical reading of the manuscript.