text:
# Obscured Asymptotic Giant Branch Variables in the Large Magellanic Cloud and the Period-Luminosity RelationThis paper is based on observations made at the South African Astronomical Observatory.

## 1 Introduction

Stars at the top of the Asymptotic Giant Branch (AGB) are of interest, primarily for two reasons. First, they are undergoing mass loss – via mechanisms which are still poorly understood. Secondly, they represent the most luminous phase in the evolution of low- and intermediate-mass stars and, as high resolution techniques enable us to examine stars in ever more distant galaxies, are therefore among the first individual objects to be isolated in a given population; furthermore, they will contribute a significant fraction of the infrared luminosity of distant unresolved populations (Ferraro et al. 1995; Mouhcine &Lançon 2002). The Large Magellanic Cloud (LMC) provides a particularly useful laboratory for the studies of luminous AGB stars at a known distance but with a range of initial masses.

Two factors complicate the studies of these luminous AGB stars - their long period large amplitude variability, which necessitates observations over periods of many years, and their high mass-loss rates which lead to thick dust shells and necessitate observations at mid-infrared wavelengths. Earlier papers in this series (Loup et al. (1997 - Paper I), Zijlstra et al. (1996 - Paper II), van Loon et al. (1997, 1998, Papers III and IV), Trams et al. (1999b - Paper V) and van Loon et al. (1999 - Paper VI)) described the identification of luminous AGB stars in the LMC, using IRAS and $`JHKL`$ photometry, as well as extensive ISO observations of a subset of them, and the derivation of their mass-loss rates.

This paper concentrates on the variability characteristics of the sample, reporting near-infrared ($`JHKL`$) photometry obtained in support of the spectroscopy and mid-infrared photometry discussed in previous Papers (I to VI). These data, together with others from the literature, are used to determine pulsation periods and amplitudes and to look for long-term trends in the behaviour of the stars. They are also combined with the observations discussed previously to determine the bolometric magnitudes and examine their relationship with the pulsation characteristics of the stars.

It is known that stars near to the end of their AGB evolution are large amplitude, reasonably regular pulsators and that stars with these characteristics form a well defined group (e.g. Whitelock 2002 and references therein). As the intention here is to consider stars very near to the tip of the AGB, we use the information at our disposal to isolate the objects of interest, removing from the sample supergiants as well as semi-regular and irregular pulsators.

Following Feast, Whitelock & Menzies (2002) we assume here that the distance modulus of the LMC is $`(mM)_0=18.60`$ mag, and, where necessary, the results of other analyses are corrected to this value.

### 1.1 Hot Bottom Burning (HBB)

Towards the end of AGB evolution, in stars with initial masses in excess of 4 or 5 $`M_{}`$ (or perhaps at even lower masses for very low metallicities), the base of the H-rich convective-envelope can dip into the H-burning shell; the introduction of fresh H-rich material into the nuclear-burning shell allows the luminosity to go above the core-mass luminosity predictions (Blöcker & Schönberner 1991). This process is known as envelope burning or hot bottom burning. Carbon is burned to nitrogen affecting the transition from oxygen-rich to carbon-rich, although the details depend on the model and particularly on how mass loss is treated. HBB may prevent C stars from forming at all until the envelope mass is depleted (Frost et al. 1998), although in some models C stars do form and HBB turns the stars back from C-rich to O-rich. Another consequence, for stars in a rather narrow mass range, is the formation of lithium and its mixing to the surface via the beryllium transport mechanism (Sackmann & Boothroyd 1992). Smith et al. (1995) surveyed luminous AGB stars in the LMC and SMC for lithium - the clearest indication that HBB is taking place.

Two stars in our sample are known to be undergoing HBB from the presence of lithium in their spectra (Smith et al. 1995), viz, the O-rich star HV 12070 and C-rich star SHV 05210–6904 (SHV F4488). Trams et al. (1999a) showed that the luminous C star, IRAS 04496–6958, has a silicate shell and suggested that it has only very recently stopped HBB. These three stars are of particular interest and are identified in the various diagrams and the discussion which follows.

There are very few LMC stars whose spectra have been examined for lithium, so none of the others are actually known for certain not to be experiencing HBB. Indeed we would guess that several other stars considered in this paper are in the same condition.

A preliminary analysis of the results from this paper, with an emphasis on the significance of the HBB stars in the period-luminosity relation, was discussed by Whitelock & Feast (2000).

## 2 Source Selection and $`JHKL`$ photometry

The coordinates and alternative names of the sources discussed here were listed in Paper I, while Paper II gives the results of the first batch of $`JHKL`$ photometry. In addition to the sources discussed in previous papers in this series we include a few AGB variables previously observed from the South African Astronomical Observatory (SAAO) for which new observations have been obtained. These stars have only thin dust shells, but provide a useful link with the work of Feast et al. (1989) who established a period-luminosity (PL) relation from $`JHK`$ observations of these stars and others like them.

Near-infrared ($`JHKL`$) photometry was carried out between 1993 and 2001 using the Mk III photometer on the SAAO 1.9-m telescope at Sutherland and transformed onto the SAAO standard system defined by Carter (1990). Earlier observations were obtained for some sources as part of other programmes and some of the data published by Feast et al. (1980), Catchpole & Feast (1981) or Glass et al. (1990). The new results are listed in Table 1 which also gives the Julian Date (JD) of the observation. Measurements were made through either 12 or 18 arcsec apertures depending on the seeing and the crowding of the individual source. A few of the older observations were made on the 0.75-m telescope through a 36 arcsec aperture – they are marked with an asterisk in the $`L`$ column. The magnitudes are accurate to better than $`\pm 0.05`$ mag, unless they are quoted to one decimal place when they are good to $`\pm 0.1`$ mag. Many of the sources were too faint to measure at $`J`$, $`H`$, or $`L`$, particularly near minimum light. It also remains possible that some of the fainter magnitudes quoted at $`JH`$ have been slightly contaminated due to the close proximity of other sources.

The $`JHKL`$ observation quoted in Paper II for IRAS 04498–6842 is actually a measure for IRAS 05003–6712. The observation for IRAS 05289–6617 was of the wrong source; only one observation was obtained for the correct object.

The star RHV 0524173–660913 was not part of the ISO programme and was observed towards the end of our programme because it appeared to be a particularly luminous C-star. The C spectral type comes from Reid, Hughes & Glass (1995) and its position in a PL relation is shown in fig. 1 of Groenewegen & Whitelock (1996) where it lies well above the standard relation. However, we identify it here with the star that Westerlund, Olander & Hedin (1981) classified as an M giant and which is known as WOH G311 (or WOH G310). It is considered to be an O-rich star in the following discussion, but another spectrum is required to settle the matter. If it really is an M star that has turned into a C star it is obviously extremely interesting and perhaps similar to IRAS 04496–6958 (Trams et al. 1999a).

## 3 Periods, Amplitudes and Mean Magnitudes

