format: text

text:
# Tight LMC massive star clusters R 127 and R 128Based on observations obtained at the European Southern Observatory, La Silla, Chile

## 1 Introduction

R 127 and R 128 (Feast et al., 1960) are two Large Magellanic Cloud (LMC) massive stars situated some 300 pc south of the famous 30 Doradus star forming factory and 150 pc west of the H ii region N158 (Henize, 1956). They are in fact the brightest members of the two adjacent tight clusters forming NGC 2055 towards the center of the OB association LH 94 (Lucke & Hodge, 1970). Although they lie in the H ii region DEM 248 (Davies et al., 1976), they are not apparently linked to any specific nebulosity.

R 127, also called Sk–69$^∘$ 220 (Sanduleak, 1970) or HDE 269858, is recognized as a Luminous Blue Variable (LBV), a very rare class of evolved massive stars. Parker (1997) listed 8 members in the LMC, while van Genderen (2001) increased the number of members and candidates to 21. In the past they were called S Dor variables, a designation introduced by Kukarkin et al. (1974), and there is a strong tendency by some workers to use the original name. These most luminous stars (log $`L/L_{}`$$``$ 5.0–6.3) are characterized by irregular photometric and spectral variations over decades and evolve from a hot (OB-type) visual minimum phase to a cooler (A-type), visual maximum (Humphreys & Davidson, 1994). In the Hertzprung-Russell (H-R) diagram they are located very close to the observed upper luminosity boundary for very massive stars, the Humphreys-Davidson limit. LBVs are characterized by extreme instability, dramatic outbursts and high mass loss. Between violent eruptions LBVs still lose mass at high rates, 10-5 to 10$`^4`$$`M_{}`$ yr-1 (Lamers et al. (2001) and references therein). The material ejected (up to several $`M_{}`$) forms a nebula around LBVs with typical sizes of 0.5–2 pc (Nota & Clampin, 1997). The properties of LBVs are reviewed by Humphreys & Davidson (1994) and more recently in the proceedings of a special workshop edited by Nota & Lamers (1997). LBVs are believed to be the precursors of Wolf-Rayet stars; however their precise evolutionary state is still poorly understood. R 127 is quite noteworthy even among this small population; it is the brightest LBV star in the LMC ($`M_{bol}=10.5`$ mag) and has the largest magnitude variation ($`\mathrm{\Delta }`$V $``$ 2.5 mag) (Stahl et al., 1983; Humphreys & Davidson, 1994; van Genderen, 2001). Due to its exceptional characteristics, R 127 has been the subject of much research in the past. R 127 was classified as OIafpe extr or WN 9-10 by (Walborn, 1977, 1982). A brightening of 0.75 mag later caused its classification as an S Dor variable (Stahl et al., 1983). It was then included in the Long-Term Photometry of Variables organized by Sterken (1983), and published in several papers (Manfroid et al., 1991, 1994; Sterken et al., 1993, 1995; van Genderen et al., 1997). The monitoring showed that R 127 had a very bright phase around the end of 1986 with $`y=9.05`$ mag, and as such was the brightest star in the LMC (Wolf et al., 1988), only surpassed in February 1987 by SN1987A. R 127 reached its visual maximum of $`V`$$``$ 8.8 mag around the year 1989 and then slowly decreased to $`V`$$``$ 9.4 mag in $``$ 1995. A new maximum of $`V`$$``$ 9.3 mag was attained in 1997 coinciding with the time of present observations. Subsequently it became fainter, with $`V`$$``$10.4 mag, around the year 2000 (Stahl 2002, private communication).

Infra-red observations have shown that R 127 is surrounded by a dust shell (Stahl et al., 1984). As this dust shell is believed to be a crucial element for understanding the evolution of massive stars, it has been extensively studied using imagery (Stahl, 1987), coronography (Clampin et al., 1993), spectroscopy (Smith et al., 1998), and polarimetry (Schulte-Ladbeck et al., 1993).

R 128, otherwise Sk–69$^∘$ 221 or HDE 269859, is a supergiant B2 Ia (Fitzpatrick, 1991). It is variable with a total range of 0.32 mag between 1983 and 1990, which is very large for its spectral type (van Genderen et al., 1998). Interestingly, this star has formerly been considered as an LBV candidate (van Genderen, 2001).

The studies so far devoted to these two interesting bright stars have mainly dealt with their individual characteristics. However, recent findings, from both high-resolution observations and theoretical works, suggest that massive stars form in groups. Therefore, knowing the characteristics of the cluster members is necessary for better understanding the formation and evolution of these brightest, probably the most massive, stars of the group. Nota et al. (1991) looked into the multiplicity of R 127 by obtaining high-resolution, ground-based images with the STScI coronograph and reported the presence of 20 stars towards R 127. Although these components are mainly field stars detached from R 127 proper, they detect a relatively close component lying 3′′.5 north-west of the bulk of R 127.

The present work is therefore devoted to the photometry and spectroscopy of the cluster members. Using high resolution imaging techniques, we aim at resolving the cores of the clusters R 127 and R 128. Moreover, high-spatial-resolution spectroscopy will allow us to study the physical properties of so far unknown members of these two tight clusters.

## 2 Observations and data reduction

### 2.1 Sub-arcsecond imaging and deconvolution photometry