Table 2 is divided into three parts, the first two containing the summary data for the O- and C-rich AGB stars, using the chemical types determined in earlier papers of this series or from Reid et al. (1995). All of the sources in the first two parts of the table are clearly large-amplitude variables with measurable periods. The Fourier-mean magnitudes ($`JHKL`$) are listed, as are the number of $`K`$ observations (n) used to determine the period ($`P_K`$) and the full peak-to-peak amplitude at each wavelength ($`\mathrm{\Delta }J,\mathrm{\Delta }H,\mathrm{\Delta }K,\mathrm{\Delta }L`$). Periods from the literature ($`P_{other}`$) are also listed and the appropriate references given in the notes. There are 5 stars whose mean light is undergoing secular or very long period variations, these are indicated with an asterisk in the $`P_K`$ column and are discussed further in section 4.

The third part of Table 2 contains all of the sources which are not large amplitude AGB variables or have insufficient data to determine periodicities. Seven of these are supergiants; they have $`K<8.0`$ mag, low amplitude variability ($`\mathrm{\Delta }K<0.4`$ mag) and will have bolometric magnitudes brighter than the tip of the AGB, i.e. $`M_{bol}<7.2`$ mag. The other 10 stars are semi-regular or irregular variables or stars without sufficient data to classify them properly (it is possible that the measurements of some of these stars were contaminated by the proximity of other IR sources). The following discussion is limited to regular AGB-variables, i.e. leaving out the supergiants etc.

The pulsation periods ($`P_K`$), listed in Table 2, were determined from Fourier transforms of the $`K`$ light curves. Data from Glass et al. (1990) was used in addition to that given in Table 1, while for stars in common with Wood et al. (1992) and/or Wood (1998) the combined $`K`$ data sets, which provide very well sampled light curves, were used for period determination. The results are generally in good agreement with values determined by others ($`P_{other}`$) except for the semi-regular and irregular variables where the agreement is often rather poor. Among the large amplitude variables, our data for SHV 05260–7011 does not fit the 436 day period determined by Hughes et al. (1990).

In the case of one O- and four C-rich stars the light curves show very long-period or secular changes in addition to the obvious periodic variations from pulsation. The periodicity of these long-term trends is very uncertain and would require many more years of measurement to characterize. The long term trends are discussed further in section 4. The peak-to-peak amplitudes of the best fitting sine curves, given in Table 2, refer only to the periodic pulsations and the full variation of those stars with long-term trends can be much larger than the values tabulated. A few of the very long period O-rich stars have rather asymmetric light curves and the tabulated amplitudes of the best fitting sine curves therefore represent lower limits to the true peak-to-peak variations.

The mean magnitudes are determined from the Fourier fitted sine curves. In the case of stars with a long term trend the mean used is that of the brightest cycles. As these means are used to determine the bolometric magnitudes, this procedure might lead to an overestimate of the total flux, depending on what is causing the long term changes.

Fourier mean intensities are brighter than the mean magnitudes quoted here by up to 0.3 mag at $`K`$, although for most stars the differences are around $`\mathrm{\Delta }K0.1`$ mag. The difference in bolometric magnitudes estimated via intensity means is much smaller (see below).

Figure 3 shows a range of amplitudes at a given period, with the amplitudes of the longer period C stars tending to be larger than those of their O-rich counterparts, possibly because the C-star temperatures are lower and therefore a steeper part of the blackbody light-curve is sampled at these near-infrared wavelengths. The bolometric amplitudes may show much less spread than do these $`\mathrm{\Delta }K`$ values, but our information on them is very limited. At a given period the stars with the larger $`K`$ amplitudes do tend to have redder colours and higher mass-loss rates, as would be expected if pulsation is driving the mass loss; compare, e.g. the two C stars IRAS 05113–6739 and IRAS 05128-6455, both of which have periods around 700 days, but pulsation amplitudes of $`\mathrm{\Delta }K`$ 1.8 and 1.1 mag, colours $`K[12]`$ 7.3 and 5.7 and mass-loss rates $`\dot{M}`$ $`1.1\times 10^5`$ and $`0.6\times 10^5`$ $`M_{}\mathrm{yr}^1`$, respectively (see Table 3). In general the amplitudes of these AGB stars are comparable to those found by Olivier, Whitelock & Marang (2001) among galactic O- and C-rich stars with similar periods. The HBB stars all have rather low amplitudes for their period, although not outside the range of the others.

In a plot of $`K[12]`$ against period or amplitude (see, e.g. Fig. 4 and Fig. 5) there is clear separation of the O- and C-rich stars, although it is interesting to note that the C-rich HBB stars fall with O-rich stars rather than with the other C stars. The separation between C- and O-types is also clear in the near-infrared colours as illustrated in Fig. 6 and Fig. 7.

## 4 long term trends

The $`K`$ light curves which show clear long term trends are illustrated in Figs. 8 to 12. The curves fitted to the data are described in the captions; all of them comprise two periods, the main pulsation period as listed in Table 2 and a significantly longer period chosen to represent the long term trend.

For IRAS 04545–7000, the only O-rich star which shows a trend, a harmonic is also used to provide a much better fit to its distinctly asymmetric light curve. For the four carbon stars a simple sinusoid provides a very good fit to the main pulsational variation.

The time base is far too short to establish the periodicity, if any, of these long term trends. Indeed it is only for TRM 72 and TRM 88 that we have any real indication that the trends might by cyclical, and in these two stars periods of $`5P_K`$ and $`6P_K`$, respectively, were used to fit the trends. However, these long periods are not well constrained and significantly different periods could provide equally good fits to the trends.

The amplitudes of the trends are large, ranging from at least 0.5 up to about 1.0 mag at $`K`$, i.e. comparable to the main pulsation amplitude. The data at other wavelengths are not so extensive, but they are adequate to show that trend amplitudes are wavelength dependent in roughly the same way as the pulsation amplitudes, i.e. greater at shorter wavelengths.

Long term trends and erratic or secular variations have been observed in various long period variables and seem to be more prevalent among C stars than O-rich stars. In the case of V Hya the $`\mathrm{\Delta }K1.7`$ mag variations over a 6160 day cycle (e.g. Olivier et al. 2001 and references therein), are best explained as the result of orbital variations in a binary system. In several other galactic C stars, most notably R For, the long term changes are erratic and best interpreted as the result of mass-loss variations (e.g. Whitelock et al. 1997), possibly akin to RCB variations which are caused by the ejection of “soot” clouds in random directions.

Winters et al. (1994) have predicted periodic mass-loss variations on time scales of a few times the pulsation period. It seems likely that the trends described here are caused by mass-loss variations, but only very long term monitoring will characterize them sufficiently well to see if they fit the Winters et al. model.

## 5 IRAS fluxes

The IRAS fluxes for most of these sources are listed in table 3 of Paper V where a comparison is also made with the fluxes measured by ISO. We use only the 12 and 25 $`\mu `$m flux here, noting the high uncertainty of measurements at longer wavelengths. The fluxes were colour-corrected only for the purpose of estimating the bolometric magnitudes. Thus for the various figures the IRAS magnitudes were calculated from: $`[12]=\mathrm{log}28.3/F_{12}`$ and $`[25]=\mathrm{log}6.73/F_{25}`$, for comparison with other papers in this series, but we note that this process differs from that used by Olivier et al. (2001), for example, who colour-correct the IRAS data used throughout their discussion.

## 6 Bolometric Magnitudes

Bolometric magnitudes, listed in Table 3, were calculated by integrating under a spline curve fitted to the mean $`J_{74},H,K,L,`$ 12- and 25-$`\mu `$m fluxes as a function of frequency, following the procedure described in section 6 of Whitelock et al. (1994).

The $`HKL`$ values used are those listed in Table 2. The $`J`$ mag was converted back to the 74 inch system for consistency with Feast et al. (1989); in fact the $`J`$ flux is negligible for almost all of these stars so the conversion makes no practical difference. The 12- and 25-$`\mu `$m IRAS fluxes (see section 5) were colour-corrected using the 12- to 25-$`\mu `$m flux ratio, as described in the explanatory supplement of the PSC.

In cases where there was no magnitude available for particular wave bands the following procedures were adopted for estimating a value to use in the integration (based on relationships deduced from sources measured in all the relevant bands):
For O-rich stars:
For the one star, SHV 05402–6956, with no $`H`$ magnitude we use equation (3) from Paper IV to derive $`(HK)`$ from $`K[12]`$.
For the eight stars without $`L`$ measures, the expression: $`(KL)=0.239+0.946(HK)`$, which was derived from a least squares fit to the colours of the other O-rich stars, was used to derive $`(KL)`$ from $`(HK)`$. This assumption is important only for SHV 05221–7025 which is moderately red with $`(HK)=1.18`$ and thus $`(KL)=1.36`$; the others all have $`(KL)<0.75`$.
For the four stars without IRAS 12-$`\mu `$m fluxes and the three with upper limits only, we use equation (3) from Paper IV to derive $`K[12]`$ from $`(HK)`$. Stars without 25-$`\mu `$m fluxes are assumed to have uncorrected $`[12][25]=0.9`$ mag.
For the C-rich stars:
For three stars without $`H`$ magnitudes we use equation (2) of Paper IV to derive $`(HK)`$ from $`K[12]`$.
For six stars without $`L`$ mags, $`L`$ is estimated from $`[12]`$ via $`K[12]=1.737+1.838(KL)`$ a relation derived from a least squares fit to the other C stars. For WBP14, $`(KL)=0.4`$ is assumed as its IRAS flux is anomalously faint.
Only one star has no 12-$`\mu `$m flux, IRAS 05291–6700, and this is estimated from equation (2) of Paper IV. The three stars without 25-$`\mu `$m fluxes are assumed to have $`[12][25]=0`$.

In all cases without a $`J`$ measure, O- or C-rich, the relation $`(JH)=1.4(HK)`$ is assumed to derive it.

Because it is only the stars which have most of their flux at long wavelengths that lack $`J`$ or $`H`$ measures, and vice versa, most of these estimates of the missing quantities have little effect on the derived bolometric luminosity. Thus it is only the assumptions about the $`L`$ magnitudes that are potentially important for the bolometric magnitude.

A comparison of Fourier mean magnitudes and intensities for the $`JHKL`$ observations, indicates that the effect of using intensity means rather than magnitude means is to brighten the bolometric magnitude by between 0.01 and 0.06 mag. Using intensities would therefore have a negligible affect on the discussion and conclusions of this paper. Bolometric amplitudes are discussed briefly in section 7.1.

### 6.1 Bolometric Correction

Figure 13 shows the bolometric correction to $`K`$ as a function of various colours. It is illustrated here because it is potentially useful for application to other stars. The $`(KL)`$ and $`(HK)`$ relationships are distinctly different for O- and C-rich stars. For $`K[12]`$, the relationship is very similar for the two types, at least for the redder colours. Note that, because the bolometric magnitude is determined by integrating under the flux curve as a function of frequency, the bolometric correction is not independent of the colours against which it is plotted.

### 6.2 Comparison with Feast et al. (1989)

Feast et al. established PL relations from 29 O-rich and 20 C-rich LMC Miras after deriving bolometric magnitudes by fitting blackbodies to fluxes from mean $`JHK`$ photometry – all that was available at the time.

Fitting blackbodies is an effective way of deriving bolometric magnitudes if the dust shells make a negligible contribution to the total luminosities, as they do for most short period Miras. It is obviously important to look for any systematic differences in the blackbody and spline-fit methods of establishing total luminosity. For this purpose we use both methods to estimate the total flux for sources with $`P<441`$ days and make a comparison in Table 4.

The O-rich stars used by Feast et al. (1989) to establish the PL relation all have $`(JK)<1.5`$ mag. The six O-rich stars in Table 4 with $`(JK)<1.5`$ mag have bolometric and spline luminosities that agree to better than 0.2 mag, five out of six agree to better than 0.04 mag. SHV 05221–7025 shows a slightly brighter $`m_{bol}`$ from the spline fit as would be expected for stronger dust shell contribution; note, however, that SHV 05221–7025 (which is classed among the SRs in Table 2) is the only source in this sample that has an IRAS flux measurement whereas the others were all estimated from the near-IR colours.

The C-rich stars used by Feast et al. (1989) to establish the PL relation all have $`(JK)<3`$ mag. Stars with similar colours in this sample have spline- and blackbody-fit magnitudes that agree to better than 0.2 mag (except for WBP 14, see section 6 and Paper V fig 2 where the ratio of the ISO to revised-IRAS flux for WBP 14 is 2.2 – the largest value found). The redder stars tend to have slightly brighter spline-fit magnitudes. In this case all but one of the stars have measured IRAS fluxes.

In summary, blackbody- and spline-fit bolometric magnitudes agree to better than 0.2 mag (usually better than 0.1 mag) for O-rich stars with $`(JK)<1.5`$ mag and for C-rich stars with $`(JK)<3`$ mag. For stars with larger values of $`(JK)`$ the blackbody method probably underestimates the total flux.

## 7 PL Relation

Figure 14 illustrates the PL relation ($`M_{bol}`$ calculated from IRAS fluxes) for the stars under discussion and compares them with extrapolations of the PL relations derived by Feast et al. (1989) for stars with periods under 420 days. The PL relation derived by Groenewegen & Whitelock (1996) for C stars goes approximately through the middle of the two lines from Feast et al. and would fit these C-star luminosities rather well. WBP 14 is not illustrated in Fig. 14 as there are reasons to think its IRAS flux is too faint. It is, however, illustrated in the IRAS/ISO comparison made below (Fig. 15).

Nishida et al. (2000) discuss the luminosities and periods of three newly discovered C-rich Miras in Magellanic Cloud clusters. The two from the LMC clusters, NGC 1783 and 1978, are illustrated in Fig. 14, where they fall among the other C stars. van Loon et al. (2003) estimate the period and luminosity for the carbon star, LI-LMC 1813, in the cluster KMHK 1603 as (very approximately) 680 days and $`M_{bol}=5.77`$. The turn-off masses for NGC 1783 and 1978 are around 1.5 $`M_{}`$ and that of KMHK 1603 is around 2.2 $`M_{}`$. Taking this into account, and by analogy with Olivier et al. (2001), we can deduce that the carbon stars in this sample will mostly have had progenitors in the 1.5 to 2.5 $`M_{}`$ range.