The R 127 and R 128 clusters were observed on 20 November 1997 using the ESO New Technology Telescope (NTT) equipped with the active optics SUperb Seeing Imager (SUSI). The detector was a Tektronix CCD (#42) with $`1024\times 1024`$ pixels of 24 $`\mu `$m (0′′.13 on the sky), and the seeing varied between 0′′.81 and 1′′.34 (fwhm).

The observations were performed in the $`uvby`$ Strömgren photometric system using the ESO filters # 715, 716, 713, and 714 respectively. We were particularly careful to keep most of the brightest stars in the field under the detector’s saturation level to have at our disposal high quality Point Spread Function (PSF) stars. This led us to adopt exposure times of 350, 175, 210 and 140 seconds in $`u`$, $`v`$, $`b`$ and $`y`$ respectively. We also used ditherings of 5′′ –10′′ for bad pixel rejection and in order to be able to use the full oversampling capabilities of the MCS deconvolution algorithm (Magain et al., 1998). Indeed when performing simultaneous deconvolution of several frames, the algorithm uses the different frame centerings as a constraint while decreasing the pixel size. We took a grid of seven dithering positions for each filter. Luckily, the targets of interest are close enough to be contained in a single SUSI field of view. Unfortunately, the $`u`$ images could not be used for the photometry due to their insufficient quality.

Photometry was derived in the Strömgren $`v`$, $`b`$ and $`y`$ filters according to the following procedure: after bias subtraction and flat-fielding, the seven frames were co-added in each of the filters. The photometry of the stars lying outside the compact clusters was performed on the resulting frames through the DAOPHOT reduction package. This yielded the photometry of 233 stars situated outside clusters R 127 and R 128 (Fig. 1).

The crowded clusters were processed with the MCS restoration algorithm. The deconvolution was performed on a 256 $`\times `$ 256 pixel region containing both R 127 and R 128. The MCS code, proposed and implemented by Magain et al. (1998), results from a new approach to deconvolution taking care not to violate the Shannon (1949) sampling theorem: the images are deliberately not deconvolved with the observed PSF, but with a narrower function, chosen so that the final deconvolved image can be properly sampled, whatever sampling step is adopted to represent the final data. For this purpose, one chooses the final, well-sampled PSF of the deconvolved image and computes the PSF which should be used to perform the deconvolution. The observed PSF is constructed from several stars close enough to the clusters in order to avoid any possible PSF variation across the field.

The restoration resolves the R 127 and R 128 clusters into 14 and 33 components respectively (Fig. 2). The deconvolution method requires a high S/N ratio for the sources in order to achieve accurate photometry and astrometry. Therefore, very faint stars too close to the bright components were excluded from the deconvolution process. Star #13 was detected but its magnitude could not be accurately measured, due to the presence of relatively strong process residuals.

A technical problem during the observing run prevented us from obtaining adequate standard star observations necessary for calibrating the photometry. The available calibration frames allowed us however to correct for the atmospheric extinction for all filters and derive the zero point for the $`y`$ filter. Results concerning $`y`$ magnitudes were checked using the photometric data available in the literature for three stars in the field, Sk–69 217, Sk-69 218 and R 128 (Ardeberg et al., 1972; Isserstedt, 1975). Intrinsic $`(by)_0`$ colors were calculated for each star of known spectral type (see Sect. 4), in two steps: first we established their effective temperatures from the calibration by Vacca et al. (1996), then the corresponding theoretical $`(by)_0`$ colors were deduced using Relyea & Kurucz (1978). An average interstellar reddening of $`E(BV)=0.15`$ mag (St-Louis et al., 1997), or $`E(by)=0.12`$ mag (Kaltcheva & Georgiev, 1992), was used to obtain the observational $`(by)`$ colors.

The final photometric results for the two clusters R 127, R 128, and the brightest field stars are presented in Tables 1, 2, and 3 which also list the $`\alpha `$-$`\delta `$ (2000.0) positions of the stars. The astrometry was obtained from the identification of several bright, isolated stars with GSC 2.2 stars. The positional accuracy is found to be around 2′′.

A word of caution seems appropriate regarding the colors. The powerful deconvolution method has allowed us to resolve the compact clusters R 127 and R 128 revealing their so far unknown components. Furthermore, the code has enabled us to perform the photometry of the tight components. However, we should underline that this photometry is relative for a number of reasons which have nothing to do with the limitations of the code. The shortcomings in the observation of standard photometric stars particularly affect the colors, and we have therefore taken care not to over-interpret them.

### 2.2 Spectroscopy with NTT/EMMI

The EMMI spectrograph attached to the ESO NTT telescope was used on 22 November 1997 (BLMD mode) to obtain several long slit spectra. The grating was # 12 centered on 4350 Å and the detector was a Tektronix CCD (# 31) with 10242 pixels of size 24 $`\mu `$m. The range was 3810–4740 Å and the dispersion 38 Å mm-1, giving fwhm resolutions of $`2.70\pm 0.10`$ pixels or $`2.48\pm 0.13`$ Å for a 1′′.0 slit. At each position we first took a short 5 min exposure followed by one or two longer 15 min exposures. The instrument response was derived thanks to observation of the calibration stars LTT 1020, LTT 1788, EG 21.

The seeing varied from 0′′.45 to 0′′.53, allowing us to obtain un-contaminated spectra. The identification of the stars along the slits was based on sketches drawn during the observations. Since the targets lie in crowded regions, some ambiguous identifications required the development of a small IRAF task, using the position angle information in the FITS headers. First, each spectrum was integrated along the dispersion axis, the result being stored into a two-pixel wide strip, which is close to the actual size of the slit. Then, the position angle and the pixel-arcsec correspondence were used to calculate the rotation matrix for the World Coordinate System. This allowed creation of a slit chart, an $`\alpha `$-$`\delta `$ calibrated two-dimensional image containing accurate slit orientations. Displaying simultaneously the slit chart besides the SUSI images and using the WCS correlations it was possible to accurately check the identity of the star in the slit.

## 3 Photometry results

### 3.1 Components of the clusters

Fig. 1 presents the two tight clusters R 127 and R 128 and their immediately surrounding field (Sect. 3.2). The central image restored by deconvolution is displayed in Fig. 2, while Tables 1, 2 and 3 list the results of the photometry.

We show that R 127 is composed of at least 14 stars. The brightest component, #7 with $`y=9.37`$ and $`(by)=0.23`$ mag, must be the real LBV (Sect. 4.1). The closest components to #7, stars #8 and #3 having $`y=14.39`$ and 12.12 mag respectively, lie at 1′′.5 and 3′′.3 from the main component. Star #3 turns out to be very interesting as shown by its spectrum (Fig. 7). R 127 was the target of coronographic imaging through a broad-band $`R`$ filter and also a narrow-band H$`\alpha `$ +N ii filter by Clampin et al. (1993) and Nota et al. (1991). They detect star #3 which they call R 127B. Those observations also show the presence of a faint bipolar nebulosity that fails to appear on our frames. On the other hand, their use of a 2′′.7-diameter masking disk centered on R 127 forbids any detection of the companion identified as star #8. Even star #10 remains undetected, possibly due to low S/N ratio. The existence of stars #8 and #10 is however confirmed by an unpublished $`R`$ frame which we obtained in 1991 with NTT+SUSI under better seeing conditions. More recently, R 127 has been observed using WFPC2 on HST (archive data, GO 6540, P.I. Schulte-Ladbeck). Again, stars #8 and #10 are visible, as well as several fainter spots, but the combination of the diffraction pattern, saturation by star #7, and the presence of nebulosity makes these spots difficult to identify.

R 128 is resolved into at least 33 stars. It is therefore more populated than its neighboring cluster R 127. The main component, #37 with $`y=10.82`$ mag, should be considered as the star usually referred to as R 128, since the next brightest star of the cluster, #24 with $`y=13.16`$ mag, is much weaker.

Fig. 8 shows the color-magnitude (C-M) diagrams of the R 127 and R 128 clusters. The position of the brightest star of the sample, #7 ($`y=9.37`$, $`(by)=0.23`$ mag), is confirmed by independent photometric results. In fact the light curve of R 127 from 1985 to 2002 (Stahl, private communication) shows that the LBV attained a brightness of $`V=9.35\pm 0.2`$ at the end of the year 1997. Similarly, the correlation between the magnitude and color of R 127 (Spoon et al., 1994) indicates a color of $`(by)=0.23\pm 0.05`$ mag corresponding to a visual brightness of $`y=9.36`$ mag.

The C-M diagram shows a main sequence centered on $`(by)`$$``$ 0.00 mag for both clusters. The intrinsic colors of the O type stars vary from –0.12 to –0.15 mag for $`T_{eff}`$ ranging from 30 000 to 50 000 degrees and log g = 3.50 to 5.00 (Relyea & Kurucz, 1978). Taking (b – y)0 = –0.14 mag as the average value and using the conversion relation $`E(by)=0.67E(BV)+0.02`$ (Kaltcheva & Georgiev, 1992), we find $`E(BV)=0.18`$ mag ($`A_V=0.6`$ mag) for the average reddening of clusters R 127/R 128. The reddest components of the R 127 cluster are #8 and #9, which we will discuss in Sect. 5. The R 128 cluster also shows a number of red stars.

The absolute magnitudes, $`M_V`$, listed in Tables 1, 2, and 3 were derived assuming equality between the Strömgren $`y`$ and Johnson $`V`$ magnitudes. Taking a distance modulus of 18.6 mag for the LMC (Whitelock et al., 1997), we used the above-mentioned relation to correct for the reddening of each star. The corresponding luminosities used in Sect. 5 were obtained from the bolometric corrections given by Vacca et al. (1996) for O type stars. There is now general agreement that Vacca et al. (1996) overestimate the temperatures and luminosities (Martins et al., 2002). However, their calibration is adequate for our purpose in this paper, especially since we use their bolometric corrections. For B types and the LBV R 127, which shows a (peculiar) A type supergiant spectrum, we used the calibration by Fitzpatrick & Garmany (1990).

### 3.2 Field stars

Photometry was obtained for 244 field stars lying in the direction of the R 127 and R 128 clusters. Among them those brighter than $`y=16`$ mag are listed in Table 3. Stars #283 and #284, situated out of our images, appear in this table because they happened to fall on the spectrograph slit and therefore have spectral classifications. The brightest stars of this sample are #61, #70 (more commonly known as Sk $`69^{}`$ 217 and Sk $`69^{}`$ 218) and then #290 and #113. The C-M diagram for the whole sample of the field stars is displayed in Fig. 8, and their detailed photometry is available upon request.

## 4 Spectral types

The spectral classification of the OB stars was performed with respect to the digital atlas of Walborn & Fitzpatrick (1990). No classification standards with the current observational setup were available, so there may be a slightly greater uncertainty in the (lower) luminosity classes than is usual in this work, but it does not exceed one class and may well be smaller.

### 4.1 R 127 cluster

Spectrograms of the two brightest stars, #3 and #7, are available. These are the “preceding” (western) and “following” (eastern) components of R127 (HDE 269858), respectively, as denoted by Feast et al. (1960), and we confirm that the eastern component, #7, is the peculiar star. The two components had similar apparent magnitudes when observed by Feast et al. (1960), but as discussed in the Introduction, #7 is currently in an extended LBV phase and is much brighter. Our observation of its spectrum is illustrated in Fig. 3. It is representative of the peculiar A supergiant LBV maximum phase, dominated by emission lines of hydrogen and singly ionized heavy metals (e.g., Fe ii, \[Fe ii\], Ti ii, Cr ii); see Wolf et al. (1988) for a detailed presentation of this kind of spectrum at high resolution.

We believe that our spectrogram of star #3 (Fig. 7) is the first spatially resolved observation of its spectrum. It has an interesting, very well defined B0 Ia type with deficient nitrogen (cf., e.g., the extreme weakness of the $`\lambda 4097`$ absorption in the blue wing of H$`\delta `$, which is near maximum strength in a normal spectrum of this type). The latter characteristic indicates that this star either is a relatively un-evolved supergiant, or was a slow rotator on the main sequence (e.g., Walborn et al. (2000) and references therein).

### 4.2 R 128 cluster

We have spectroscopy for 11 stars in this compact cluster, which is shown in Figs. 4–7. Most of them are late-O or early-B giants (i.e., the main sequence was not reached spectroscopically), but one (#23, Fig. 5) is classified O9 Ib, and #37 is R 128 itself, the brightest star in the cluster, which we classify B1.5 Ia (N weak). R 128 was classified B1 Ia by Fitzpatrick (1988) and B2 Ia (N weak) by Fitzpatrick (1991). We prefer the intermediate type on the basis of the strength of Si iv$`\lambda 4089`$ and concur on the anomalous deficiency of nitrogen; e.g., N ii$`\lambda 3995`$ reaches its maximum strength at B2 Ia and would be comparable to He i $`\lambda 4026`$ in a normal spectrum. Current explanations of such anomalies are as cited above for star #3 in the R 127 cluster.

### 4.3 Field stars

Finally, 6 stars in the surrounding field were also observed spectroscopically. The most interesting are the relatively bright stars #61 (Sk $`69^{}\mathrm{\hspace{0.17em}217}`$) and #70 (Sk $`69^{}\mathrm{\hspace{0.17em}218}`$) shown in Fig. 4. The Of nature of Sk $`69^{}217`$ was discovered by Walborn et al. (1991) and the spectral classification of O7 Iaf from the present, higher quality material is in perfect agreement with the earlier result; this is the hottest star in the present sample. On the other hand, the earlier classification of ON9 Ib for Sk $`69^{}\mathrm{\hspace{0.17em}218}`$ is here revised to O8 Iab(f), with no nitrogen/carbon anomaly. This discrepancy can be understood in terms of the lower quality of the earlier spectrogram, which led to a later type at which the N/C line strengths would be anomalous, but they are not at the earlier type derived here. The recently identified S iv$`\lambda \lambda \mathrm{\hspace{0.17em}4486}`$, 4504 emission lines (Werner & Rauch, 2001) are prominent in both of these spectra.

## 5 Discussion

The brightest members of the R 127/R 128 clusters are massive evolved stars as indicated by their spectral types. There are, however, stars residing still on the main sequence (Fig. 8), but they have been missed in our spectroscopy due to their relatively low brightness. In order to look into the evolutionary states of the clusters’ stars, we have used the theoretical models of Meynet et al. (1994) for a metallicity of z = 0.004 and high mass loss rates. Fig. 9 presents the relevant isochrones and evolutionary tracks on which are overlaid the stellar positions. The filled dots are based on the physical parameters given by Vacca et al. (1996)’s calibration for the derived spectral types. The crosses indicate the corresponding luminosities obtained from our photometry. We see that for a number of stars the observed luminosities are much smaller than those predicted by Vacca et al. (1996)’s calibration. We will briefly discuss this problem below.

The positions of the filled dots suggest that the hottest stars of the sample, #61 and #70, are $``$ 3 Myr old, and have evolved from initial masses of $``$ 80 $`M_{}`$ into late-type O supergiants. Likewise, the oldest stars, #39 and #35, are B0 supergiants with an age of $``$ 6 Myr and a ZAMS mass of $``$ 25 $`M_{}`$. It seems therefore that in the R 127/R 128 region lower mass stars have formed prior to massive ones. Moreover, the positions of the crosses, taken at face value, suggest masses as low as $``$ 15 $`M_{}`$ with ages up to $``$ 8 Myr.

The LBV member of the sample, star #7, does not appear in the H-R diagram since its relatively low effective temperature puts it outside the plot. It is well-known that LBV stars change their spectral type, temperature and radius, becoming cooler and larger during the visual maximum phase. However, their luminosity remains approximately constant, with log $`L/L_{}`$ = 6.1 for R 127, as reported by several workers (Stahl et al., 1983; Wolf, 1989; Lamers et al., 1998; van Genderen, 2001). An effective temperature of log $`T_{eff}`$ = 3.954 is calculated by Lamers et al. (1998) for R 127 corresponding to its phase of visual maximum.

Star #7 is the most evolved member of the R 127/R 128 clusters, presumably because it was the most massive star of the group. An initial mass of $``$ 85 $`M_{}`$ can be attributed to the LBV progenitor from Fig. 9. A present mass of 46$`{}_{3}{}^{}{}_{}{}^{+17}`$ $`M_{}`$ is derived for the LBV by Lamers et al. (1998). This means that star #7 has lost $``$ 45% of its mass since its formation $``$ 3 Myr ago, corresponding to a constant mass loss rate of $`\dot{M}`$ = 1.3 $`\times `$ 10-5 $`M_{}`$ yr-1. However, mass loss is not a regular process and LBVs undergo drastic episodes during which they lose mass dramatically. For example, model calculations by Langer et al. (1994) show that during the first 1.5 Myr of life of a 60 $`M_{}`$ star the mass loss rate is constant, $``$ 5 $`\times `$ 10-6 $`M_{}`$ yr-1. Then it becomes much stronger, $``$ 3 $`\times `$ 10-5 $`M_{}`$ yr-1, during the next 2 Myr. The peak of mass loss during the LBV phase reaches $``$ 3 $`\times `$ 10-3 $`M_{}`$ yr-1 when the star has an age of $``$ 3.4 Myr. On the other hand, association in a close binary system has been suggested for accelerating the evolution of massive stars and creating the LBV phenomenon (Tutukov & Yungelson (1980); Gallagher (1989); Humphreys & Davidson (1994)). However, so far no observational support exists for this LBV being a binary system.

Fig. 9 shows a discrepancy between the luminosities suggested by the spectral type calibration schemes (Vacca et al., 1996) and those derived from the photometry. While the smaller deviations, like those of stars #61, #70, #3, #45, etc., can be attributed to measurement uncertainties, a discrepancy as large as that for star #23 called for verification. This trend was first interpreted as due to a sort of bias. However, despite our efforts, no systematic errors were found in our photometry. In fact, our examination of several published papers on other regions showed that $`\mathrm{\Delta }M_V`$ (the difference between the absolute magnitudes derived from spectral types and those from photometry) increases with magnitude, although with different extents (Fitzpatrick, 1988; Walborn & Blades, 1997; Walborn et al., 1999; Parker et al., 1992; Massey & Johnson, 1993). Such trends are expected from the defining relations, if there are dispersions in the actual absolute magnitudes or distances, or systematic effects in the spectral classification, absolute-magnitude calibration, or photometry. Uncorrected differential extinction or relative errors in the calibration will broaden the distribution about the trends.

Aperture photometry observations previously performed on R 127 used diaphragms of 10 ′′ to 15′′ in diameter. As a result, the global magnitude of a small group of stars was mistaken for that of a single source. Assuming that we can simulate such a measurement by co-adding the stars found inside a circle of 10 ′′ centered on star #7 (which are stars #8, #3, #10, #9 and #6), we find that the measured $`y`$ magnitude would be 9.27 instead of 9.37. This small difference shows that the contribution of the weak companions identified in this work is not significant for the luminosity of the main object.

As for the brightest member of the R 128 cluster, star #37, it is a supergiant B1.5 Ia (N wk) with log $`L/L_{}`$ = 5.76 and log $`T_{eff}`$ = 4.24 (van Genderen, 2001) and therefore lies outside our H-R diagram (Fig. 9). This star is known to be variable and has even previously been dubbed as an LBV candidate (van Genderen, 2001). Indeed an initial mass of $``$ 50 $`M_{}`$, as derived from Fig. 9, does not rule out the possibility for star #37 to become an LBV (Maeder, 1997; Langer et al., 1994).

As mentioned in Sect. 3.1, some of the cluster components show rather red colors. In particular, star #8, the closest neighbor of the LBV, lying 1′′.5 away, has the reddest color. This reminds one of the LMC transition star R 84 which has a close evolved component of type M2 Ia lying $``$ 1′′.7 from it (Heydari-Malayeri et al., 1997). Note, however, that R 84 is a lower mass object belonging to a B-type cluster. Anyhow, presently we cannot determine whether star #8 is evolved, and the red color of this star may have a totally different origin. We are aware of the color accuracies involved (Sect. 3.1), nevertheless if star #8 is physically associated with the cluster, and not a sight-line coincidence, its color could be due to contamination by the ejecta from the LBV. Clampin et al. (1993) measure a size of 8′′.0 $`\times `$ 9′′.0 (1.9 $`\times `$ 2.2 pc) for the circumstellar nebula lying at a position angle of 165. Even a smaller nebula, as found by Stahl (1987), with a size of 3′′.2 $`\times `$ 4′′.4 (0.8 $`\times `$ 1.1 pc), can affect the color measurements of star #8. A larger nebula would also explain the red color of star #9.

## 6 Concluding remarks

The decomposition of the LMC compact star clusters R 127 and R 128 into at least 14 and 33 components respectively and medium-dispersion spectroscopy carried out in very good seeing conditions have allowed us to study the stellar contents of these interesting objects. This work is essential since the R 127 cluster harbors the most prominent LBV object in the LMC. We resolve the core of R 127 into four components and clearly identify the LBV object, which shows spectral features typical of a visual maximum phase. Moreover, we present the first spatially resolved spectrum of the brightest neighbor of the LBV, star #3 lying at a separation of 3′′.3. It turns out to be a supergiant B0 Ia (N wk).

We also show that the star currently known as R 128 is in fact the component number #37 of an adjacent cluster. A good quality spectrum of this supergiant leads to a revised classification of B1.5 Ia (N wk).

The two clusters are composed of evolved massive stars. The oldest members are $``$ 6–8 Myr old and the most massive one has an initial ZAMS mass of $``$ 80 $`M_{}`$. These age estimates are in agreement with the fact that no H ii regions are associated with the clusters. The most massive stars have had enough time to disrupt the molecular cloud and the nebula, and the earliest spectral types are no longer present.

This spectroscopic study misses the main sequence stars of the clusters, which probably have lower initial masses with respect to the evolved ones. It will therefore be interesting to observe them spectroscopically in order to determine their real status. It will also be attractive to study the dynamics of the clusters by obtaining accurate radial velocities of the members.

###### Acknowledgements.

We are particularly indebted to Dr. Pierre Magain, Institut d’Astrophysique et de Géophysique de l’Université de Liège, Belgium, and his group for their warm hospitality offered to F.M. during his stay in Liège for learning to use the MCS deconvolution algorithm. Without their patience, availability, and help the training could not be fruitful. Our thanks go also to Dr. Frédéric Courbin for discussions and advice at several opportunities during his trips to Paris. We are also grateful to the referee, Dr. A.M. van Genderen, for his helpful comments. We would like also to thank Dr. Otmar Stahl, Landessternwarte, Heidelberg, for discussions and information about the monitoring of R 127 and Dr. Jesús Maíz-Apellániz, STScI, about the absolute-magnitude trends.


text:
# I Introduction

## I Introduction

Nonlinear optical fibers and waveguides and arrays composed of them furnish a unique example of a medium in which various solitary-wave patterns and their complexes can be easily observed in a real experiment, and described, with a very high accuracy, by relatively simple models, a benchmark example being the nonlinear Schrödinger equation . Besides commonly known bright and dark solitons, solitary-wave structures in the form of domain walls (DWs) were also predicted in a fiber with normal group-velocity dispersion (GVD) which carries two different waves with orthogonal polarizations, circular or linear, that interact through the cross-phase modulation (XPM) induced by the Kerr nonlinearity . Similar structures can also be expected to exist in planar nonlinear optical waveguides . DW patterns are distinguished by the property that, asymptotically, they contain a single polarization, with a switch between two of them in a localized region. In fact, solutions for the optical DWs were constructed, following the pattern of earlier known solutions of the DW type for a system of coupled Ginzburg-Landau equations that describe interactions between roll patterns with different orientations in a convection layer .

GVD in the fiber must be normal in order to prevent the modulational instability (MI) of the DW’s uniform background fields. Nevertheless, it is known that a uniform two-component field, unlike single-component ones, may be subject to MI even in the case of normal GVD . Loosely, DWs are related to the MI of the two-component uniform state the same way as the usual bright soliton is related to MI of the single-component field in the case of anomalous GVD .

The same model (which is presented in detail below) applies to an altogether different physical system, namely, a string of drops of a binary (two-component) Bose-Einstein condensate (BEC) trapped at minima of a periodic potential, which can be readily induced by an optical interference pattern . In this connection, it is relevant to mention that stable DW configurations have been predicted in a continuous quasi-one-dimensional (cigar-shaped) binary BEC .

Optical DWs in nonlinear fibers have been observed in direct experiments , including high-repetition periodic DW trains . On the other hand, recent experimental achievements in the observation of discrete spatial optical solitons in arrays of waveguides in the spatial domain  suggest that observation of DW-like structures in waveguide arrays may be quite feasible too. Additionally, solitons of the DW type may be a new species of solitary waves in the discrete nonlinear Schrödinger (DNLS) equations, which have recently attracted a great deal of interest (for a recent review see e.g., Ref. ). Given that the solutions to coupled DNLS equations have been examined for a considerable while now (the first relevant results appeared about 20 years ago in ), it appears that DWs may be one of the few types of DNLS solitons that have not been studied yet.

An objective of the present work is to introduce this type of discrete solitons and study their stability by means of precise numerical methods. An exact formulation of the model, together with estimates of relevant physical parameters, are given in section 2. In section 3, we study in detail DWs in the case of strong XPM, which corresponds to circular polarizations. We find several different types of DWs, the simplest one being stable for all values of the inter-site coupling constant $`C`$ (which is the ratio of the propagation length determined by the Kerr nonlinearity to the linear-coupling length, in terms of the optical-waveguide array). DWs of a different type exhibit a rather unusual stability behavior, being stable at small values of $`C`$, and stable at larger $`C`$, i.e., in weakly and strongly coupled arrays, respectively. The existence of stable bound states of DWs is also demonstrated. Other types of DWs turn out to be either completely unstable, or stable only at small values of $`C`$. In section 4, we consider DWs in the model with weak XPM, which corresponds to linear polarizations. In this case, DWs are unstable. However, the instability growth rate of the simplest (fundamental) DW becomes vanishingly small for large values of $`C`$, so that the pattern becomes (marginally) stable in the continuum limit, $`C\mathrm{}`$. In all the cases, the predicted stability of DWs is tested in direct numerical simulations, and in those cases when DWs are expected to be unstable, the instability development is simulated too. Section 5 summarizes the paper.

## II Formulation of the model

The model of an array of nonlinear optical fibers carrying fields $`\varphi _n(z)`$ and $`\psi _n(z)`$, which correspond to two orthogonal polarizations of light, has the form

$`i{\displaystyle \frac{d}{dz}}\left(\psi _n\right)`$ $`=`$ $`C\left(\psi _{n+1}+\psi _{n1}2\psi _n\right)(|\psi _n|^2+\beta |\varphi _n|^2)\psi _n,`$ (1)
$`i{\displaystyle \frac{d}{dz}}\left(\varphi _n\right)`$ $`=`$ $`C\left(\varphi _{n+1}+\varphi _{n1}2\varphi _n\right)(|\varphi _n|^2+\beta |\psi _n|^2)\varphi _n,`$ (2)
where $`z`$ is the propagation distance along the fiber and $`n`$ is the index of the lattice site. Equations (1) and (2) are written in a rescaled form, in which the constant $`C`$ of the linear coupling between adjacent fibers has a straightforward physical meaning: it is a ratio of the characteristic propagation length $`L_{\mathrm{nonlin}}`$ along the waveguide, determined by the Kerr nonlinearity (self-phase modulation, SPM), to the coupling length $`L_{\mathrm{coupl}}`$ determined by the linear interaction between adjacent waveguides. In the most typical experimental situation, with the power of light beams $`5001000`$ W, the nonlinearity length in AlGaAs waveguides (whose nonlinearity is $`500`$ time stronger than that of fused silica glass) takes values $`12`$ mm, and the coupling length is $`1`$ mm . In fact, $`L_{\mathrm{coupl}}`$ may vary in broad limits, as it exponentially depends on the separation $`h`$ between the waveguides: for instance, $`L_{\mathrm{coupl}}`$ decreases by a factor of $`1.6`$ as $`h`$ increases from $`9`$ $`\mu `$m to $`11`$ $`\mu `$m . Results displayed below (such as a width of DWs in terms of the number of lattice sites) suggest that stable DWs in optical-waveguide arrays should be observable in the same experimental setups, and in essentially the same range of the power-per-waveguide, where bright discrete solitons have been found in works .

In Eqs. (1) and (2), the SPM coefficient is normalized to be $`1`$, and $`\beta `$ is a relative coefficient of the XPM nonlinear coupling, which is $`2`$ or $`2/3`$, in the case of linear or circular polarizations, respectively. In the latter case, we neglect, as usual, four-wave-mixing (FWM) nonlinear coupling terms  Strictly speaking, for the short propagation distance ($`10`$ mm), relevant to the recent experiments , it may be necessary to keep the FWM terms in the case $`\beta =2/3`$; we do not consider this issue in the present work because, as it will be shown below, most interesting results are found in the case $`\beta =2`$. We will refer to the two cases with $`\beta =2/3`$ and $`\beta =2`$ as those with the weak and strong XPM coupling, respectively. The value $`\beta =2`$ applies also to the case when the two modes refer not to polarizations, but rather to light signals carried by different wavelengths , which is another plausible realization of the present model in terms of nonlinear optics.

If the nonlinearity in the waveguides is induced by the usual self-focusing Kerr effect, and the light is launched so that its Poynting vector is oriented parallel to the array, the linear coupling between fibers (discrete diffraction) corresponds to $`C<0`$ in Eqs. (1) and (2). In this case, the system gives rise to two-component discrete bright solitons, which were recently studied in detail, including their generalization to the two-dimensional lattice . However, a newly developed experimental technique, based on launching an oblique beam into the array , makes it possible to implement anomalous discrete diffraction, which corresponds to $`C>0`$ in Eqs. (1) and (2). In this case, bright solitons do not exist.

As it was mentioned above, the same system of Eqs. (1) and (2) may also be regarded as a normalized model of a string of binary-BEC drops trapped in an optical lattice. In this case, assuming the usual situation with the positive scattering length (i.e., repulsion between atoms), one has $`C>0`$, which is precisely what is necessary to generate DW solutions. Unlike the realization in terms of nonlinear optics, the value of the coefficient $`\beta `$ may be arbitrary (but positive). The separation between the drops is of the order of the lattice-generating light wavelength, i.e., $`1`$ $`\mu `$m. A promising candidate for the binary condensate being a mixture of 85Rb and 87Rb, the corresponding positive scattering lengths being $`10`$ nm . The necessary temperature and densities of the drops can be estimated as $`10^510^6`$ K and $`10^{11}10^{12}`$ cm-3, which can be readily achieved by available experimental techniques

In this work, we examine discrete domain walls (DWs) in the system (1) and (2), which may be stable only if $`C>0`$. We will construct DW solutions, starting from the anti-continuum (AC) limit with $`C=0`$ , and using the continuation in $`C`$ to extend the solutions up to bifurcation points where they lose their stability, or up to the continuum limit ($`C\mathrm{}`$) when possible. Thus, we seek solutions of the form

$`\psi _n`$ $`=`$ $`\mathrm{exp}(i\mathrm{\Lambda }z)u_n,`$ (3)
$`\varphi _n`$ $`=`$ $`\mathrm{exp}(i\mathrm{\Lambda }z)v_n,`$ (4)
arriving at the stationary equations

$`F(u_n,v_n)C\mathrm{\Delta }_2u_n+(|u_n|^2+\beta |v_n|^2)u_n\mathrm{\Lambda }u_n`$ $`=`$ $`0,`$ (5)
$`G(u_n,v_n)C\mathrm{\Delta }_2v_n+(|v_n|^2+\beta |u_n|^2)v_n\mathrm{\Lambda }v_n`$ $`=`$ $`0.`$ (6)
Once a solution to Eqs. (5) and (6) has been found (by means of a Newton-type numerical method), we perform the linear stability analysis around it, looking for perturbed solutions as

$`\psi _n`$ $`=`$ $`\mathrm{exp}(i\mathrm{\Lambda }z)\left[u_n+ϵa_n\mathrm{exp}(i\omega z)+ϵb_n\mathrm{exp}(i\omega ^{}z)\right],`$ (7)
$`\varphi _n`$ $`=`$ $`\mathrm{exp}(i\mathrm{\Lambda }z)\left[v_n+ϵc_n\mathrm{exp}(i\omega z)+d_n\mathrm{exp}(i\omega ^{}z)\right]`$ (8)
(the asterisk stands for the complex conjugation), and solving the ensuing matrix eigenvalue problem.

## III Numerical results for the strong-XPM model

In this section, we consider the model based on Eqs. (1) and (2) with $`\beta =2`$. We start the examination of DW structures by looking for stationary solutions in a natural form which is taken, in the AC limit, as

$$u_n=(\mathrm{}0,0,0,V_1,1,1,1,\mathrm{}),v_n=(\mathrm{}1,1,1,V_2,0,0,0,\mathrm{})$$
(9)
\[see Eqs. (3) and (4)\], where $`V_1`$ and $`V_2`$ belong to one (central) site of the lattice and will be defined below. We will, in particular, consider the steady states with $`\mathrm{\Lambda }=1`$ in Eqs. (3) and (4), but results will be generally true if $`1`$ in AC ansatz (9) is replaced by $`\sqrt{\mathrm{\Lambda }}`$ for an arbitrary positive value of $`\mathrm{\Lambda }`$. Notice that in the present study, we fix $`\mathrm{\Lambda }`$ and vary $`C`$, modifying, essentially, in this way the degree of localization or equivalently the peak power. One can instead always rescale $`C`$ to the value $`C=1`$, and equivalently vary the propagation constant $`\mathrm{\Lambda }`$.

Numerical calculations demonstrate that, for $`\beta =2`$ (which corresponds to the circular polarizations of light in the fibers), the AC pattern (9) with

$$V_1=1,V_2=0$$
(10)
(or vice versa) generates the most structurally robust and stable solutions. In particular, this solution was found to be stable for all values of the coupling constant $`C`$, up to the continuum limit $`C\mathrm{}`$. Examples of such a DW for cases of weak ($`C=0.034`$) and strong ($`C=3.5`$) inter-site coupling (recall $`C`$ may be realized as the ratio $`L_{\mathrm{nonlin}}/L_{\mathrm{coupl}}`$ in the optical-waveguide array) are shown in Fig. 1. The complete stability of DWs belonging to this branch of the solutions has also been verified by direct simulations of the full nonlinear equations (1) and (2), with initial conditions containing a small noisy component.

The next case to consider is when

$$V_1=V_2=\sqrt{\mathrm{\Lambda }/(1+\beta )}$$
(11)
in Eq. (9). These values imply that, in the AC limit, the central waveguide in the array carries a vectorial state, in which both polarities have equal amplitudes. The numerical investigation reveals an unusual feature of this solution branch: it is unstable for small values of $`C`$, starting from the AC limit ($`C=0`$), but becomes stable at $`C=0.61`$, and remains stable thereafter up to $`C=\mathrm{}`$. This property is opposite to the common scenario, when a sufficiently strong discreteness is responsible (through the effective potential energy barrier that it creates) for the stabilization of various solitary-wave lattice patterns (for instance, two-dimensional pulses  or vortices ) that are unstable in the continuum limit.

Solutions generated by Eq. (11) are shown, for the same values of $`\beta `$, $`\mathrm{\Lambda }`$, and $`C`$ as in Fig. 1, in two upper rows of Fig. 2 (top panel). Notice the presence of an unstable (imaginary) eigenvalue pair in the panel pertaining to $`C=0.034`$. The imaginary part of the unstable eigenvalue is displayed, as a function of the coupling constant, in the lower panel of Fig. 2, showing the transition from instability to stability at $`C0.61`$.

The predictions for the stability of these solutions were also checked against direct simulations of the full equations. In the case when the DW is predicted to be stable, it is indeed found to be completely stable (not shown here). In the case when it is expected to be unstable, the simulations show (see Fig. 3) that the unstable DW emits a packet of lattice phonons and rearranges itself so that the field in one component at one site ($`n=50`$) in Fig. 3) becomes equal to the field in the other component at the adjacent site ($`n=51`$, in Fig. 3). Comparison with the expressions (9) and (10) clearly suggests that the result of the instability development is the rearrangement of the DW into a stable one belonging to the solution branch generated by Eq. (10) in the AC limit.

A solution generated by the AC pattern (9) with $`V_1=V_2=0`$ (i.e., with both field components equal to zero at the central lattice site in the AC limit) was also examined. Unlike the solution branches considered above, this one does not reach the continuum limit. Instead, it terminates at $`C=0.15143`$ through a turning-point bifurcation. This is obvious in Fig. 4, which displays the $`L^2`$ norm of one component of the solution (the other component behaves similarly),

$$P\left(\underset{\mathrm{}}{\overset{+\mathrm{}}{}}u_n^2\right)^{1/2},$$
(12)
vs. the coupling constant $`C`$ (note that $`P`$ is finite due to the finiteness of the computational domain). The dependence $`P(C)`$, shown in Fig. 4, is not a completely invariant characteristic, as values of $`P`$ depend on the size of the domain. However, the turning-point structure is invariant.

In the lower panel of Fig. 4, the present solution and its linear stability eigenvalues are shown for three different values of $`C`$ ($`C=0.06`$, $`C=0.11`$, and $`C=0.15143`$ in the in the top, middle, and bottom subplots). As is seen from this part of the figure, the first solution is stable, while the latter two are not. The solutions of the present type are stable for $`C<0.0725`$, and the spectrum of their linear-stability eigenvalues contains two continuous bands. As $`C`$ increases, two pairs of eigenvalues bifurcate from the outer band and move towards the inner one. Their first collision with the inner band occurs at $`C0.0725`$, and the second collision takes place at $`C0.08`$. Each collision generates an instability-bearing quartet of eigenvalues, manifesting the so-called Hamiltonian Hopf bifurcation . Subsequently (for larger $`C`$), the two quartets move towards the imaginary eigenvalue axis. Very close to the turning point, the collision of the first quartet with the axis induces a symmetry breaking, which results in one pair of eigenvalues moving towards the origin, while another pair moves upwards along the imaginary axis. Finally, as a result of the collapse of the lower imaginary pair onto the origin of the spectral plane, the turning-point bifurcation occurs and the branch terminates, as is shown in the upper panel of Fig. 4.

Direct simulations again show that the solutions which are stable according to the linearization are indeed stable in full simulations. An example of direct simulations of the unstable solution belonging to this branch is given in Fig. (5). As is seen, the instability does set in, but its growth is extremely slow.

We also considered the situation with the zero-field state occupying, in the AC limit, two (or more) lattice sites. The resulting scenario turns out to be similar to that demonstrated above for the patterns generated by the AC limit (9) with $`V_1=V_2=0`$: the solution is stable at very small $`C`$, then becomes unstable due to the bifurcation of four (instead of two in the previous case) pairs of eigenvalues from the outer band and their collision with the inner band (not shown here). Finally, the branch terminates at a turning point at $`C0.17`$, similarly to what was shown in Fig. 4. Note that this turning point is found at a larger value of $`C`$ than in Fig. 4. We continued the analysis, starting with several empty (zero-field) sites in the center of the pattern at $`C=0`$, which produced quite a similar picture, an unexpected feature being that the value of $`C`$ at the turning point increases with the increase of the number of the initially empty sites.

Contrary to what was described above for the case of one or several empty sites, in the case when the initial pattern has two sites occupied by the vectorial state (11), the solution behaves quite differently from its counterpart in the case when only one site was initially occupied by the vectorial state. In fact, not only does the branch terminate - in this case, at $`C0.1085`$ (the solution found at the turning point and its stability eigenvalues are shown in Fig. 6) - but it is also found to be always unstable: for every value of $`C`$ starting from the AC limit $`C=0`$, there are two imaginary eigenvalue pairs. The branch terminates when one of them passes through the origin.

An example of directly simulated evolution of an unstable DW of this type is shown in Fig. 7. In this case, the instability development is essentially faster than in the case shown in Fig. 3, and, as a result, a stable DW appears.

Furthermore, stable DW solutions found above, such as the ones generated by the AC-limit pattern (9) with $`V_1`$ and $`V_2`$ chosen as per Eqs. (10) or (11), can form stable complexes (higher-order DWs). Obviously, the complex must contain an odd number of fundamental DWs. An example is shown in Fig. 8 for $`C=0.15`$. Direct simulations confirm that these complexes are completely stable.

Recall that, if the underlying system of Eqs. (1) and (2) is applied to the string of BEC drops, the value of $`\beta `$, which is determined by the values of the cross-scattering length, is arbitrary. In the case $`\beta >1`$, the results are quite similar to those described in detail above for $`\beta =2`$. For $`\beta <1`$, the results are very different, as is described in the next section.

## IV The weak-XPM model

The underlying pattern (9), used to construct all the DW states considered above, was suggested by solutions found (in the temporal, rather than spatial, domain) in the continuum model of the single nonlinear optical fiber  or cigar-shaped BEC . Solutions of this type in the continuum model exist only in the case of sufficiently strong XPM, namely, for $`\beta >1`$.

Unlike the continuum model, in the discrete system (1), (2) DW patterns of the type (9) can also be found for $`\beta <1`$ (weak XPM coupling), including the case $`\beta =2/3`$, which is specially relevant for the applications to nonlinear optics. However, in this case properties of the DW solutions are drastically different from those described above for $`\beta =2`$.

For $`\beta =2/3`$, all solutions of the DW type were observed to become staggered when continued from the AC limit (i.e., the phase difference between the fields at adjacent sites of the lattice tends to be $`\pi `$, rather than $`0`$). As a result, the outer continuous band of the stability eigenvalues has the Krein signature  opposite to that of the inner band, and as soon as the bands collide (which happens at $`C0.011`$, i.e., still in the case of very weak coupling), numerous quartets of complex eigenvalues arise. An example, corresponding to the solution generated by the AC pattern (9) with the central site occupied by the vectorial state (11) (with $`\beta =2/3`$) is displayed in Fig. 9. Similar results were obtained for the solution branches generated by the pattern (9) with $`V_1=1`$ and $`V_2=0`$ (also with $`\beta =2/3`$). Thus, discrete DW solutions of this type may exist in the weak-XPM case, unlike the continuum limit, but they are stable only at very small values of the coupling constant. These conclusions, concerning the existence and (in)stability of the DW patterns for the weak-XPM case, are corroborated by direct simulations; however, the instability is extremely slow, therefore it is not shown here.

In the regime of weak XPM, another type of DW solutions is suggested by analogy to the findings of Ref.  for the continuum version of the model. This type of solution exists only for the weak-XPM case ($`\beta <1`$) in the continuum model, and, quite naturally, turns out to be primarily relevant in the same case in the discrete model. In the continuum limit, this solution is

$$u(x)=\sqrt{\mathrm{\Lambda }}\mathrm{cos}(\chi ),v(x)=\sqrt{\mathrm{\Lambda }}\mathrm{sin}(\chi ),$$
(13)
where $`\mathrm{\Lambda }`$ is the same as in Eqs. (3) and (4), and, in the first approximation (which corresponds to a broad DW),

$$\chi (x)=\frac{\pi }{4}\mathrm{tan}^1\left[\mathrm{exp}(s\sqrt{ϵ\mathrm{\Lambda }}x)\right],ϵ1\beta .$$
(14)
The corresponding AC limit of such solutions is (for $`\mathrm{\Lambda }=1`$)

$$u_n=(1/\sqrt{2},\mathrm{},1/\sqrt{2},1/\sqrt{2},\mathrm{},1/\sqrt{2}),v_n=(1/\sqrt{2},\mathrm{},1/\sqrt{2},1/\sqrt{2},\mathrm{},1/\sqrt{2}).$$
(15)
Starting from Eq. (15), the solution was extended all the way from $`C=0`$ to the continuum limit, which shows that the solution exists for all values of $`C`$. The upper panel of Fig. 10 shows two examples of this solution, for $`C=0.034`$ and $`C=3`$. The lower panel of the figure clearly shows that this DW state is unstable at all values of $`C`$, but gets stabilized in the continuum limit, $`C\mathrm{}`$. For instance, for $`C=3`$ the growth rate of the relevant instability is already extremely small, $`10^6`$, hence this state may seem a practically stable one in the discrete case too, provided that the coupling constant is large enough.

These predictions for the DW patterns generated by the ansatz (15) were verified in direct simulations of the full nonlinear system with $`\beta =2/3`$. Figure 11 shows the development of the relatively strong instability in the case $`C=0.1`$. A remarkable result, which makes this case drastically different from those considered above, is that the instability transforms the quiescent DW into a moving one, which is accompanied by emission of quasi-linear lattice waves.

It has also been verified that, in full accordance with the prediction of the linear-stability analysis, the formally unstable DWs of the present type seems to be virtually stable ones. A typical example of that is displayed in Fig. 12.

Even though the continuum-limit analytical expressions (13) - (14) for these solutions are only valid for $`\beta <1`$, one can seek for similar solutions in the discrete strong-XPM system, with $`\beta =2`$. The same way as solutions of the type (9) can exist in the discrete system with $`\beta =2/3`$, despite the fact that they do not exist in the continuum limit for the weak-XPM case (see above), the solutions of the type (13) - (14) have a chance to exist in the discrete strong-XPM model. We have found that such solutions indeed exist for very small values of $`C`$, but they exhibit very strong instabilities, with a part of the imaginary eigenvalue axis being populated by the continuous spectrum. In Fig. 13, a solution of this type for $`\beta =2`$ and $`C=0.12`$ is compared with a “natural” one existing at $`\beta =2/3`$ and $`C=0.2`$.

The nonlinear evolution of this strong instability was directly simulated, showing a quick transition to a state of “lattice turbulence”, see Fig. 14. Thus, the lattice admits the existence of the “unnatural” solutions in both strong- and weak-XPM cases, but it never allows them to persist ad infinitum.

## V Conclusions

In this work we have studied, by means of numerical methods, the structure and stability of domain-wall (DW) solutions in the system of two discrete nonlinear Schrödinger equations with the coupling of the cross phase modulation (XPM) type. The consideration of this problem is suggested by the analogy with known DW solutions in a standard model of a nonlinear optical fiber carrying two polarizations of light or two different wavelengths, as well as in the quasi-one-dimensional binary Bose-Einstein condensate (BEC). The results directly apply to an array of fibers of this type, with the anomalous intrinsic diffraction controlled by the direction of the light beam, or to a string of BEC drops trapped in an optical-lattice potential; in the latter case, the generic case with the positive scattering lengths is that which may give rise to DW patterns.

Using Newton-type methods and continuation from various initial patterns, starting from the anti-continuum (AC) limit, we have found a number of different stationary solutions of the DW type (while the continuum model admits a single type of the DW solution). Different stability scenarios and transitions to instability were identified for these solutions. In the case of strong XPM coupling, corresponding to two circular polarizations or two different wavelengths in the optical-waveguide array, natural DW configurations contain only one polarization at each end of the chain. The most fundamental solution of this type, generated by the simplest AC pattern, has been found to be always stable. Another solution, generated by the AC pattern that includes a vectorial state in the central site of the lattice, generates a behavior which is unusual for nonlinear dynamical lattices: it is unstable for small values of the coupling constant $`C`$ (which is the ratio of the nonlinear propagation length to the coupling length in the waveguide array), acquiring stability and remaining stable at larger values of $`C`$. Stable bound states formed by stable DWs were also found. A number of DW configurations generated by more sophisticated AC patterns were obtained, but they were either found to be completely unstable, or to be stable only at very small values of $`C`$.

In the case of weak XPM, which corresponds to linear polarizations in optics, the solutions become staggered and are subject to oscillatory instabilities. In this case, a more natural DW solution is the one with a combination of both polarizations, with the phase difference between them being $`0`$ and $`\pi `$ at the opposite ends of the lattice. This solution is unstable at all values of $`C`$, but the instability is very weak for large values of $`C`$, corresponding to stabilization in the continuum limit.

The robustness of all the patterns that were predicted to be stable in direct simulations was corroborated in direct simulations of the full nonlinear system. The evolution of unstable patterns was simulated too. In some cases, this instability is extremely weak; in other cases, a stronger instability leads to rapid rearrangement of the unstable DW into a stable one, including the possibility of generation of a moving DW in the weak-XPM model; a very strong instability may even induce “lattice turbulence”.

Estimates of physical parameters necessary for the formation of DWs in the optical-waveguide array and BEC string were given too. In particular, discrete optical DWs are expected to be found in the same region of parameters were bright discrete solitons have been already observed.

The consideration of DW patterns in dynamical lattices can be continued in several directions. A topic of direct interest concerns the mobility of DWs across the lattice. Note, in particular, that a mechanical twist applied to an optical fiber may give rise to a driving force acting on DWs in it , which can support the motion. On the other hand, a weak symmetry-breaking deformation of the waveguides will induce linear mixing between the two orthogonal polarizations, which will drastically affect the DWs. Another interesting problem is the interaction between DWs (in Ref. , it was shown that a bound state of two DWs with opposite polarities is possible in a nonlinear optical fiber in the presence of the twist). These issues will be considered elsewhere.


text:
# Observation of Long-lived Vortex Aggregates in Rapidly Rotating Bose-Einstein Condensates

## Abstract

We study the formation of large vortex aggregates in a rapidly rotating dilute-gas Bose-Einstein condensate. When we remove atoms from the rotating condensate with a tightly focused, resonant laser, the density can be locally suppressed, while fast circulation of a ring-shaped superflow around the area of suppressed density is maintained. Thus a giant vortex core comprising 7 to 60 phase singularities is formed. The giant core is only metastable, and it will refill with distinguishable single vortices after many rotation cycles. The surprisingly long lifetime of the core can be attributed to the influence of strong Coriolis forces in the condensate. In addition we have been able to follow the precession of off-center giant vortices for more than 20 cycles.

The phenomenon of quantized vortex formation is a unifying feature present in many quantum mechanical systems. Vorticity is intimately connected to superfluidity and thus is a premier means for fundamental and comparative studies of different “super” systems. In dilute gas Bose-Einstein condensates, vortex experiments have ranged from the study of individual or few vortices  and vortex rings  to the first creation of vortex lattices  and the study of systems containing large amounts of vorticity .

Using a refinement of our experimental technique, we have now been able to generate a giant vortex. We will use the term “giant vortex” to denote a region, containing multiple units of vorticity, in which the density is completely suppressed such that the individual vortices are no longer discernible. We note that some authors reserve this term for higher order phase singularities. The possible existence of stable giant vortices under certain conditions like supersonic flow in potentials that have strong quartic terms has been theoretically predicted (Ref. , see also Ref. ), and in Ref.  it is shown that a pinning potential can lead to stable multiquantum vortices. In our case, however, the giant vortex formation arises as a dynamic effect. Nevertheless the lifetime of our giant vortices can extend over many seconds, which we attribute to a stabilization of density features in a rapidly rotating condensate due to strong Coriolis forces. The influence of Coriolis forces can also induce oscillations of the giant vortex core size in the early stages of its evolution.

Our experimental starting conditions are similar to the ones described in our previous papers . We start with a cloud of uncondensed 87Rb atoms in the $`|F=1,m_F=1`$ state, held in a magnetic TOP trap and precooled to a temperature slightly above the critical temperature for Bose-Einstein condensation. By elliptically deforming the magnetic trap in the horizontal xy plane and suddenly changing the angle of the deformation, we impart angular momentum around the vertical z-axis to the cloud. We then restore axial symmetry and perform a quasi one-dimensional evaporation in a trap with frequencies $`\{\omega _\rho ,\omega _z\}=2\pi \{8.3,5.4\}`$Hz. This evaporation preferentially removes atoms close to the axis of rotation (z-axis), thus spinning up the cloud. At the end of this step, we routinely achieve fast rotating condensates consisting of $`3\times 10^6`$ atoms with a rotation rate of 95% of the radial trap frequency, with no detectable thermal cloud. This rotation rate is deduced from the aspect ratio, which changes from the value of 1.53 for a static cloud to 0.48 because of the centrifugal forces created by the rotation. These condensates typically contain 180 or more vortices as seen in Fig. 1(a) or Fig. 2(a). After careful optimization of the trap roundness and in the presence of a quasi 1-D rf-shield we are able to observe the BEC rotation for times exceeding 5 minutes. While we do lose atoms from the condensate over this time scale, the rotation rate remains at its initial value.

The evaporative spin-up mechanism is limited because eventually all the uncondensed atoms have been removed or condensed during the quasi one-dimensional evaporation. However, the angular momentum per particle in the condensate can be further enhanced by selectively removing atoms with low angular momenta by a resonant, focused laser beam sent through the condensate on the axis of rotation. Since the rotating BEC typically has a Thomas-Fermi radius of $`66\mu m`$, the width of the laser beam of about $`16\mu m`$, stated here as the full width at half maximum (FWHM) of the Gaussian intensity profile, is small enough to provide sufficient selectivity. The frequency of the laser is tuned to the $`F^{\prime \prime }=1F^{}=0`$ transition of the D2 line, and the recoil from spontaneously scattered photons blasts atoms out of the condensate.

Using this new technique we are able to substantially increase the vorticity in the BEC and can now routinely obtain condensates containing more than 250 vortices \[Fig. 1(b)\]. Since these vortex cores are too small to be imaged when the BEC is held in trap, the images in Figs. 1, 2, 3, and 6 are taken after having let the BEC expand. For this, a rapid adiabatic passage technique is employed to change the hyperfine state of the atoms from the original trapped $`|F=1,m_F=1`$ state to the antitrapped $`|F=2,m_F=1`$ state. Subsequently the minimum of the magnetic trapping field is rapidly moved below the position of the atoms so that the atoms are supported against gravity while at the same time they are radially expelled . After a chosen expansion time the magnetic fields are switched off and the atoms are imaged along the original axis of rotation. This technique allows us to expand the condensate to an adjustable diameter that can exceed 1.4 mm .

If the removal of central atoms is driven strongly enough, a giant vortex appears. The core is surrounded by a ringshaped superflow. The circulation of this superflow is given by summing up all the “missing” phase singularities in the core and can assume very large values. Using the fact that for large amounts of circulation the rotation of the velocity field can be obtained classically, we determine that in some cases the BEC supports a ring-shaped superflow with a circulation of up to 60 quanta around the core. We can easily control the amount of this circulation by changing the intensity or the duration of the laser beam that removes atoms along the axis of rotation.

In our experiment the formation of the giant vortex comes about in a sequence of very distinct stages as shown in Fig. 1. For this expansion image sequence a rapidly rotating BEC is first formed by our evaporative spin-up technique. Then the atom-removal laser is applied with a fixed power of 8 fW for a variable amount of time, followed by a 10 ms in-trap evolution time and a 45 ms expansion in the antitrapping configuration described above. During the atom removal an rf-shield is left on to constantly remove uncondensed atoms that tend to decelerate the condensate rotation. Figure 1(a) shows the result of only the evaporative spin-up. This particular condensate contains 180 vortices and has a Thomas-Fermi radius of $`63.5\mu m`$ when held in the trap. When the atom removal laser is applied for 14 s as in Fig. 1(b), the number of vortices is increased to 250 and the Thomas-Fermi radius to $`71\mu m`$. The rotation rate, determined from sideview aspect-ratio images, has increased from $`0.94\omega _\rho `$ to $`0.97\omega _\rho `$. After atom removal times of 15 to 20 s, the vortex lattice becomes disordered \[Fig. 1(c,d)\], and the giant vortex core starts to develop in the center \[Fig. 1(e,f)\]. An enlarged view of the core region in this early stage of giant vortex core formation is seen in Fig. 1(i), which nicely demonstrates how the individual vortices in the center consolidate. For the larger removal times shown in Figs. 1(g) and (h), a clear elliptical deformation is observed. This deformation is due to residual static asymmetries of the trapping potential that can nearly resonantly drive quadrupolar surface modes when the rotation frequency of the condensate approaches the radial trapping frequency. We have verified this explanation by deliberately changing the roundness of the trapping potential and observing corresponding changes in the ellipticity and in the direction of the ellipse. If the trap roundness is not optimized, the ring-shaped superflow has a tendency to fragment into typically 5 to 8 blobs that continue to rotate around the core.

If the atom removal laser is applied only for a limited amount of time, the giant vortex will fill in and the vortices refreeze  into a well ordered lattice again after a sufficiently long evolution time. An example is shown in Fig. 2. Here, the starting point is the condensate in Fig. 2(a), containing 190 vortices and rotating at $`0.95\omega _\rho `$. A pulse from the atom removal laser creates the core seen in Fig. 2(b). Crudely assuming that all atoms, and only those atoms, originally within a cylindrical volume of one-third of the original cloud radius are removed, we can perform an integral over the Thomas-Fermi density profile with a rigid-body-rotation velocity distribution, and find that 23% of the atoms have been removed but only 4% of the angular momentum. Since in equilibrium the total number of vortices should be linear in the average angular momentum per atom , we would expect that the original 190-vortex cloud should reequilibrate to form a cloud with 236 vortices, in rough agreement with the value of 260 vortices observed in the reequilibrated cloud shown in Fig. 2(d).

In our case the formation of the giant vortex is not connected to a repulsive conservative optical potential produced by the central laser beam, but only due to the removal of atoms from the axis of rotation by spontaneous photon scattering. To prove this, we have varied the laser frequency over the range from -6 MHz to +6 MHz around the $`F^{\prime \prime }=1F^{}=0`$ resonance and in all cases have been able to generate a giant vortex.

Since theoretical studies predict the appearance of a giant vortex if the trapping potential has a sufficiently strong admixture of quartic terms, it is important to quantify the strength of the quartic terms in our case. For this we have calculated the term of the TOP trap potential that depends on the radial position as $`r^4`$. At the outer edge of the condensate this quartic term produces a correction to the potential of less than $`10^3`$, compared to the effective $`r^2`$ term, even though the latter is much weakened by the centrifugal force. Higher order corrections contribute even less. Empirically, the refilling and refreezing we observe (Fig. 2) indicates that the giant vortex is not an equilibrium configuration but has instead only dynamical stability.

We suggest an intuitive, classical picture describing the formation as a dynamical effect: the removal of atoms from the center of the condensate produces a pressure gradient due to mean field energy that tries to drive atoms from the outer regions into the center, so as to close the hole. Due to Coriolis forces, however, atoms moving radially towards the center are deflected and assume a fast azimuthal motion around the core rather than filling the core, thus creating the giant vortex.

Some exotic features of the early stages of core formation are revealed when we apply only a very short, weak atom removal pulse and observe the subsequent evolution as in Fig. 3. Here, the atom removal laser has a power of 2.5 pW and a fixed pulse length of 5 ms. Its FWHM of $`16\mu m`$ is approximately twice the lattice spacing, $`7\mu m`$. By varying the delay time between the laser pulse and the expansion, we see (Fig. 3) that the core formation clearly lags behind the laser pulse. Fig. 3(e) shows a zoomed-in view of the core region of Figure 3(c). It is very interesting to observe that for these short, weak pulses the density depression appears to develop in discrete steps, and the step boundaries follow the hexagonal lattice pattern.

When applying stronger laser pulses, we detect clear, damped oscillations of the core area as shown in Fig. 4. We study these oscillations by analyzing in-trap images taken after a variable evolution time between the end of the atom-removal pulse and the expansion, such as the ones shown in the inset of Fig. 4. The oscillation frequency depends on the initial conditions; for the case of Fig. 4 with initially $`2.2\times 10^6`$ atoms at a rotation rate of $`0.9\omega _\rho `$ we obtain a frequency of $`3.5\omega _\rho `$. Decreasing the initial rotation rate of the condensate leads to faster core-oscillation frequencies and increases the amplitude of this oscillation; for rotation rates below $`0.9\omega _\rho `$ in-trap images even show a near complete closure of the core, followed by the core opening again (see inset of Fig. 4). Presumably these oscillations arise when a sudden removal of atoms leaves forces from the density gradient and Coriolis forces initially out of equilibrium. The core oscillation is not observed to be related to an overall breathing mode of the condensate.

In addition we can observe the precession  of an off-center giant vortex. For this, the atom removal laser is deliberately offset from the center of cloud rotation. The duration of the removal pulse is kept short (10 ms) in comparison to the initial vortex lattice rotation period ($`2\pi /(0.95\omega _\rho )=126ms`$) so as to create only a local hole. By using different laser powers, we can vary the size of the core from a small hole as in Fig. 5(a) all the way to the extreme case of a big hole that only leaves a crescent segment of the BEC as shown in Fig. 5(c), and the precession of these holes can be monitored by applying variable evolution times of the trapped BEC after the atom removal pulse is finished. In all three cases the measured precession frequency is approximately $`\omega _\rho `$. In the case of Fig. 5(b) and (c) we are able to follow the precession for more than 20 cycles. If, instead, the laser is left on for approximately a full lattice rotation cycle, a complete ring can be cut out of the condensate \[Fig. 6(a)\], which eventually breaks up into many individual blobs \[Fig. 6(c)\], each of which supports remnants of the original lattice. Both the long-lived core precession as well as the ring structure in Fig. 6 are impressive demonstrations of the stability of density features due to Coriolis forces in rapidly rotating condensates.

Finally, we also observe that the sudden extraction of atoms can launch transverse modes in the lattice , which will be the subject of a future publication.

The work presented in this paper was funded by NSF and NIST. PE acknowledges support by the Alexander von Humboldt Foundation.