Feast et al. (1989) pointed out that stars with periods in excess of 420 days lay above the period luminosity relation, while Hughes & Wood (1990) specifically derived a PL relation for stars with periods above 400 days – this is illustrated as a dotted line in the various PL diagrams shown here.

### 7.1 Comparison between bolometric magnitudes from ISO and IRAS data

The differences between the IRAS and ISO fluxes for these LMC sources were discussed in detail in Paper V, where the ISO fluxes were shown to be systematically weaker than those measured by IRAS, although systematic selection effects probably contribute to the differences. The fluxes for these sources are close to the observational limit for the IRAS experiment, while there is some doubt about the calibration of the ISO photometry. Thus it is not clear which measurement should be regarded as more reliable.

A comparison of the long period ($`P>1000`$ days) O-rich sources in the two parts of Fig. 15 shows the bolometric magnitudes from IRAS ($`\overline{M}_{bol}=7.03`$) to be systematically brighter than those from ISO ($`\overline{M}_{bol}=6.76`$). In contrast the IRAS C-star magnitudes are clearly fainter than those derived from ISO data. The difference in the behaviour of the two types of star arises from the application of the colour-correction to the IRAS data. The measured IRAS flux is divided by a factor, $`f`$, which is a function of the ratio of the 12 to 25$`\mu `$m flux. This correction changes the O-rich fluxes very little because the average value of $`f`$ is 1.01 and 1.22 for 12 and 25$`\mu `$m respectively, while for C stars the mean $`f`$ is 1.19 and 1.31 for 12 and 25$`\mu `$m, respectively. It is not possible to derive these corrections with any accuracy given that we have IRAS fluxes only at 12 and 25$`\mu `$m for most of the sources, so this correction procedure is far from ideal.

An interesting feature of the ISO observations relates to the four stars for which measurements were obtained on different dates. The bolometric magnitudes derived from these observations are shown in the upper part of Fig. 15 where observations of the same star are connected with vertical bars. These illustrate that bolometric amplitudes of one magnitude are quite possible (see also Papers IV and VI) and suggest that some of the scatter seen in the various PL diagrams must be attributed to variability.

The greater scatter apparent in the upper (ISO) diagram in comparison to the lower (IRAS) one can be understood in terms of the way the bolometric magnitudes were derived for the two data sets. For the ISO diagram they were derived for the specific epoch at which the particular ISO observation was obtained, while the IRAS diagram makes use of mean magnitudes, both for $`JHKL`$ and from IRAS itself (although the IRAS mean is obviously not a true value as IRAS observed for rather less than a year in total). Fig. 16 illustrates the PL relation derived using the ISO data after making a correction to the mean magnitude (phases calculated from dates given in Paper V). This correction was estimated by first calculating the bolometric amplitude assuming $`\mathrm{\Delta }M_{bol}=\mathrm{log}P2.0`$, then assuming that the bolometric variations are sinusoidal and determining the phase from the $`K`$ light curve. While the corrections are obviously very uncertain, the procedure does result in a very obvious reduction in the scatter over that seen in Fig. 15 (top).

If the process of deriving means for Fig. 16 was perfect then all of the observations for a particular star would coincide in this diagram, and that is obviously not the case. In particular the amplitudes of the HBB stars, which are probably low like their $`K`$ amplitudes (see section 9) are overestimated as can be seen for IRAS 04496–6958.

It will be seen from the above discussion that it is not possible to make a sensible estimate of the uncertainty associated with the bolometric magnitudes. Variability is clearly significant and values derived from single observations may differ from the mean by up to 0.5 mag, but there are also systematic effects, possibly amounting to as much as a few tenths of a magnitude, which depend on the way colour-corrections are treated.

### 7.2 Comparison with Wood et al. (1992) and Wood (1998)

Wood et al. (1992) obtained near-infrared photometry, measured periods and determined bolometric magnitudes for six of the long period O-rich stars discussed here, while Wood (1998) did the same for seven of the C stars and one O-rich star (in fact that O-rich star, IRAS 05329–6708, appears in both papers).

The periods reported in these papers are very similar to ours; we determine better values only by combining Wood’s data with our own, as discussed in section 3. However, the luminosities determined by Wood are distinctly different from ours as can be seen from Fig. 17. Both Wood et al. (1992) and Wood (1998) derived bolometric magnitudes, as we do, from a combination of near-infrared and IRAS fluxes, although the later paper also applies an approximate correction to the IRAS flux to scale it to the mean light. For the O-rich stars the mean from Wood’s data for six stars, $`M_{bol}=6.72\pm 0.17`$, compares to our value of $`M_{bol}=6.91\pm 0.12`$, which is not significant, while for seven C stars the difference is larger with Wood’s data giving $`M_{bol}=5.14\pm 0.17`$ compared to our $`M_{bol}=5.45\pm 0.13`$. For IRAS 05329–6708 there is a difference of 0.43 in the bolometric magnitudes quoted by Wood et al. (1992) and the fainter, phase corrected, value quoted by Wood (1998), both points are illustrated in Fig. 17.

Note that three of the C stars illustrated in Fig. 17 (IRAS 05300–6651, TRM 72, TRM 88) are among those with long term trends discussed in section 4. Nevertheless, the differences between our magnitudes and those of Wood are not systematically influenced, in any significant way, by our approach to establishing the mean near-infrared magnitudes for these stars. In fact our bolometric magnitude ($`5.56`$) for TRM 88, which has long term variations of around $`\mathrm{\Delta }K1.0`$ mag, is slightly fainter than that derived by Wood ($`5.69`$).

### 7.3 PL Relation - Summary

In conclusion we find that, omitting the stars undergoing HBB (see below), both the C- and the O-rich variables tend to follow the trends suggested by an extrapolation of the PL relation determined for shorter period ($`P<420`$ days) stars. There is no clear evidence that these sources with dust shells deviate from the PL relation. They are, however, less luminous than HBB stars with comparable periods. Owing to the uncertainty in the exact mean values of the mid-infrared flux of these stars and therefore in their bolometric magnitudes, it seems inappropriate to attempt to use these data to better define the PL relation.

Wood (1998) particularly emphasizes that the luminosities of the dust enshrouded stars (those with $`500<P<900`$ which we now know are C stars) fall below the extrapolation of the short period PL relation. This is, however, inconsistent with our findings discussed above.

## 8 Mass-loss rates

Mass-loss rates were derived in Paper VI (and are listed in Table 3) for many of the stars under discussion and these are used here to examine the mass-loss dependence on the parameters discussed above. For IRAS 04496–6958 we follow Paper VI in adopting the mass-loss rate derived by assuming the dust is pure carbon.

The dependence of mass-loss rate on luminosity was discussed in Paper VI. Given that these stars obey a PL relation we would anticipate a dependence on period, as is indeed seen in Fig. 18. A very similar dependence is seen for Galactic AGB variables (e.g. Whitelock et al. 1994; Olivier et al. 2001). The relation is the same for O- and C-rich stars, just as the PL relation derived from ISO data (Fig. 15 (top) and Fig. 16) is similar for the two groups. The stars suspected of HBB have lower mass-loss rates than other stars with similar periods. This is of course well known, in that the $`P>420`$ day, luminous, O-rich AGB stars, which we now think are undergoing HBB, were discovered visually, while the C stars of comparable period have thick shells and were only discovered via their infrared luminosity.

Given that pulsation is thought to drive mass loss for these stars, some dependence of mass-loss rates on amplitude is to be expected. However, the pulsation amplitude and period are correlated (Fig. 3), making the effects difficult to separate. Figure 19 shows separate sequences for the O- and C-rich stars which cross around $`\mathrm{\Delta }K1`$ mag and $`\mathrm{log}\dot{M}5.5`$. At large amplitudes the O-rich stars have both longer periods (higher luminosities) and higher mass-loss rates than do their C-rich counterparts. The HBB stars follow the same trend as the others in this figure.

$`K[12]`$ is a measure of the optical depth of the dust shell via the relative emission of the star (which dominates the $`K`$ light) and the shell (which dominates ), and is closely dependent on the mass-loss rate (Whitelock et al. 1994). Figure 20 shows rather different relationships for the O- and C-rich stars in the LMC. Omitting for a moment the clump of four stars at small $`K[12]`$, most of the C- and O-rich stars follow qualitatively similar relations to those found for high mass-loss AGB stars in the solar neighbourhood (Olivier et al. 2001) although the O-rich stars in the LMC reach higher mass-loss rates, $`\dot{M}>10^4M_{}\mathrm{yr}^1`$, at lower $`K[12]`$ than do their Galactic counterparts. Galactic stars with $`P>1000`$ days have much redder colours as has been discussed elsewhere (e.g. Wood et al. 1992).

The group of four stars with low $`K[12]`$ and relatively high $`\dot{M}`$ comprises two of the three HBB stars, HV 2446 which we would suspect is also undergoing HBB and WBP 14 which has too faint an IRAS flux and is probably there erroneously. It is, perhaps, not surprising that the stars undergoing HBB have a different mass-loss vs $`K[12]`$ relation from the other AGB stars.

## 9 HBB Stars

We note that the light-curves for HV 12070, which is known to be HBB, and for HV 2446 and RHV 0524173–660913 (see Fig. 1), which lie in a similar part of the PL diagram (Fig. 15), all have relatively low amplitudes (see Fig. 3). Furthermore, they have similar shapes, exhibiting flat minima (the minima of these particularly bright stars cannot be a consequence of field star contamination).

HV 12070, although O-rich, has an MS spectral type and a rather weak 10$`\mu `$m silicate feature (Paper V). It is therefore likely that it is very close to the O- and C-rich transition, although it is not possible to say in which direction it is evolving.

The three stars thought to be experiencing, or to have recently experienced, HBB (see section 1.1) are more luminous than most sources of comparable period. There are many more O-rich S-type stars that fall in the same part of the PL diagram (Feast et al. 1989; Hughes & Wood 1990), i.e. close to the PL relation derived by Hughes & Wood for stars with $`P>400`$ days. If the relatively high luminosities for these stars are a consequence of HBB (see also Marigo, Girardi & Bressan 1999) then they may also show high lithium abundances and should be examined spectroscopically. It would also be worth examining the O-rich stars with $`P>1000`$ days for lithium.

It is interesting to note that there is a candidate HBB star in IC 1613 which was discovered by Kurtev et al. (2001) and discussed by Whitelock (2002) in the context of AGB stars in Local Group galaxies. It is O-rich, with a spectral type of M3e, has a period of 641 days and is clearly more luminous than an extrapolated PL relation. Its spectrum has not been examined for lithium. Such stars are expected to be a feature of irregular galaxies with significant intermediate age populations.

## 10 Conclusions

Obscured AGB variables with $`P>420`$ days have luminosities close to the extrapolation of the PL relation derived for stars with $`P<420`$ days. We cannot rule out the possibility that there are fainter obscured AGB stars which would not have been found to date because of the limited sensitivity of infrared surveys of the LMC.

We suggest that the apparent change in slope of the PL relation around 400-420 days, for stars with thin dust shells, is caused by the inclusion of variables with luminosities brighter than the predictions of the core-mass luminosity relation, owing to excess flux from HBB.

If RHV 0524173–660913 is confirmed to have changed from an M to a C spectral type (see section 2) then it is of considerable interest. It is presumably undergoing third dredge-up following a helium-shell flash, possibly because hot bottom burning has just terminated.

## Acknowledgments

We are grateful to Robin Catchpole, John Menzies, Ian Glass and Brian Carter who made some of the $`JHKL`$ observations published here.


format: text

text:
# Connectivity of Turing structures

## I Introduction

Nature presents a fascinating diversity of patterns in plants, animals and other natural formations arising often from complex physico-chemical processes Ball (2001). Alan Turing was the first to show that a simple system of coupled reaction-diffusion equations for two chemicals could give rise to spatial patterns due to a mechanism called diffusion-driven instability Turing (1952). These so-called Turing patterns have since been proposed to account for pattern formation in many biological systems, e.g. patterns on fish skin Kondo and Asai (1995); Barrio et al. (1999), butterfly wings Sekimura et al. (2000) and lady beetles Liaw et al. (2002) to mention a few. However, the first experimental evidence of a Turing structure was reported by Castets et al. Castets et al. (1990), who observed a sustained standing nonequilibrium chemical pattern in a single-phase open reactor with chloride-iodide-malonic acid (CIMA) reaction. Currently, there is an increasing interest to develop simple and plausible mathematical models that could describe, at least qualitatively, these pattern formations Nicolis and Prigogine (1977); Haken (1977); Meinhardt (1982); Murray (1993). Although there is a variety of models that could potentially produce similar patterns, a Turing system is perhaps the simplest one.

The forms and variations of patterns generated by Turing systems have been studied by investigating the conditions for instability Szili and Tóth (1993), assuming inhomogeneous diffusion coefficients Barrio et al. (2002), and by introducing domain curvature Varea et al. (1999) and growth Varea et al. (1997). In addition, symmetries in Turing systems are of great interest, since they might have biological relevance, see e.g. Refs. Barrio et al. (2002); Hasslacher et al. (1993); Aragon et al. (2002). Recently, we have studied the effect of dimensionality by simulating three-dimensional Turing systems, which displayed complex pattern formation Leppänen et al. (2002). While in two dimensions one obtains spots, stripes or labyrinthine patterns, in three dimensions complex shapes of, e.g. lamellae, spherical droplets and their combinations appear.

Previously, studies of Turing patterns have typically concentrated on reaction kinetic and stability aspects (see e.g. Judd and Silber (2000); Lengyel and Epstein (1992)), while the issue of pattern structure and its connectivity has received less attention. Here, we focus on connectivity of Turing patterns and its dependence on the parameters of the system. We present a simple way to quantitatively characterize Turing structures and their connectivity. These methods are not only able to characterize the structures, but also to explain some of the effects of nonlinearities behind Turing patterns. We also demonstrate the existence of a “connectivity transition” when the nonlinear interactions of the Turing model are varied.

It was reported earlier that in generic 2D (3D) Turing models nonlinear cubic interactions favor stripes (lamellae) with both morphogens being connected through the system Barrio et al. (1999); Leppänen et al. (2002). On the other hand, quadratic interactions favor spots (spherical structures), in which case only one of the morphogens is connected. Thus, a connectivity transition appears when the nonlinear interactions are varied from being cubic to predominantly quadratic or vice versa. In order to characterize this transition between separated and connected structures, we define a dimensionless control parameter and a scaling function for the number of clusters. This is tested for a variety of system sizes and unstable modes. To characterize the resulting morphologies, we calculated the spatial Fourier spectrum to get the reciprocal space representation of the chemical concentrations and followed its evolution through the connectivity transition.

This paper is organized as follows. Next, we briefly describe the reaction-diffusion model of the Turing kind. Then, we discuss the concept of connectivity in these systems and the methods for characterizing the transition between different patterns. In Sec. IV, we present results of comprehensive numerical simulations, which is followed by Sec. V with a reciprocal space analysis. Then in Sec. VI we draw conclusions.

## II The Model

A Turing system models the evolution of the concentrations of two chemicals, or morphogens, and is in general represented by the following reaction-diffusion equations

$`U_t`$ $`=`$ $`D_U^2U+f(U,V)`$
$`V_t`$ $`=`$ $`D_V^2V+g(U,V),`$ (1)
where $`UU(\stackrel{}{x},t)`$ and $`VV(\stackrel{}{x},t)`$ are the morphogen concentrations, and $`D_U`$ and $`D_V`$ the corresponding diffusion coefficients setting the time scales for diffusion. The reaction kinetics is described by the two nonlinear functions $`f`$ and $`g`$.

In this study we focus on the generic Turing model introduced by Barrio et al. Barrio et al. (1999), in which the reaction kinetics was developed by Taylor expanding the nonlinear functions around a stationary solution $`(U_c,V_c)`$. If terms above the third order are neglected, the system reads as follows

$`u_t`$ $`=`$ $`D\delta ^2u+\alpha u(1r_1v^2)+v(1r_2u)`$
$`v_t`$ $`=`$ $`\delta ^2v+v(\beta +\alpha r_1uv)+u(\gamma +r_2v),`$ (2)
where $`u=UU_c`$ and $`v=VV_c`$ are the concentration fields. The parameters $`r_1`$ and $`r_2`$ set the amplitudes of the nonlinear cubic and quadratic terms, respectively, $`D`$ is the ratio of the diffusion coefficients of the two chemicals, and $`\delta `$ acts as a scaling factor fixing the size of the system. Setting $`D1`$ is a necessary but not a sufficient condition for the diffusion-driven instability to occur. For details about the instability and the linear stability analysis of the model we refer the reader to Barrio et al. Barrio et al. (1999); Leppänen et al. (2002).

From the linear stability analysis Barrio et al. (1999); Leppänen et al. (2002) we obtain the dispersion relation and the conditions for the diffusion-driven instability as the region in $`k`$-space with positive growth rate, i.e., eigenmodes $`u=u_0e^{\lambda t}`$ and $`v=v_0e^{\lambda t}`$ with eigenvalues $`\lambda (k)>0`$. In addition, one can analytically derive the modulus of the critical wave vector

$$k_c^2=\frac{1}{\delta }\sqrt{\frac{\alpha (\beta +1)}{D}},$$
(3)
which was here determined for the case $`\alpha =\gamma `$ (set to obtain only one stable state at $`u=v=0`$ in the absence of diffusion). In a discretized three-dimensional cubic system, the wave number is of the form

$$|\stackrel{}{k}|=\frac{2\pi }{L}\sqrt{n_x^2+n_y^2+n_z^2},$$
(4)
where $`L`$ is the system size and $`n_x`$, $`n_y`$, $`n_z`$ are the wave number indices (in a two-dimensional system $`n_z=0`$). By adjusting the parameters and allowing only a few unstable modes, one can obtain several different parameter sets. As in our earlier work Leppänen et al. (2002) we chose the parameters $`D=0.516`$, $`\alpha =\gamma =0.899`$, $`\beta =0.91`$ and $`\delta =2`$ corresponding to a critical wave vector $`k_c=0.45`$, and $`D=0.122`$, $`\alpha =0.398`$, $`\beta =0.4`$ and $`\delta =2`$ corresponding to $`k_c=0.84`$.

In this study, we vary the parameters $`r_1`$ and $`r_2`$ in Eq. (2), since they control the appearance of stripes or spots. By gradually changing these parameters we observe a transition from spotty (2D) or spherical droplet (3D) patterns to striped (2D) or lamellar (3D) patterns. In order to investigate this transition, numerical simulations were carried out by discretizing the spatial dimensions into a square or cubic cell lattice and calculating the Laplacians in Eq. (2). In all of our simulations we used $`dx=dy=dz=1.0`$. The equations of motion were iterated in time using the Euler scheme with time step $`dt=0.05`$. The boundary conditions were chosen to be periodic, and initially both chemicals were distributed randomly over the whole system.

## III Connectivity

In the numerical simulations of Eq. (2) one deals with two concentration fields with characteristic wave lengths. In order to visualize this, the concentration of only one of the chemicals is typically plotted with a gray scale, since in this type of a system the fields are in anti-phase, i.e., if there is a large amount of chemical $`U`$ in some sub-domain, the concentration of chemical $`V`$ would be low there. These concentration fields vary continuously having diffuse boundaries. What do we mean by connectivity in patterns of chemicals?

To answer this question one can define sub-domains dominated by either chemical $`U`$ or $`V`$, provided that the amplitudes of the patterns is large enough. If we define the boundary as the interface between sub-domains dominated by different chemicals, we can easily locate the boundaries, since the concentrations change rapidly, typically within one or two lattice sites. Now, if two points belong to the same domain, i.e., are not separated by a boundary, they are considered connected. The definition of the boundaries in this way is conceptual in the sense that in the $`U`$-dominated domains the concentration of $`V`$ does not have to be zero, only much less than the concentration of $`U`$.

In Fig. 1 we show changes in the concentration fields, i.e., $`u`$ and $`v`$, of a 2D system for different values of nonlinear parameters $`r_1`$ and $`r_2`$ in Eq. (2). When the cubic term ($`r_1`$) dominates, the resulting stationary pattern is striped with a small number of imperfections, see Fig. 1A. These imperfections can be considered as topological defects, or dislocations, which could serve as nucleation sites for spots. More dislocations appear (see Figs 1B-C) when the strength of the quadratic term is made larger. As the quadratic term increases, more spots nucleate and they arrange themselves to triangular structure and at the same time getting rid of the remaining stripes (see Figs. 1F-H). Finally when the cubic term is diminished even further, only spots remain. In this sequence of simulations the strength of the cubic term was changed relative to the quadratic term by using a single control parameter $`P`$, which we elaborate below. Nevertheless, the transition from striped to spotty pattern seems to happen quite abruptly in $`P`$. In the last frame of this sequence (Fig. 1I) we see a fully stabilized spotty pattern with almost perfect triangular symmetry in two different orientations such that there is a mirror plane or twin boundary between them.

Now, let us discuss patterns in Fig. 1 from the clustering point of view. In order to simplify this, and without loss of generality, we can assign zeros and ones to the whole lattice based on the chemical which dominates a given domain. With this mapping we can consider the number of clusters, which can be calculated using the well-known Hoshen-Kopelman algorithm Hoshen and Kopelman (1976) as in typical percolation problems. However, before applying this cluster algorithm let us first visually inspect possible shapes of extended clusters (Fig. 1). In Fig. 1A one can see that in the case of stripes the number of $`U`$\- and $`V`$-dominated clusters is almost the same, and both types are extended dominantly in one of the dimensions. On the other hand, in the case of a spotty structure (Fig. 1I), chemical $`U`$ appears as separate round clusters or spots, whereas chemical $`V`$ forms only one connected cluster. Between these two limiting cases there is the transition region, depicted in Figs. 1D-F, where $`U`$-dominated clusters appear as spots and stripes in the form of a “string-of-pearls”. A question arises whether the transition from stripes to spots could be characterized in a reasonable way by the number of clusters? With the Hoshen-Kopelman algorithm we can directly determine the number of clusters ($`N`$) in the system, but we also need to define a single control parameter that $`N`$ is a function of.

One should bear in mind that many different parameter pairs of nonlinear terms ($`r_1`$ and $`r_2`$) result in similar patterns. To obtain some general insight into the dynamics, we apply dimensional analysis Barenblatt (1987) and derive a dimensionless control parameter for describing this behavior. The relevant dimensions are as follows

$$[\alpha ]=\frac{1}{s},[r_1]=\frac{1}{[c]^2},[r_2]=\frac{1}{[c]s},$$
(5)
where $`s`$ denotes seconds and $`[c]`$ is an arbitrary unit of concentration. Thus, the control parameter $`P`$ can be written in the following dimensionless form

$$P=r_1(\frac{\alpha }{r_2})^2.$$
(6)
The two limits of this parameter, i.e., $`P\mathrm{}`$ ($`r_20`$) and $`P0`$ ($`r_10`$) yield stripes and spots (Figs 1A and 1I), respectively, as the two extremes in 2D. The same holds for 3D. Between these two limits a certain combination of striped (lamellar) and spotty (spherical droplet) patterns is expected to coexist (see Figs 1D- 1F).

## IV Simulation Results

In order to study the connectivity in two- and three-dimensional Turing patterns, we have performed extensive simulations using system sizes up to $`5\times 10^5`$ lattice cells and up to $`2\times 10^6`$ time steps to reach a stationary state. For all cases the results are taken as statistical averages of at least 20 separate runs. Using $`P`$ (Eq. (6)) as a control parameter is plausible, since the transition from a striped (lamellar) pattern to a spotty (spherical droplet) pattern took place in a very narrow region of $`P`$ regardless of the individual values of the nonlinear parameters ($`r_1`$ and $`r_2`$), the unstable mode (dictated by $`\alpha `$ and $`\beta `$) or the system size. This is clearly seen in Fig. 2, where the number of $`U`$\- and $`V`$-dominated clusters is shown as a function of the control parameter $`P`$ for two different cases in 2D. In our studies the $`P`$-space was scanned by keeping either $`r_1`$ or $`r_2`$ constant and by carrying out separate simulations for each value of $`P`$. From Fig. 2 it is clearly seen that the mode transition from spots to striped structure is quite sharp in $`P`$. The control parameter $`P`$ serves the purpose of a unique transition variable, having a well-defined value for the transition to occur.

The most unstable mode for the smaller ($`100\times 100`$) system is with $`k_c=0.84`$, while for the larger ($`200\times 200`$) system it is with $`k_c=0.45`$. The larger system has smaller wave vector and thus the wave length, i.e., characteristic length of the pattern is larger. This is why the values on the vertical $`N`$-axis, are of the same order of magnitude. In order to compare the numbers of clusters one can normalize it by dividing with $`N_c^d`$, where $`N_c=k_cL/2\pi `$, $`L`$ denoting the linear system size (square or cube) and $`d`$ the spatial dimension. $`N_c^d`$ is the maximum number of spherically symmetric clusters in a $`d`$-dimensional system if the clusters were uniformly distributed and the effect of boundaries was neglected. Due to the periodicity of the chemical structure, the number of clusters in the actual $`d`$-dimensional system can be estimated to be $`N^d=(N_c+1/2)^d`$. However, an additional correction is required to take into account the effect of boundaries. One can estimate the number of additional partial clusters due to boundaries by estimating the length (area) of the boundary and the number of clusters within this domain ($`dN^{d1}`$). As a result of this discussion we propose the normalization function for the number of clusters to be

$$C_d(N(P),N_c)=\frac{N(P)}{N_c^d}\left(1\frac{d}{N_c+\frac{1}{2}}\right),$$
(7)
where $`N(P)`$ is the calculated number of clusters for control parameter $`P`$. By revising the Hoshen-Kopelman algorithm one could have directly calculated the number of clusters by taking periodicity into account, in which case Eq. (7) reduces to $`C_d(N(P),N_c)=N(P)/N_c^d`$. However, this approach was not implemented and the normalization was carried out by using the Eq. (7).

We have also studied the effect of the system size on the mode transition in 2D. This is shown in Fig. 3, where the normalized quantity $`C_2(N(P),N_c)`$ is plotted for $`U`$-clusters against $`P`$ for three different system sizes ($`L=`$ 100, 175 and 500). Neglecting the number of $`V`$-clusters does not affect our conclusions, since the curves would be symmetrical as can be seen from Fig. 2. On the other hand, one can clearly see in Fig. 3 that the control parameter $`P`$ succeeds in capturing the essential features of the transition. In addition, it can be seen that the normalization function of Eq. (7) scales the number of clusters such that it collapses onto the same curve with only very small deviations outside the transition.

If, on the other hand, one carries out the simulations for very small systems, finite-size effects can be observed. For small system sizes the $`C_2(N(P),N_c)`$ curve becomes very steep in the transition region. This would suggest that in the limit of small systems, the transition would become almost discontinuous. However, the system cannot be made infinitely small since the (periodic) boundary conditions start to affect the behavior of the system. As discussed earlier the spots tend to nucleate from topological defects, or dislocations, of the striped pattern, i.e., from the points where the stripes coincide (Fig. 1). In the case of a small system even one dislocation can affect the morphology of the whole system and thus quickly transform stripes into a lattice of spots. In a larger system many dislocations have to appear at various sites to give rise to spots which in turn make the appearance of more spots favorable.

So far we have discussed our simulation results in 2D systems. We have also studied the connectivity transition extensively in three dimensions. In this case stripes and spots become lamellae and spherical droplets, respectively, and the structure seems more complicated especially in the transition region. For illustrations of three-dimensional Turing structures we refer the reader to Ref. Leppänen et al. (2002). In 3D the transition does not occur at the same point with respect to $`P`$ as in 2D since the third dimension gives to the clustering process one more degree of freedom, and thus it is easier for the structures to connect.

This is indeed what one finds. Figure 4 depicts the normalized number of clusters for three different system sizes ($`L=`$ 40, 50 and 75). One can see that the behavior of the system is different from two dimensions. Now, the transition occurs at a value of $`P`$ which is a decade smaller than in 2D, since a smaller cubic nonlinear coefficient ($`r_1`$) favoring lamellar structures is sufficient for increasing connectivity in three-dimensional space. The significance of this large drop of the critical $`P`$-value remains unanswered. In addition, unlike in 2D we did not observe any finite-size effects for the smallest possible system sizes.

## V Reciprocal space analysis

Next, we exploit another way to study the connectivity transition, namely the discrete Fourier transform of the concentration data in the wave vector space, $`\stackrel{}{k}`$, i.e.,

$$\widehat{\rho }(\stackrel{}{k})=\rho (\stackrel{}{r})e^{i\stackrel{}{k}\stackrel{}{r}}.$$
(8)
This approach has been used earlier, e.g. in connection of reaction-diffusion systems Karttunen et al. (1998), Turing patterns Yang et al. (2002), and to characterize the evolution of patterns Karttunen et al. (1999). The quantity $`\widehat{\rho }(\stackrel{}{k})`$ corresponds to a diffraction pattern. Figure 5 shows a sequence of the original concentration fields and their diffraction patterns in a two-dimensional system for different control or transition field parameters, $`P`$.

In Fig. 5A, the diffraction intensity is predominantly to $`k_x`$-direction due to more stripes in y-direction, while the smaller diffraction intensity around the $`k_y`$-axis is due to the stripes in x-direction. The distance of these diffraction peaks from the origin gives the length of the wave vector of the unstable mode, while its width in the perpendicular direction describes deviations of stripes from the principal directions due to stripe tilting and bending. In Fig. 5B the diffraction intensity becomes more spread around the $`(k_x,k_y)`$ diagonal, as a result of appearing dislocations and nucleating spots. After that in Figs. 5C-D the diffraction intensity starts splitting into separate peaks due to more spots forming, then developing into six separate equidistant (from the origin and from each other) intensity peaks as evident in Figs. 5E-F. This hexagonal symmetry in the reciprocal space is because the system evolves towards regular spotty pattern with predominantly triangular symmetry in the real space. However, these diffraction peaks are somewhat spread and in fact each peak split in two, which indicates that the triangular symmetry is not perfect over the whole system.

If one examines the real space picture of spots in Fig. 5F, one might think that there are fairly large disordered regions that deviate from the ideal triangular perfect lattice of spots. However, the corresponding reciprocal space picture reveals that this is not so, since it does not show diffuse diffraction pattern or a ring due to orientational disorder, but it does show a double peak structure as can be seen in crystalline materials due to twinning or low-angle grain boundaries. The existence of this kind of boundaries is common place in any nucleation mechanism and usually it is associated to the fact that twinning or low-angle grain boundaries require only small local displacements at low energy cost. In the present case, however, there is no energetics involved and in Fig. 5F it is difficult to find clear twinning, which we found in Fig. 1I. There are other more general ways of establishing order, namely the appearance of eutactic local structures defined in the study of regular polytopes Coxeter (1973). It has been recently proposed Torres (2002) that eutacticity is a very important property exhibited not only by crystalline and quasicrystalline lattices, but also by the geometrical forms of some biological systems. The presence of order in these Turing patterns, as revealed by their clear diffraction patterns, suggests that the principles governing the preference for eutactic structures also apply to the present case.

Here, we have done the reciprocal space analysis for 2D systems, but the same analysis can easily be extended to 3D systems by fixing the orientation of one of the spatial vectors, thus obtaining an in-plane diffraction graph. These diffraction patterns, however, are expected to look more or less the same as in 2D, which can be understood on the basis that almost every cross-section of a lamellae or spherical droplet pattern in 3D would look as striped or spotty pattern, respectively.

## VI Discussion

In this study, we have investigated the connectivity of spatial patterns generated by the reaction-diffusion mechanism both in 2D and 3D. This was done by cluster analysis for the dominating chemical. Since several different combinations of parameters in a Turing system can produce similar patterns, we derived a dimensionless control parameter to investigate the effects of nonlinearities. The numerical simulations were consistent with the predictions drawn from the dimensional analysis, and the system showed a transition in the proximity of the same $`P`$-value irrespective of the individual system parameters. With help of an analytically derived normalization function, the number of clusters collapsed onto the same master curve independently of the system size or the unstable mode.

From biological perspective, the characteristics of the Turing system presented here could be important. The fact that the transition from striped or lamellar to spotty or spherical droplet structures does not depend on the size of the system or the specific parameters involved, but just the combination of the nonlinearities ($`P`$), makes the mechanism very general. In addition, the observation that irregular structures are not as stable as regular structures could be used to explain the steepness of the transition: Unfavorable structures corresponding to the transition domain occupy minimal volume of the phase space (possible structures).

The presence of eutacticity is well known in atomic crystal structures. The fact that similar geometric principles seem to be acting also in the formation of Turing patterns could be an important feature in favor of their application in the study of the geometrical forms of some simple living organisms. This could also support the idea that Nature prefers geometrical forms generated from eutactic stars Torres (2002). This very simple analysis of Turing patterns and their morphological transitions due to nonlinearities could open up a new way of investigating universal features of patterns obtained by very different mechanisms.

The approach we have taken here to study pattern formation shares common features with percolation and could be studied to certain extend as such. Percolation of morphogens to different directions is not as interesting a measure as the number of clusters, since it holds only a binary value and does not capture features of the morphological structure effectively. However, if one considers the field of biological applications, where Turing systems are often used, percolation behavior is interesting and has been used for studying some biological problems Sahimi (1994).

If one considers the patterns in Fig. 1 one can observe that in the first frame both morphogens have percolated in the vertical direction. On the other side of the transition one of the morphogens has not percolated while the other has done so in both directions. Should one of the chemicals be favorable for diffusion, these two opposite stages obtained by a small variation of nonlinear parameters could for example correspond to chemical concentrations in tissue as signaling via diffusion is either enabled or disabled. However, one should always be careful in making suggestions for biological applications due to complexities of nature not yet understood.

###### Acknowledgements.

We wish to thank János Kertész for helpful discussions, and one of us (R. A. B.) wishes to thank the Laboratory of Computational Engineering at Helsinki University of Technology for its hospitality. This work has been supported by the Academy of Finland through its Centre of Excellence Program (T. L. and K. K.).