text:
# Molecular Gas and Star Formation in Bars of Nearby Spiral Galaxies

## 1 Introduction

Barred spiral galaxies are ideal laboratories for the study of star formation because they host a variety of environments with distinctive star formation activity and gas dynamics. These environments include the circumnuclear, inner and outer rings (see discussion of rings in Buta & Combes 1996; Regan et al. 2002), the bar ends, and the bar region itself, located in between the bar ends and the nucleus. From region to region, the star formation activity can vary dramatically: bars have star formation rates of $``$0.1–0.4 $`M_{\mathrm{}}`$yr-1 (e.g., Martin & Friedli 1997), whereas in circumnuclear rings, in an area 10-30 times smaller, star formation rates may be as high as 1 $`M_{\mathrm{}}`$yr-1 (e.g., Buta et al. 2000). Since star formation occurs in molecular clouds, comparative studies of the distribution of molecular gas and star formation activity in different environments can shed light on how star formation may be induced or inhibited.

Along the spiral arms, for example, observations find that the molecular gas and HII regions are usually not co-spatial; HII regions are preferentially offset towards the leading side of the molecular gas or dust lanes (Vogel, Kulkarni & Scoville, 1988; Rand, 1993; Knapen & Beckman, 1996; Loinard et al., 1996). This offset is interpreted as evidence of star formation induced by a spiral density wave (Vogel, Kulkarni & Scoville, 1988; Rand, 1993). However such offsets are not universally present. In M100, Sempere & Garcia-Burillo (1997) find that the offset is absent, or even inverted, along the spiral arms. Still, in other cases, the offset is more pronounced between the HII regions and the dust lanes as the molecular gas and dust lanes diverge (e.g., Lord & Kenney 1991; Rand, Lord & Higdon 1999). Rand, Lord & Higdon (1999) attribute the divergence to heating of the gas by young stars or cosmic rays, or to a two-phased molecular gas medium, but admit that neither of these explanations works satisfactorily.

In barred spirals, most previous studies comparing the molecular gas distribution and star formation activity have focused on the highly active circumnuclear region (e.g., Eckart et al. 1991; Roy & Belley 1993; Kenney, Carlstrom, & Young 1993; Sakamoto et al. 1995; Benedict, Smith, & Kenney 1996). In this paper, we focus on the unique region between the bar ends and the nucleus. This region is unlike any other in the galactic disk because it is dominated by highly elliptical stellar orbits (see reviews by Sellwood & Wilkinson 1993; Athanassoula 1992a, b). Throughout this paper, we refer to it simply as the bar.

Only one previous study has studied the location of HII regions relative to the stellar bar (Martin & Friedli, 1997); they found that the “H$`\alpha `$ bar” was usually offset towards the leading side of the stellar bar with misalignments as large as 15. Studies of molecular gas or dust in bars have found that the gas and dust are also on the leading side of the stellar bar (e.g., Ondrechen 1985; Handa et al. 1990; Regan & Vogel 1995; Downes, Reynaud, Solomon, & Radford 1996; Sheth et al. 2000). But the relationship between the molecular gas (or dust) and star formation has only been studied in a few cases. In M101, for example, Kenney, Scoville, & Wilson (1991) note that the molecular gas and H$`\alpha `$ emission are at the same position angle. In contrast, Sheth et al. (2000) found that the HII regions in the bar of NGC 5383 are offset towards the leading side of the bar dust lanes. It is unclear whether such offsets are common in bars. Since bars have distinctive (and well-studied) gas kinematics, a study of the relative distribution of the gas and stars can further elucidate the complex phenomenon of star formation. With this goal, we have studied six barred spirals from the recently completed BIMA Survey of Nearby Galaxies (SONG) (Regan et al., 2001; Helfer et al., 2002b).

The sample selection is discussed in §1.1, and the observations and data reduction in §2. Using the CO (J=1-0) emission line to trace the molecular gas, and the ionized hydrogen (H$`\alpha `$ ) line to trace recent star formation activity, we compare the relative distribution of the two in §3.1. In all six galaxies we find that the majority of the H$`\alpha `$ emission is offset towards the leading side of the molecular gas.

We quantify the offset using one and two dimensional cross-correlation analysis in §3.2. The CO-H$`\alpha `$ offset constrains how stars may form depending on the gas flow into the dust lane. We discuss the results of the cross-correlation analysis in the framework of the two main classes of bar gas flow models (the N-body sticky particle and hydrodynamic models) in §4, and summarize our results in §5.

### 1.1 Sample Selection

In BIMA SONG (see details of the survey in Regan et al. 2001; Helfer et al. 2002b), we detect CO emission in 27 of the 29 barred spirals (e.g., Figures A1 and A2 in Sheth 2001). Typically, the emission is detected in the circumnuclear region, where it is usually the brightest. In some bars, CO emission is also detected at the bar ends, along the bar, and even in inner rings (Regan et al., 2002). Since we are mainly interested in studying star formation in the bar region we limit ourselves to those BIMA SONG galaxies in which both CO and H$`\alpha `$ emission are clearly detected over a significant portion of the bar. Six galaxies satisfy this criterion: NGC 2903, NGC 3627, NGC 4321, NGC 5457, NGC 6946 and IC 342. Global properties of these six are listed in Table 1.

Though small, our sample spans a range of Hubble types with 1 Sb, 2 Sbc and 3 Scd galaxies. This range may be important because most of the observed differences in star formation activity in bars occur between early and late Hubble type galaxies (e.g., Elmegreen & Elmegreen 1985; Ohta et al. 1986). For instance, early Hubble type bars have low star formation activity along the bar and high star formation activity at the bar ends, whereas late Hubble type galaxies have higher star formation activity in the bar, but have a gap in star formation at the bar ends (Phillips, 1996).

All the galaxies in our sample are classified as SAB in the RC3 catalog. However this does not mean that they are all weak bars because the Hubble classification of SAB types is not rigorous; in fact, a recent analysis of infrared data has shown that the true fraction of strong bars may be as high as 56%, higher than the typically quoted 33% in the RC3 (Eskridge et al., 2000). The strength of a bar may be correlated with its Hubble type. Elmegreen & Elmegreen (1985) concluded that early Hubble type galaxies have stronger bars because these bars are longer, relative to their disks, and have flat profiles; these galaxies also have strong spiral arm patterns. However, bar strength is a difficult parameter to quantify; other structural properties such as the bar ellipticity and bulge size are also important (see discussion in Buta & Block 2000).

A good indicator of the bar strength is the shape of the bar dust lanes because it reflects the gas response to the stellar bar. Athanassoula (1992b) showed that strong bars have relatively straight dust lanes whereas weaker bars have curved dust lanes. Extending her analysis to our sample, we infer that NGC 2903 and NGC 3627 are strongly barred because of their straight dust lanes. NGC 5457 and NGC 4321 with slightly curved dust lanes are weaker bars. IC 342 has an even more curved dust lanes indicating an even weaker bar. In NGC 6946, the northern dust lane appears to be straight but the southern dust lane is curved so it is difficult to classify this bar solely on the shape of the dust lanes. Regan & Vogel (1995) suggest that if there is bar in NGC 6946, it is rather weak. Using bar dust lanes as a measure of bar strength, we find that the three latest Hubble types in our sample are classified as relatively weak bars. In summary, our sample contains a broad range of bar strengths, from strong bars in NGC 2903 and NGC 3627, to intermediate strength bars in NGC 4321 and NGC 5457, to relatively weak bars in NGC 6946 and IC 342.

## 2 Observations and Data Reduction

### 2.1 Molecular Gas Data

All six galaxies were observed in the CO (J=1–0) emission line with the BIMA (Berkeley-Illinois-Maryland Association) array and the NRAO111The National Radio Astronomy Observatory is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc. telescope at Kitt Peak 12m single dish telescope as part of the BIMA SONG key project. The 44 galaxy SONG sample was chosen with the following criteria: heliocentric velocity, VHEL $`<`$ 2000 km s-1, declination, $`\delta >`$ -20o, inclination, $`i<`$ 70o, and apparent magnitude, B\_T$`<`$ 11. The typical data cube has a synthesized beam of 6$`\mathrm{}`$, a field of view of 3$`\mathrm{}`$. In a 10 km s-1 channel, the typical noise is $``$58 mJy beam-1. Further details of the data acquisition and reduction can be found in Regan et al. (2001) and Helfer et al. (2002b).

### 2.2 H$`\alpha `$ Data

We observed NGC 2903, NGC 3627 and NGC 4321 at the 0.9m telescope at Kitt Peak222Kitt Peak National Observatory, National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc. (AURA) under cooperative agreement with the National Science Foundation on the nights of 4–6 April 1999, with the T2KA 2048$`\times `$2048 CCD camera in f/13.5 direct imaging mode. In this mode the camera has a 13$`\stackrel{}{\mathrm{.}}`$1 field of view with 0$`\stackrel{}{\mathrm{.}}`$384 pixels. In the broad band R-filter, we took three exposures of 180s each, and in the H$`\alpha `$ filter ($`\lambda _0`$ = 6573 Å, $`\mathrm{\Delta }\lambda `$ = 67 Å) we took three exposures of 100–140s each. We divided each frame by a normalized flat-field frame and combined the resulting frames into one single image for each corresponding filter. Then from these images we subtracted the mean sky brightness, and removed cosmic rays using standard routines in the NOAO/IRAF333IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation software package. We also corrected all images for atmospheric extinction. Finally, we registered foreground stars in each image with the Hubble Guide Star Catalog (GSC) and determined astrometric solutions for each image. The residuals in determining absolute positions were smaller than 0$`\stackrel{}{\mathrm{.}}`$2; however, systematic errors in the GSC prevent us from achieving an accuracy $`<`$ 1$`\mathrm{}`$. We used the R-band image to subtract the underlying continuum from the H$`\alpha `$ image.

We observed IC 342 at the 1.5m telescope at Palomar444Observations were made on the 60 inch telescope at Palomar Mountain, which is jointly operated by the California Institute of Technology and the Carnegie Institution of Washington. on the night of 19 November 1999 with the Maryland-Caltech Fabry-Perot camera. These data were reduced and calibrated using the procedure described in Vogel et al. (1995). The H$`\alpha `$ image for NGC 6946 is also a similar velocity-integrated H$`\alpha `$ Fabry-Perot map obtained previously at Palomar. A broadband continuum image and a continuum-subtracted H$`\alpha `$ image for NGC 5457 were supplied to us by R. Kennicutt; these images have 2$`\stackrel{}{\mathrm{.}}`$6 pixels compared to 0$`\stackrel{}{\mathrm{.}}`$38 pixels for the other optical images.

## 3 Results

### 3.1 Comparing the Distribution of Molecular Gas and Star Formation

We first consider the relative distribution of the molecular gas and star formation activity in the six galaxies by overlaying the CO emission on the continuum-subtracted H$`\alpha `$ images (Figures 16). Our goal is to learn where the stars may be forming by comparing the relative locations of molecular clouds and HII regions. Since any image is only a snapshot in time it is difficult to assign a star forming region to its parent molecular complex. We make the usual assumption that CO traces regions most prone to star formation, and H$`\alpha `$ traces young stars most likely to have formed in molecular gas complexes traced by the CO.

In every image we show the extent of the stellar bar, as determined by Sheth (2001), with an ellipse. The ellipse corresponds to the outermost bar isophote and was determined from fitting ellipses to optical or near-infrared images following a method devised by Regan & Elmegreen (1997). The width of a bar is a difficult parameter to measure accurately because contamination by bulge light tends to make the fitted ellipses wider (Sheth, 2001). Since the minor axis length is not relevant for the study presented here, we chose a representative minor axis that is one-quarter that of the major axis for all bars except IC 342, which has a small, oval bar for which we chose a minor axis that is one-half that of the major axis.

In each image (Figures 16) we show CO contours overlaid on a grey-scale image of the continuum-subtracted H$`\alpha `$ emission. We describe the distribution of H$`\alpha `$ and CO using approximate position angles because their morphology is often curved. Compact H$`\alpha `$ structures are labeled with a prefix H, and CO features are labeled with a prefix C. We also show the direction of rotation of the bar with a curved arrow symbol in the lower right hand corner; the direction is based on the assumption that spiral arms trail in galaxies. All the galaxies in our sample rotate counter-clockwise except for NGC 6946.

In all six galaxies we find that the CO emission is brightest on the leading side of the stellar bar. Overlays of the CO maps on optical-infrared color maps (not shown) confirm that the molecular gas generally coincides with the bar dust lanes (Sheth et al., 1998). In five out of the six galaxies, weak spurs of CO emission are seen on the trailing (upstream) side of the dust lanes. At the bar ends, bright CO emission is present on the leading and trailing sides; this gas is most likely associated with inner rings, which often encircle bars (Buta & Combes, 1996; Regan et al., 2002). The H$`\alpha `$ emission is distributed in bright, compact structures and weak, diffuse emission. The differences in bar star formation activity noted by Phillips (1996) are not always present. For instance, both the early Hubble type bar in NGC 3627 and the late Hubble type bar in NGC 5457 have relatively little star formation activity, whereas strong activity can be seen in the bar of NGC 2903. But such differences across Hubble type are difficult to quantify with this small sample. In these six bars, while there are a few cases where H$`\alpha `$ is coincident with, or on the trailing side of CO, all six show that majority of the H$`\alpha `$ emission is offset towards the leading side of the CO emission. We describe these results for each individual galaxy in detail below:

NGC 2903 (Figure 1): The stellar bar in NGC 2903 is oriented at a position angle (PA) of 24 east of north. The molecular gas is distributed towards the leading side of the stellar bar, with a PA of $``$30. The H$`\alpha `$ emission is offset even farther towards the leading side of the CO emission.

In the northern half of the bar, the CO emission is brightest in a large complex, labeled C1. A bright H$`\alpha `$ complex, H1, coincides with part of the C1 complex. 3–4$`\mathrm{}`$ south of C1 is a narrow ridge of CO, labeled C2, which bridges C1 with the bright circumnuclear gas. Several discrete H$`\alpha `$ features (H2, H3, H4, and H5) are seen on the leading side of the CO. These structures are connected by a diffuse, concave ridge of H$`\alpha `$ emission. On the trailing (upstream) side, some diffuse H$`\alpha `$ emission is present but there is a conspicuous lack of HII regions; this is not an artifact of dust extinction. In the dust lanes, the dust extinction is usually insufficient to hide HII regions (e.g., Regan & Vogel 1995). Also, the gas density is expected to be dramatically lower on the trailing side of the dust lane as the gas piles up in the bar dust lane (e.g., see model density distributions in Combes & Gerin 1985; Athanassoula 1992b; Piner, Stone & Teuben 1995).

In the southern half of the bar, there are two bright CO complexes, C3 and C4. As in the northern half, H$`\alpha `$ complexes (H6, H7, and H8) are all on the leading side of the CO emission. A diffuse ridge of H$`\alpha `$ emission extends south of H7, with a position angle of $``$5, on the leading edge of the CO emission. In this half of the bar, we also see a weak CO spur, labeled C5, east of C4 on the trailing side of the dust lane. Diffuse H$`\alpha `$ emission is present at the trailing edge of C5. In two places, near H7, and north of H6, there is weak CO emission on the leading side of the dust lane.

NGC 3627 (Figure 2): The bar in NGC 3627 is oriented at a PA of 161. The molecular gas is mostly on the leading side the stellar bar, with a PA of $``$164 , but note that at the southern end of the bar, the molecular gas morphology turns sharply towards the trailing side of the bar. The H$`\alpha `$ emission is offset towards the leading side of the CO. The CO intensity decreases from the inner ends of the bar towards the outer ends, i.e. from C3$``$C2$``$C1 in the northern half, and from C4$``$C5$``$C6 in the southern half. This behavior is consistent with the predictions of the hydrodynamic models (Sheth et al., 2000).

In the northern half of the bar, there is only one bright, compact H$`\alpha `$ structure, labeled H1, located 30$`\mathrm{}`$ north of the nucleus, and 3–5$`\mathrm{}`$ east, and towards the leading side of the CO ridge connecting C1 and C2. Diffuse H$`\alpha `$ emission connecting H1 to the bar end is also present on the leading edge of the CO. As in NGC 2903, there is a paucity of H$`\alpha `$ emission on the trailing side of the CO. Additional CO and H$`\alpha `$ emission associated with an inner ring (shown by the dashed segments) circumscribing the bar, is also present; this inner ring is discussed in detail in (Regan et al., 2002).

In the southern half of the bar, all three compact H$`\alpha `$ structures (H2, H3, and H4) are also on the leading edge of the CO emission, 3–4$`\mathrm{}`$ west of the CO ridge connecting C4, C5, and C6. As in the northern half, diffuse H$`\alpha `$ emission connects H2, H3, and H4; most of this emission is also on the leading edge of the CO.

NGC 4321 (Figure 3): The bar in NGC 4321 is oriented at a PA of 102. The molecular gas is mostly on the leading side of the stellar bar with a PA$``$90 ; this is more clearly seen in the western half of the bar where the CO emission runs continuously down the entire half of the bar. In the western half, the H$`\alpha `$ emission is offset another 3-4$`\mathrm{}`$ towards the leading side, with approximately the same PA. In the eastern half, the H$`\alpha `$ emission is weak and no clear offset is seen.

In the western half of the bar, a relatively narrow lane (10$`\mathrm{}`$ in width) of CO emission emerges from the bright, central 40$`\mathrm{}\times `$ 20$`\mathrm{}`$ oval-shaped distribution of CO emission. This lane extends to the west (C1$``$C2$``$C3) and then curves inwards to the spiral arms from C3 to C4. As in NGC 3627, the CO intensity decreases away from the central region. A spur-like feature, labeled C5, is seen 5$`\mathrm{}`$ south of C1 extending almost 10$`\mathrm{}`$ ($``$780 pc) on the trailing side of the bar. These spurs may play an important role in forming stars because they have high gas density but are located in an area with lower shear (Sheth et al., 2000).

The discrete H$`\alpha `$ structures (H1, H2, and H3) are evenly-spaced ($``$10$`\mathrm{}`$), and are on the leading edge of the CO emission. A narrow ridge of diffuse H$`\alpha `$ emission connecting these structures is offset north from the CO by about 4$`\mathrm{}`$. Farther west, towards the bar end, a bright concentration of molecular gas, labeled C4, is coincident with the brightest H$`\alpha `$ emission in this half of the bar, labeled H4.

In the eastern half of the bar, the two main concentrations of CO, C6 and C7, are on the leading side of the stellar bar. Two HII regions, H5 and H6, are offset $``$5$`\mathrm{}`$ west of C7 and C6 respectively. Unlike the western half, there is no clear leading offset of HII regions from the CO emission. At the eastern bar end, bright CO and H$`\alpha `$ emission are present in the regions labeled C8 and H7.

CO and H$`\alpha `$ emission are distributed on the trailing and leading sides at the bar ends (10$`\mathrm{}`$ north of C4, and 12-15$`\mathrm{}`$ south of C8 and H7). Such a morphology seems to be common in bars and is most likely associated with inner rings (Buta & Combes, 1996) that encircle bars.

NGC 5457 (Figure 4): The bar in NGC 5457 is oriented at a PA of 80. The molecular gas is primarily on the leading side of the stellar bar with a PA$``$95 . The H$`\alpha `$ emission tends to be towards the leading side of the CO emission.

In both halves of the bar the CO intensity decreases outwards from the central concentration of gas, to the east towards C1, and to the west towards C6. In the western half of the bar, there is one bright CO complex, labeled C2 in Figure 4. Generally coincident with but offset slightly towards the leading side of C2 is a bright HII region, labeled H1. Diffuse H$`\alpha `$ emission connecting the nuclear region and H1 is seen mostly on the leading side of the CO emission in the region in between C1 and C2. To the south of C2, on the trailing side, there is a 10$`\mathrm{}`$ spur-like CO feature, labeled C3.

In the eastern half, the CO morphology is very similar to that seen in the western half. C7 is a weak CO peak, $``$30$`\mathrm{}`$ east of the circumnuclear region. Associated with C7 is weak H$`\alpha `$ emission, labeled H3. Generally the H$`\alpha `$ emission in this half of the bar is characterized by a continuous ridge extending from the circumnuclear emission to H4, on the leading side of the CO. H4, like H2, is on the leading side of the bar and is probably the beginning of an inner ring.

Towards both bar ends, the CO emission deviates away from the leading side of the bar and connects with the CO emission in the spiral arms of NGC 5457. The CO region labeled C4, and the H$`\alpha `$ region labeled H2, are probably associated with the bar end, but the two bright CO and H$`\alpha `$ peaks, 6–10 arcseconds north of H2, are associated with the spiral arm.

NGC 6946 (Figure 5): The stellar bar in NGC 6946 is oriented at a PA of 19. In the bar, the molecular gas is distributed asymmetrically about the center. To the north of the circumnuclear region, the CO emission covers a rather broad 40$`\mathrm{}\times `$40$`\mathrm{}`$ area, with a sharp ridge of CO emission on the leading side of the bar with a PA of $``$0. In the southern half of the bar, the broadly distributed CO emission is truncated to the south, extending only 20$`\mathrm{}`$ south of the circumnuclear area. While the CO and H$`\alpha `$ emission overlap considerably, the brightest H$`\alpha `$ emission tends to be on the leading edge of the CO.

In the northern half of the bar, a ridge of CO emission extending from the circumnuclear region towards C1 gradually decreases in intensity. Along this ridge there is an HII region, 10$`\mathrm{}`$ north of the nucleus labeled H1, a CO peak, 10$`\mathrm{}`$ farther north labeled C1, and two H$`\alpha `$ peaks, labeled H2 and H3 even farther north of C1. We find a significant amount of CO emission on the trailing side of this bar. There is a concentration of CO emission, labeled C2, curving east from C1. An HII region is coincident with it. On the trailing side of the bar, there is CO emission C3, C4, and C5, and H$`\alpha `$ peaks at H4 and H5.

In the southern half a ridge of CO emission extends from the nuclear region towards C7. Associated with C7 is the brightest H$`\alpha `$ emission in this half of the bar. Farther south of C7, the CO emission narrows into a ridge, labeled C9, and then ends in a broader peak at C10. On the trailing side of C9 is an HII region labeled H8. Two HII regions, H9 and H10, are seen in the leading side of the bar, a few arcseconds northeast of C10. There is also an trailing extension of CO emission to the east of H7, labeled C8.

IC 342 (See Figure 6): IC 342 has a small, oval stellar bar. The bar ellipse shown in Figure 6 is at a PA of 28. The CO emission appears to be in the center of the bar with a PA of $``$0–10. This central distribution of gas is unusual because it is expected only when the bar is extremely strong and inner Lindblad resonances are absent (Athanassoula, 1992b). Such a situation is extremely short-lived in the models. It is more likely that the bar parameters for this galaxy are misidentified because the bar resides in a bright starburst circumnuclear region. A recent study by (Crosthwaite et al., 2001) suggests that although the bar is most prominent in the central 2$`\mathrm{}`$, CO kinematics indicate a gaseous bar as long as 4$`\stackrel{}{\mathrm{.}}`$7 with a PA$``$0. In that case, the CO emission would be on the leading edge of the bar.

As in other bars, the CO emission along both dust lanes decreases in intensity outwards from the circumnuclear region towards C2 to the north, and C4 to the south. In the northern half of the bar there is an extension of CO, labeled C1, towards the trailing side. A few arcseconds north of it is an HII region, labeled H1, which also lies along the leading edge of the CO emission. Farther north near the bar end, there is a CO emission peak, labeled C3. A few arcseconds north and east of C3 is another H$`\alpha `$ peak, labeled H2.

In the southern half of the bar there is a concentration of CO emission, labeled C5, on the trailing side of the main lane of CO emission. There is also a ridge of CO emission which extends southeast of the nucleus toward the bar end, ending in a CO peak, labeled C6. In the bar there is a lot of diffuse H$`\alpha `$ emission but the few HII regions (H3 and H4) are on the leading side of the CO ridge.

In summary, CO emission in all six bars is brightest on the leading side of the stellar bar in the bar dust lanes. The CO intensity generally decreases towards the bar ends. Weak spurs of CO emission are sometimes seen on the trailing side of the bar. At the bar ends, there is often bright CO and H$`\alpha `$ emission on both the trailing and leading sides; this emission is probably associated with inner rings. In a few instances H$`\alpha `$ emission is coincident with, or on the trailing side of the CO emission. But in all six bars, the majority of the H$`\alpha `$ emission is on the leading side of the molecular gas.

### 3.2 Quantifying the Offsets

The typical distance between molecular gas complexes and HII regions sheds light on when, where and how stars form. To accurately measure this distance, we need to trace the orbit of an HII region and measure the separation along this path. In a typical galaxy disk, this path is usually circular and therefore a measurement of the azimuthal displacement is sufficient. However non-axisymmetric perturbations such as bars and spiral arms can introduce streaming motions in the gas (e.g., Roberts, Huntley & van Albada 1979; Athanassoula 1992b and references therein). In these cases, an orbit has both an azimuthal and a radial component. The radial displacement depends on the strength of the bar or spiral arm, and it may vary with radius. We expect that the radial streaming motions will be most important at the inner ends of the bar and spiral arms (Roberts, Huntley & van Albada, 1979). From §3.1, we infer that, to first order, the displacement is mostly in the azimuthal direction. Therefore we first quantify the distance between the CO and H$`\alpha `$ by doing a one-dimensional cross-correlation in azimuth (as described in the next paragraph), and then we perform a two-dimensional (radial + azimuthal) cross-correlation over the entire bar region.

#### 3.2.1 Cross-Correlation Analysis

We re-gridded each H$`\alpha `$ image to the CO image using standard MIRIAD routines (Sault, Teuben & Wright, 1995), deprojected both images using ZODIAC (Miyashiro, 1982; Shopbell, 1997) routines developed by ourselves and Gruendl (1996), and then mapped the images on to a polar grid. The deprojected CO (top left) and H$`\alpha `$ (bottom left) images are shown in the left column of Figures 712. For the one-dimensional cross-correlation analysis, we cross-correlated 1$`\mathrm{}`$ annuli. The result of these analyses are shown in the vertical panel on the right hand side of Figures 712. The star symbols indicate the angular offset corresponding to the maximum of the cross-correlation for that radius. Note that this value averages the offset between the CO and H$`\alpha `$ in the two halves of the bar. The solid circles on the deprojected H$`\alpha `$ and CO images are the boundaries within which the cross-correlation is done. The inner radius is chosen to best avoid the circumnuclear emission and the outer radius is at the end of the bar (corresponding to the bar ellipse in Figure 16). The dashed circles are typically spaced at integral multiples of 10$`\mathrm{}`$ and the corresponding radius is shown with a solid, horizontal line in the one dimensional cross-correlation diagrams in the right hand panel.

The cross-correlation result can be affected by systematic dust extinction. For example, consider a barred spiral where the H$`\alpha `$ emission is distributed evenly on the trailing and leading sides. The cross-correlation peak of such a distribution should not show any azimuthal offset. However, if the H$`\alpha `$ on the trailing side of the dust lane is obscured by a screen of dust, then the cross-correlation diagram would show that the H$`\alpha `$ is preferentially on the leading side of the dust lane. But as noted earlier in 3.1 this is unlikely because the dust extinction on the trailing side is very low as gas piles up in the dust lanes. Note that any extinction suffered by the H$`\alpha `$ in the dust lane only increases the observed offsets; so the CO-H$`\alpha `$ offsets calculated here are upper limits.

It is important to understand whether the cross-correlation peak depends on small, bright structures, like the H$`\alpha `$ and CO peaks, or on large, diffuse, and dim features. We performed two tests to address this question. First, we flagged all low level emission in the CO and H$`\alpha `$ maps of NGC 2903 and cross-correlated the images in radius and azimuth. We found that the cross-correlation peak was essentially unchanged. The main difference was that a background of small but non-zero correlation values, present in the original cross-correlation diagram, was now removed. The diffuse emission, therefore, is analogous to a constant offset contributing to the correlation diagram, but at no particular lag. This test indicates that the cross-correlation peak is determined primarily by the bright peaks in the maps.

As a second test, we set all emission above an arbitrary bright signal value to that signal value in the H$`\alpha `$ and CO maps for NGC 2903, approximating the removal of bright peaks. In the resulting correlation diagram there was a trend of H$`\alpha `$ offset azimuthally towards the leading side from the CO but the cross-correlation peak moved significantly from its original location and became broader. The larger width of the peak was expected because we removed the sharpness inherent in the maps being correlated. The shift, once again, showed that the location of the correlation peak is primarily determined by bright H$`\alpha `$ and CO peaks. Assuming that these peaks trace regions most prone to star formation, and recently formed OB associations respectively, the correlation peak effectively quantifies the offset between parent molecular gas complexes and newly formed stars.

#### 3.2.2 Results of the One-Dimensional Cross-Correlation

In every bar (right hand panels in Figures 712), we find that the peak of the cross-correlation shows H$`\alpha `$ azimuthally offset towards the leading side of the CO emission, consistent with the discussion in §3.1. In these panels, we have annotated the correlation peaks with the notation used in Figures 16 for the corresponding CO and HII regions.

The offset between the annotated regions varies from 0–800 pc with a range of offsets in any given bar. We find that the largest offsets are generally seen in the strongest bars. For instance, in NGC 2903, offsets of $``$800 pc are seen in the northern half between C2 & H5, and between H4, H5 & the northern part of C1. In the same bar, smaller offsets are also present between C3 & H6 ($``$400 pc), and between C1 & H1 ($``$200 pc). In the slightly weaker bars, NGC 4321 and NGC 5457, the HII regions are typically offset 200-500 pc (e.g., C1 & H1 in Figure 9), and in the weakest bars, NGC 6946 and IC342, the offsets are even smaller, typically less than 400 pc.

At some radii, the cross-correlation peak is at larger distances, exceeding a kiloparsec in some cases. These are instances where only weak, and diffuse H$`\alpha `$ emission is being correlated with the CO emission. An example of this is seen between 20–25$`\mathrm{}`$ in NGC 3627 (Figure 8). In a few cases we get non-physical values for the cross-correlation peaks; this occurs when we correlate either noise features or weak emission that is well outside the bar. The arrows point towards the non-physical value of the cross-correlation peak, off the scale shown in these figures.

Though there appears to be a relationship between the bar strength and the CO-H$`\alpha `$ offset, there is no systematic pattern in the offset as a function of the galacto-centric radius. We often see a trend of increasing offsets of H$`\alpha `$ from CO over radial extents of 5–10$`\mathrm{}`$, e.g., the radial bins between 38–42$`\mathrm{}`$ (C2 & H5) in Figure 7. This trend is due to beam smearing and subsequent deprojection of the images. These results are discussed further in §4.

#### 3.2.3 Results of the Two-Dimensional Cross-Correlation

The two-dimensional correlation results are plotted in the bottom panel in Figures 712. In these panels, the intersecting straight lines indicate the 0,0 lag position. Positive (negative) correlation values indicate H$`\alpha `$ offset radially inwards (outwards) or azimuthally leading (trailing) the CO. Unlike in the one-dimensional analyses, where we only correlated 1$`\mathrm{}`$ annuli, the location of the cross-correlation maximum (or maxima) in the two-dimensional analyses is susceptible to bright emission near the circumnuclear or bar-end regions because we sum over the whole bar region. While we did our best to choose an inner boundary to avoid most of the circumnuclear emission, the outer boundary sometimes included bar-end/spiral arm emission (for instance, in NGC 2903). This is one of the reasons why the two dimensional cross-correlation results are less reliable than the other analyses presented here. Another is the possibility of getting incorrect results because we may correlate physically unrelated peaks. Suppose the brightest HII region is at the outer bar end whereas the brightest CO peak is at the inner bar end. Then the cross-correlation peak will be at the radial lag at which these peaks overlap. But these regions are most likely unrelated. Though we restricted the radial lag space to avoid such extreme situations, one must treat the two dimensional cross-correlation results with caution.

In all six galaxies, the H$`\alpha `$ is leading the CO emission in the azimuthal direction. In other words, the correlation values are higher on the right hand side of the bottom panel. This is consistent with the overlays (§3.1) and the one-dimensional cross-correlation analyses (§3.2.1). We also find that, in almost all cases, there are higher cross-correlation values for H$`\alpha `$ offset radially outwards from the CO. In four galaxies (NGC 2903, NGC 3627, NGC 6946 and IC 342), the maximum of the cross-correlation values is located at outward radial lags of $``$7-10$`\mathrm{}`$. These large offsets are not obvious in the overlays of CO and H$`\alpha `$ images, and may not be physically meaningful because, as discussed above, the correlation is between H$`\alpha `$ and CO peaks along the dust lane.

## 4 Discussion

How can one explain the CO-H$`\alpha `$ offset? If there was no relative motion of the gas and stars, then the CO and H$`\alpha `$ should be co-spatial. However, since the majority of the H$`\alpha `$ emission is preferentially on the leading side of CO in all six bars, gas dynamics must play a role. Hence, we consider the two main classes of bar gas flow models: a) N-body/sticky particle models which treat the gas as a collection of sticky particles, reacting to a gravitational field, and exchanging angular momentum and energy upon collisions (e.g., Combes & Gerin 1985), and b) grid-based hydrodynamic models which treat the gas as an isothermal, ideal fluid obeying standard hydrodynamic equations (e.g., Athanassoula 1992a, b; Piner, Stone & Teuben 1995). We use the framework provided by these models to understand the observed CO-H$`\alpha `$ offsets.

### 4.1 An Important Caveat: Detection Limits on the Molecular Gas Data

The addition of single dish (NRAO 12m) On-The-Fly data to the BIMA interferometric observations greatly reduced the problem of spatial filtering. As discussed by Helfer et al. (2002a) and Helfer et al. (2002b), the combined data, such as those presented here, recovers 80-90% of the flux density for all sources of size $`<`$ 10$`\mathrm{}`$, with the recovery improving for smaller scales. At the mean distance to these six galaxies (8.25 Mpc), 10$`\mathrm{}`$ subtends a linear scale of 400 pc, much larger than a typical giant molecular cloud ($``$40–50 pc, Scoville 1990), and larger than most giant molecular complexes or associations ($``$100-400 pc, Rand 1993). So these data do not suffer greatly from the spatial filtering of an interferometer.

The sensitivity of these observations is discussed in detail in Helfer et al. (2002b). We summarize the main points here. Given our typical rms noise of 58 mJy beam-1, we have a rms column density sensitivity of 4.6 $`M_{\mathrm{}}`$ pc-2, using a CO/H2 conversion factor of 2$`\times `$1020 cm-2 (K km s-1)-1 in a 10 km s-1 channel. The data quality is such that the edges of extended emission are probably reliable at 2$`\sigma `$ or 9.1 $`M_{\mathrm{}}`$ pc-2, and isolated points are reliable at 3$`\sigma `$ or 13.7 $`M_{\mathrm{}}`$ pc-2. For comparison, the mean surface density of a Galactic giant molecular cloud in the solar neighborhood is 50–100 $`M_{\mathrm{}}`$ pc-2 (Blitz, 1993), within an area of 2.1 $`\times `$ 103 pc2. Assuming that this surface density is typical for molecular gas, we are sensitive to clouds with typical masses $``$ 106 $`M_{\mathrm{}}`$ in a single beam, or in other words, we would detect large molecular clouds in a single beam but our observations would not detect a 105 $`M_{\mathrm{}}`$ cloud. Hence it is possible that we are not detecting molecular gas from small, weak clouds that may be associated with the observed HII regions. This caveat applies throughout our discussion.

Our analysis is focussed on the observed offset between the CO and H$`\alpha `$ emission. The basic assumption is that these observations trace the highest surface density CO emission, i.e., the regions most prone to star formation activity. In the next two sections, we discuss how star formation may occur given the gas flow in the N-body/sticky particle and hydrodynamic models.

### 4.2 Gas flow and Star Formation in the N-body/sticky Particle Models

The N-body/sticky particle simulations (e.g., Combes & Gerin 1985) are designed to emulate the behavior of bound, molecular clouds, the basic organizational unit of the molecular interstellar medium. However, these models parameterize the physics involved by using ad hoc equations for simulating the collisions and energy exchange. Though they have been successful at reproducing curved dust lanes and rings (Combes, 1996), they are unable to reproduce the sharp velocity jumps and straight dust lanes observed in strongly barred spirals.

In these models, the bar dust lane results from orbit crowding, similar to a spiral arm dust lane (see Figure 10 in Regan, Sheth & Vogel 1999). Gas clouds in these models crowd together and probably form giant molecular complexes in the dust lane; they eventually leave the dust lane diverging outwards on the leading side of the bar. If stars form in the molecular complexes, as suggested for the spiral arm of M51 (Vogel, Kulkarni & Scoville, 1988), then we naturally expect H$`\alpha `$ emission to be both in the dust lane and downstream from it. Thus, these models naturally fit the observed CO-H$`\alpha `$ offset seen in our sample. They may also explain the presence of H$`\alpha `$ emission on the trailing side of the dust lane by arguing that, as the clouds gradually converge on the trailing side they may collide, form complexes, and form stars. Also as noted by (Vogel, Kulkarni & Scoville, 1988), the molecular gas clouds are difficult to detect on the leading (downstream) side of the dust lane because the average gas surface density decreases as the clouds diverge outwards from the dust lane. These models cannot easily explain the radially outwards H$`\alpha `$ offset seen in the two dimensional cross-correlation diagrams because, over most of the bar, the cloud orbits have a radially inwards component as clouds converge towards the dust lane.

### 4.3 Gas flow and Star Formation in Hydrodynamic Models

In contrast to the N-body/sticky particle models, hydrodynamic simulations use standard fluid equations to model the gas but the gas is usually an isothermal, ideal fluid. Though self-gravity and magnetic fields have not been included in the past, newer models have begun to incorporate these as well (e.g., Kim & Ostriker 2002). In contrast to the N-body/sticky particle models, the hydrodynamic models have been successful, not only at reproducing the observed gas kinematics (Regan, Vogel & Teuben, 1997; Regan, Sheth & Vogel, 1999), but also the entire range of observed bar dust lane shapes, from curved to straight (Athanassoula, 1992b).

The gas kinematics in these models are dramatically different (see Figure 9 in Regan, Sheth & Vogel 1999). The gas streamlines diverge as they approach the dust lane (Regan, Vogel & Teuben, 1997). At the dust lane, all of the gas is redirected inwards by a hydrodynamic shock such that none of the gas crosses the dust lane. Historically, the bar dust lanes, despite their high gas density, are regarded as inhospitable environments for star formation because of high shear (Athanassoula, 1992b, and references therein). Regan, Vogel & Teuben (1997) have also argued that the diverging streamlines on the trailing side of the dust lane can tear apart molecular clouds, or at least prevent molecular complexes from forming (but note that self-gravity can counteract this effect, see discussion in Sheth et al. 2000). These studies argue effectively for reduced star formation efficiency in the dust lanes.

Even if stars form in the dust lane, the gas flow dictates that stellar clusters travel down the dust lane. Thus these models cannot naturally explain the azimuthal offset between the H$`\alpha `$ and the CO. In this model gas flow, the only way for the HII regions to appear on the leading side of the dust lanes is to form stars on the trailing side while the gas still has tangential motion.

Such a mechanism was proposed in a previous study of the barred spiral NGC 5383. In that bar, Sheth et al. (2000) noted that the limited number of HII regions, all on the leading side of the dust lane, were directly across from faint dust spurs. They proposed that, in this bar, stars were forming in spurs because spurs were regions of high gas density and low shear. Spurs are faint and difficult features to detect. Though we see some spur-like CO features, the correlation between HII regions and spurs requires more sensitive CO observations. A more practical approach would be to first trace the spurs with high quality (resolution and sensitivity) optical-infrared color maps to understand their locations, structure, and frequency. Then a targeted study with sensitive CO observations can help us better understand the relationship between spurs, CO and H$`\alpha `$ . Though possible, we feel that it is unlikely that all star formation in bars occurs in spurs. However we note that, if stars form in the spurs, then the radially outwards offset can be naturally explained because the gas streamlines diverge outwards on the trailing side of the dust lanes in this model.

In summary, both gas flow models can explain the CO-H$`\alpha `$ offsets but neither is satisfactory, especially when we consider these and previous observations of bars. One possible improvement, especially to the hydrodynamic models, could be the use of a multi-phased molecular medium, consisting of a diffuse, gravitationally unbound phase (e.g., Elmegreen 1993; Rand, Lord & Higdon 1999; Hüttemeister et al. 2000), and the usual dense, bound phase made of giant molecular clouds (e.g., Scoville 1990; Young & Scoville 1991). In a bar, diffuse gas may be formed by the disruption of clouds by the tidal field of a bar, or by off-center cloud-cloud collisions (Hüttemeister et al., 2000). At least in one bar, NGC 7479, tentative evidence of diffuse gas has been presented (Hüttemeister et al., 2000). The diffuse gas has a higher velocity width than the bound phase, and so it may dominate the CO kinematics. Perhaps this is why studies of gas kinematics in the dust lane have been consistent with the hydrodynamic models (e.g, Regan, Vogel & Teuben 1997; Regan, Sheth & Vogel 1999). And it may be that the star forming component of the gas, the dense, bound clouds, behaves more like the sticky particles in N-body simulations.

### 4.4 Time Scales and Cloud Speeds in Context of Model Gas Flows

Previous spiral arm studies found H$`\alpha `$ offset from the CO by a few hundred pc (Vogel, Kulkarni & Scoville, 1988; Rand, 1993; Knapen & Beckman, 1996; Loinard et al., 1996). Vogel et al. (1988) immediately recognized that this offset was too large given that typical “drift” velocity of a cloud across a spiral arm is $``$ 10 km s-1. They suggested that stars may be forming after a “gestation” period of a few million years following the compression of clouds by the spiral density wave.

The situation is different in the bar environment. In the N-body/sticky particle model for bars, the orbit crowding is much more severe than for a spiral arm shock, and it is non-trivial to calculate an equivalent “drift” velocity across the bar dust lane. The offsets observed along most of the bar dust lanes (0–500 pc, §3.2.2) are similar to those observed along spiral arms. Therefore if the drift velocity in bars is similar to that in spiral arms, similar “gestation” periods are necessary in the framework of the N-body/sticky particle gas flow model. However in one or two locations along bars in NGC 2903 and NGC 3627, we measure offsets as large as $``$ 800 pc. Then for the star formation scenario presented here in the framework of N-body/sticky particle gas flow model, either unusually large gestation periods and/or large drift velocities are required. One observational test may be to measure speeds of newly formed clusters using stellar absorption lines to better understand the dynamics. Another possibility is that for these particular cases, our basic assumption that these stars formed in the dust lanes is not applicable (see caveat in §4.1).

In the framework of the hydrodynamical models there is no “drift” velocity because all of the gas is redirected inwards by the shock. Suppose the newly formed stars inherit the speed of their parent molecular clouds. Then to traverse the entire range of offsets seen in these bars (0–800 pc, §3.2.2), the velocity for a stellar cluster must be 0–80 km s-1, assuming a time scale of 10 Myr for the H$`\alpha `$ peaks. Individual HII regions have lifetimes $``$ 3 Myr (Spitzer, 1978) but since we do not resolve individual HII regions at these distances, we assume a larger lifetime for the H$`\alpha `$ peaks because they are most likely a collection of HII regions.

The relative speed of clouds entering the bar pattern is high. Though difficult to quantify, a few direct measurements of bar pattern speeds indicate values between 50–65 km s-1 kpc-1 (Merrifield & Kuijken, 1995; Dehnen, 1999; Gerssen, Kuijken, & Merrifield, 1999). Hydrodynamic models which try to reproduce observations use a bar pattern speed of 30–35 km s-1 kpc-1 (Piner, Stone & Teuben, 1995), but these bars are longer than typically observed. In any case, both the models and observations indicate that bars end at $``$80% of their co-rotation radius, consistent with the findings of Athanassoula (1992a, b).

Using this empirical relationship, and deprojected bar lengths (Table 5 from Sheth et al. 2002), the bar pattern speeds in 5/6 galaxies were measured using CO rotation curves by Das et al. (2001) to be: 68 km s-1 kpc-1 (NGC 2903), 55 km s-1 kpc-1 (NGC 3627), 35 km s-1 kpc-1 (NGC 4321), 71 km s-1 kpc-1 (NGC 5457), and 55 km s-1 kpc-1 (NGC 6946). These values are consistent with the direct measurements of bar pattern speeds in other galaxies. The pattern speed of the bar in IC342 is difficult to estimate because the bar parameters are not well-known. We estimated a bar length of $``$58$`\mathrm{}`$ (0.6 kpc). For this bar length, the pattern speed from the CO rotation curve is 196 km s-1 kpc-1, an unrealistically high value. If instead we use a length of 4$`\mathrm{}`$ as suggested by Crosthwaite et al. (2001), and assume that the rotation curve remains flat at $``$ 141 km s-1, then the bar pattern speed would be a more reasonable 81 km s-1 kpc-1. We do not use IC 342 in the ensuing discussion.

The speed at which a molecular cloud enters is the difference between the disk and bar pattern speeds. Hence, at the mid-points of these bars, cloud speed varies from 49–99 km s-1. So in the context of the hydrodynamical models, if a stellar cluster forms from clouds upstream of the shock, it would inherit a velocity of this magnitude. The velocity of the clouds upstream of the shock is primarily tangential to the bar dust lane. Hence, while the gas is redirected inwards down the dust lane, newly formed stellar clusters continue on the original orbits, and the observed range of offsets is consistent with the expected cloud speeds. However, we emphasize that the range of CO-H$`\alpha `$ offsets in a bar is not correlated with the galacto-centric radius. This suggests that star formation is not a simple function of the parent cloud speed, or alternatively a function of the shock strength. Consistent with this conclusion is the range in offsets observed in a single bar which suggests that star formation depends on other factors, e.g., gas surface density, and/or local shear.

## 5 Conclusions

We have investigated the distribution of molecular gas and star forming regions in the bars of six spirals from the BIMA Survey of Nearby Galaxies. Our main conclusions are as follows:

1. The CO emission is brightest along the leading edges of the bar. Weak spurs of CO emission are seen on the trailing side of the dust lanes. These spurs may play a role in star formation upstream of the dust lane. At the bar ends, strong CO and H$`\alpha `$ emission are seen on the trailing and leading sides; these may be the beginnings of inner rings.

2. The H$`\alpha `$ emission is distributed in compact and diffuse structures. There are a few instances where the H$`\alpha `$ is coincident with, or on the trailing side of the CO emission. But the main result is that in all six cases, the majority of the H$`\alpha `$ emission is offset towards the leading side of the CO.

3. We quantify the offsets using a cross-correlation analysis and find a range of 0–800 pc, with larger offsets in stronger bars. However, in a given bar there is a range of offsets and there is no systematic pattern as a function of the galacto-centric radius.

4. In the two dimensional cross-correlation analysis, there is a tendency for the H$`\alpha `$ emission to be offset radially outwards from the CO emission. However, this trend is less significant than the azimuthal offsets because correlations in two-dimensions can be between physically unrelated HII and CO regions.

5. The observed CO-H$`\alpha `$ distributions may be explained by either the N-body/sticky particle models or the hydrodynamic models with different, plausible, star formation scenarios. In the context of the N-body simulations, stars may form via cloud-cloud agglomeration in the dust lanes. In the context of the hydrodynamic models, the stars could form in dust spurs on the trailing side of the dust lane. We suggest that addition of a two-phased or multi-phased molecular medium can improve the agreement between these and previous observations, and gas flow models in bars.

Acknowledgements

This work would not have been possible without the rest of the SONG team members (T. Wong, T. Helfer, L. Blitz and D. Bock) and the dedicated observatory staff at Hat Creek and at the Laboratory for Millimeter-wave Astronomy at the University of Maryland. We thank S. Aalto, S. Hüttemeister, J. Kenney, E. Ostriker, N. Scoville, E. Schinnerer, J. Stone, and T. Treu for invaluable and insightful discussions about gas kinematics and star formation. We are grateful to N. Reddy and E. Schinnerer for their careful reading and helpful comments which significantly improved this paper.

Research with the BIMA array is supported by NSF grant AST-9981289. Support for the Laboratory for Millimeter-wave Astronomy is also provided by the state of Maryland. This research is also partially funded by NSF grant AST-9981546.


format: text

text:
# Ising magnets with mobile defects

## I Introduction

Low–dimensional magnetism in high-$`T_c`$ superconductors has attracted much interest, both theoretically and experimentally Review . In particular, striped structures in magnets derived from the $`La_2CuO_4`$ compound have been discussed rather extensively. Motivated by related analyses and, more specifically, by recent experiments Buch  on $`(Sr,Ca,La)_{14}Cu_{24}O_{41}`$, we shall introduce a novel class of quite simple two–dimensional Ising models, mimicing $`Cu^{2+}`$ ions by spin-1/2 Ising variables and holes by non–magnetic defects ($`S`$= 0). Of course, the aim of the present study is not to offer a full, or partial explanation of the experimental subtleties. Indeed, beyond the experimental motivation, the model is hoped and believed to show various intriguing properties being of genuine theoretical interest.

locatWe consider the situation where the spins are arranged in chains, with antiferromagnetic interactions, $`J_a`$, between adjacent chains, and a ferromagnetic coupling, $`J`$, between neighbouring spins in the chains, augmented by an antiferromagnetic coupling, $`J_0`$, between next–nearest spins in the same chain with a defect in between them. The defects are allowed to move through the crystal, along the chains.

The defects tend to form stripes, perpendicular to the chains, which, for increasing temperature, first meander and then become unstable. To identify the impact of the defect mobility on the stripe instability, a ’minimal model’ is proposed by assuming indefinitely strong interactions in the chains, $`J`$ and $`J_0`$. This model is studied analytically, at low and high temperatures, and, for a wide range of temperatures, by using standard Monte Carlo techniques Binder . Contact will be made to well known descriptions of wall instabilities in two dimensions, as have been put forward, for instance, in the context of incommensurate superstructures in two dimensions Pokro  and step roughening on vicinal surfaces Villain . Deviations from these standard scenarios will be discussed.

To study the robustness of the properties of the minimal model and to identify other possibly interesting aspects of this class of Ising models as well, we also considered cases with finite couplings in the chains. In addition, both for the minimal model and the ’full model’, the effect of an external magnetic field has been investigated.

The paper is organized accordingly. In the next section, we introduce the model and elucidate its experimental background. Then, we present our results on the minimal model, followed by a discussion on properties of the full model. A short summary concludes the article.

## II Models and methods

We consider Ising models on a square lattice, setting the lattice constant equal to one. Each lattice site $`(i,j)`$ is occupied either by a spin, $`S_{i,j}=\pm 1`$, or by a defect corresponding to spin zero, $`S_{i,j}=0`$. The defects are assumed to be mobile along one of the axes of the lattice, which will be called the chain direction in the following. Neighbouring sites of the same chain, $`(i,j)`$ and $`(i\pm 1,j)`$, are not allowed to be occupied both by defects (short range repulsion between defects). The mobility of a defect may be influenced by a pinning potential, which, however, will be disregarded in most of the following analysis. Even for vanishing pinning, the defect cannot diffuse freely, in general, because an elementary move, characterized by exchanging a defect and a, possibly flipped spin at neighbouring sites, is affected by the magnetic interactions along and perpendicular to the chain direction. We assume a ferromagnetic coupling, $`J>0`$, between neighbouring spins, $`S_{i,j}`$ and $`S_{i\pm 1,j}`$, along the chain augmented by an antiferromagnetic interaction, $`J_0<0`$, between those next–nearest spins of the same chain, which are separated by a defect. Spins in adjacent chains, $`S_{i,j}`$ and $`S_{i,j\pm 1}`$, are coupled antiferromagnetically, $`J_a<0`$. Accordingly, the Hamiltonian of the model may be written as

$`={\displaystyle \underset{ij}{}}(JS_{i,j}S_{i\pm 1,j}+J_0S_{i,j}S_{i\pm 2,j}(1S_{i\pm 1,j}^2)`$
$`+J_aS_{i,j}S_{i,j\pm 1})H{\displaystyle }_{ij}S_{i,j}`$ (1)
where we included a field term. Note that the defects, $`S_{ij}=0`$, are separated, along the chains, by at least one spin. We shall assume that the number of defects is the same in each chain, determined by the defect concentration $`\mathrm{\Theta }`$, denoting the total number of defects divided by the total number of sites.

The model describes, among others, the thermal distribution of defects, leading to rather intriguing properties, as will be shown below. In particular, at low temperatures, the defects tend to form stripes perpendicular to the chains which become unstable at higher temperatures. The model is believed to be of genuine theoretical interest.

The theoretical model may be motivated by rather recent experimental findings for the so called telephone number compound $`(Sr,Ca,La)_{14}Cu_{24}O_{41}`$ which contains two magnetically one-dimensional elements. One subsystem is a sheet like arrangement of $`Cu_2O_3`$ two-leg ladders, which is irrelevant in the context of the present paper. The second subsystem is an array of $`CuO_2`$ chains formed by edge sharing $`CuO_4`$ plaquettes. For this bond geometry with a $`CuOCu`$ bond angle close to 90 degree the Goodenough-Kanamori-Anderson rules predict a ferromagnetic exchange between nearest neighbour $`Cu`$ ions with spin $`S=\frac{1}{2}`$. This is confirmed by neutron diffraction studies of the magnetic structure in the ordered state which, in addition, show an antiferromagnetic coupling of the spins perpendicular to the chains Matsuda98 . To our knowledge the absolute value of the ferromagnetic coupling constant $`J`$ in the undoped chains has not been determined yet from inelastic neutron data. A mean field treatment of the magnetic susceptibility suggests a coupling constant of several meV Carter ; Buch .

A long range magnetic order of the $`Cu`$ spins is only observed for certain compositions of the telephone number compound. For many compositions the chains contain a large number of hole-like charge carriers. These holes imply non-magnetic $`Cu`$ sites Carter ; Ammerdimer ; KataevPRB  which inhibit the formation of a long range ordered magnetic state. For example, in the stoichiometric compound $`Sr_{14}Cu_{24}O_{41}`$ about 60 percent of the $`Cu`$ sites in the chains are non-magnetic Osafune97 ; Regnault ; Nucker00  and the remaining spins form nearly independent dimers Regnault ; Ammerdimer . The analysis of this dimer state shows a rather large antiferromagnetic coupling $`|J_0|11`$ meV between the two $`Cu`$ spins adjacent to a non-magnetic $`Cu`$ site Regnault ; Ammerdimer ; Remark1 . This coupling is about one order of magnitude larger than the antiferromagnetic coupling between $`Cu`$ ions in adjacent chains Regnault .

The experimental data mentioned so far yield at least three relevant magnetic coupling constants in the telephone number compounds: a ferromagnetic coupling $`J`$ of several meV between nearest neighbour $`Cu`$ ions in the chain, the antiferromagnetic $`|J_0|11`$ meV for $`Cu`$ spins adjacent to holes Remark1  and an antiferromagnetic coupling $`J_a`$ between $`Cu`$ spins in adjacent chains which is of the order of 1 meV Regnault . For undoped chains the magnetic properties depend mainly on $`J`$ and $`J_a`$, whereas the behaviour for large hole content is determined by a single coupling constant, the antiferromagnetic exchange $`J_0`$. For small hole concentrations all three magnetic interactions and their interplay should be relevant. Experimentally such a situation is realized in $`La_5Ca_9Cu_{24}O_{41}`$ where the hole content in the chains amounts to about 10 percent Nucker00 . Studies of this compound reveal a very unusual suppression of the magnetic order in external fields which can not be explained in terms of conventional spin models Buch . It is tempting to attribute this strange behaviour to a magnetic field induced movement of the charge carriers which frustrates the antiferromagnetic interchain coupling.

As will be shown below our numerical results for the model Eq. (1) indeed reveal a movement of holes in external fields. We mention that the treatment within an Ising model is also related to experimental findings for the telephone number compound. Different experimental data for lightly doped chains show a strong Ising-like anisotropy Buch ; KataevPRL  which was predicted by two independent theoretical treatments for spin chains formed by edge sharing $`CuO_4`$ plaquettes Aha .

To study the above Ising model with mobile defects, Eq. (1), we applied analytical low and high temperature considerations, and Monte Carlo techniques monitoring various thermodynamic and microscopic properties. In the simulations, we took into account flips of single spins as well as hops of a defect to a neighbouring site in the chain leaving a spin at the former defect site. Of course, simulations are performed on finite lattices with $`L\times M`$ sites ($`L`$ refers to the chain direction). Usually, we employed full periodic boundary conditions. In a few selected cases free boundary conditions perpendicular to the chain direction were applied, corroborating the results for periodic boundary conditions to be presented in the following. To investigate finite size effects, the linear dimensions, $`L`$ and $`M`$, were varied from 20 to 160. Typically, runs of at least $`10^6`$ Monte Carlo steps per spin were performed, averaging then over a few of such realizations to estimate error bars. The concentration of defects, $`\mathrm{\Theta }`$, ranged from zero to 15 percent. In most cases, we set $`\mathrm{\Theta }=0.1`$.

Physical quantities of interest include the specific heat, $`C`$, determined from energy fluctuations and the temperature dependence of the energy, the magnetization per site, $`m`$, and the correlation functions parallel, $`G_1(r)=(\underset{ij}{}S_{i,j}S_{i+r,j})/LM`$, and perpendicular, $`G_2(r)=(\underset{ij}{}S_{i,j}S_{i,j+r})/LM`$, to the chain direction. We also calculated other microscopic quantities describing the stability of the defect stripes and the ordering of the spins and defects in the chains. In particular, we computed the average minimal distance, $`d_m`$, between each defect in chain $`j`$, at position $`(i,j)`$, and those in the next chain, at $`(i^{},j+1)`$ (i.e. $`\mathrm{min}|ii^{}|`$, dividing this sum over the defects by their number), the cluster distribution, $`n_d(l)`$, denoting the probability of a cluster with $`l`$ spins of equal sign in a chain (in analogy to the distribution of cluster lengths in percolation theory Stauffer ), and the normalized number of sign changes of neighbouring spins, $`n_c`$, in the chains. Finally, it turned out to be quite useful to visualize the microscopic spin and defect configurations as encountered during the simulation.

In case of the magnets mentioned above, the absolute values of both $`J_0`$ and $`J`$ are large compared to the coupling between chains, $`J_a`$. To describe then the behaviour of the Hamiltonian, Eq. (1), at low temperatures one may consider a simplified model in which the spins form intact clusters in the chains between two consecutive defects changing the sign at the defect, i.e. $`J`$ and $`|J_0|`$ are assumed to be indefinitely strong. Quantities can now be expressed in terms of $`k_BT/|J_a|`$. The thermal excitations are due to motion of the defects in the chains. Again, defects are separated by at least one spin. The analysis of this ’minimal model’ will be presented in the next section.

## III Properties of the minimal model

### III.1 Monte Carlo results

In the ground state of the minimal model with vanishing external field, $`H=0`$, the defects form straight lines perpendicular to the chains and separating antiferromagnetic domains of spins. The ground state is highly degenerate, with the degeneracy depending on the concentration $`\mathrm{\Theta }`$ of defects. Each arrangement of straight defect lines, with a separation distance between the lines of at least two lattice spacings, has the same, lowest possible energy, resulting in a fast decay of the correlations $`G_1`$ parallel to the chains, while the spins are perfectly correlated perpendicular to the chains. The degeneracy may be (partly) lifted, for instance, by introducing a pinning potential or by applying an external field as will be discussed briefly below.

Increasing the temperature, $`T>0`$, the defects are allowed to move so that the stripes start to meander and finally break up, as exemplified in typical Monte Carlo configurations depicted in Fig. 1. It seems plausible that the destruction of the defect stripes by thermal fluctuations is accompanied by singular behaviour of thermodynamic quantities, like the specific heat, and various correlations functions. This suggestion is, indeed, supported by the numerical evidence discussed below. The effect of both phenomena, meandering and breaking up of the stripes, on various physical quantities are shown in Figs. 2 to 6, for the case $`\mathrm{\Theta }=0.1`$. Note that in most of the figures we did not include error bars being, typically, not larger than the size of the symbols.

At low temperatures, deviations from the straight stripes may be characterised by kinks and kink–antikink pairs Pokro ; Villain ; Pers ; Selke. Actually, the minimal model resembles closely a terrace–step–kink (TSK) model describing step fluctuations on vicinal surfaces. The energies of the elementary excitations may be readily calculated. For instance, a kink with depth of one lattice spacing costs $`J_a`$, a kink–antikink pair, created by moving a single defect by one site away from the perfect stripe, costs $`2J_a`$, for a further diffusion of that defect by another lattice spacing away from the stripe an additional $`4J_a`$ is needed, etc. The fact that consecutive defects in a chain cannot be closer than the minimum distance of two lattice spacings leads to the well known phenomenon of ’entropic repulsion’ Mullins ; Pers  between meandering neighbour stripes. Due to the entropic repulsion, the meandering stripes tend to approach their average distance as given by the concentration of defects, $`\mathrm{\Theta }`$ (here, at $`\mathrm{\Theta }=0.1`$, the average distance is ten lattice spacings). This feature is seen, e.g., in the thermal behaviour of the cluster distribution, $`n_d(l)`$, where the distance between two neighbouring defects in a chain is, in the minimal model, equal to $`(l+1)`$. As shown in Fig. 2 for the moderate size $`L=M=40`$, at zero temperature, $`n_d(l)`$, calculated numerically by taking into account all degenerate ground states, decreases monotonically with the separation distance $`l`$. When turning on the temperature, the maximum of $`n_d(l)`$ moves towards equidistant spacing of defects, here $`l=9`$, and the shape of $`n_d(l)`$ may be approximated by a Gaussian or Wigner function, as has been discussed recently in the context of terrace width distributions for vicinal surfaces Einstein . Increasing the temperature furthermore, $`T\mathrm{}`$, the defects take random positions in the chains, and $`n_d(l)`$ acquires, of course, again the same Poissonian form as at $`T=0`$ for random distribution of defect lines. The destabilization of the stripes is indicated, for example, by a rather rapid decrease, with temperature, of the probability to find defects at their average spacing, i. e. $`n_d(l=9)`$ for $`\mathrm{\Theta }=0.1`$, see Fig. 5.

As depicted in Fig. 3, meandering and breaking up of the stripes may be also observed in the correlation function parallel to the chains, $`G_1(r)`$. Again, the behaviour at zero temperature has been determined numerically without difficulty, for fairly small system sizes, by averaging over all ground states. The correlations are seen to decay rapidly. The oscillations in $`G_1(r)`$, already hardly visible for $`L=60`$, as shown in Fig. 3, become less and less pronounced when enlarging $`L`$ at fixed $`\mathrm{\Theta }(=0.1)`$. The asymptotics of $`G_1(r)`$, for large $`L`$, may be determined analytically, as discussed in the following subsection. Raising now the temperature, $`T>0`$, the correlations first become stronger, reflecting the ordering tendency which favour equidistant stripes due to the entropic repulsion, and then decrease quite drastically due to the thermal destabilization of the stripes. In fact, the perpendicular correlations, $`G_2`$, fall off rather rapidly in the same range of temperatures. Finally, when $`T\mathrm{}`$, one encounters again the behaviour at zero temperature, with the defects at random positions.

The destruction of the defect stripes is detected directly in the average minimal distance between defects in adjacent chains, $`d_m`$. Obviously, $`d_m`$ is equal to zero at $`T=0`$, and $`d_m2\mathrm{exp}(|J_a|/k_BT)`$ at low temperatures, $`k_BT|J_a|`$. In Fig. 4, data for various system sizes, $`L=M`$ ranging from 20 to 80, are displayed. While at low temperatures, $`d_m(T)`$ does not depend, in fact, significantly on the system size, it starts to rise rapidly at some characteristic temperature, with the height of the maximum in the temperature derivative of $`d_m`$ increasing strongly with larger system size. The location of the maximum, at $`T_d^{max}`$, signalling the breaking up of the stripes, moves to lower temperatures as $`L`$ gets larger. The quantitative behaviour is quite similar to the one of the specific heat, to be discussed below.

One possible reason for the destabilization of the stripes are effectively attractive interactions between neighbouring defects or lines, mediated by the spins. Indeed, such an interaction may occur, for instance, for strongly fluctuating stripes so that three consecutive defects in one chain, $`j`$, are in the cage formed by pairs of defects in the adjacent chains, $`j\pm 1`$. Consequently, two of the three defects tend to form a pair of next–nearest neighbouring defects, as may be checked easily. In any event, the probability of such pairs of defects is obviously given by $`n_d(l=1)`$. Its temperature and size dependence is depicted in Fig. 5 (together with that of $`n_d(l=9)`$, as mentioned above), showing a drastic increase close to the characteristic temperature of the breaking up of the stripes, $`T_d^{max}`$. Note that this type of stripe instability is not included in the standard descriptions of wall instabilities in two dimensions Pokro ; Villain ; Pers ; Halperin , where either the number of walls is not fixed, giving rise to incommensurate structures, or dislocations play a crucial role, in the context of melting of crystals. Also the bunching of steps in TSK models with attractive step–step interactions Weeks  or instabilities in polymer filaments due to attractive couplings Lipo ; Burk  are quite different from the loss of stripe coherency we observe here. Of course, the breaking up of the stripes has to be distinguished from their meandering which may result in their roughness, driven by capillary wave excitations Villain .

Meandering and destabilization of the stripes also show up in the specific heat, $`C`$, see Fig. 6 for systems with $`L=M`$ sites, $`L`$ ranging from 20 to 80. For each size, $`C`$ exhibits two maxima. The maximum at the lower temperature is almost independent of the system size, and it is related to the kink excitations of the stripes. The upper maximum, occurring at $`T_C^{max}(L)`$, signals the instability of the defect stripes. Its height increases with increasing system size, indicating possibly a phase transition in the thermodynamic limit, $`L\mathrm{}`$. To estimate the transition temperature, we plotted $`T_C^{max}`$ versus $`1/L`$, with $`L`$ going up to 160, see Fig. 7. From a linear extrapolation we obtain approximately $`k_BT_C^{max}(L=\mathrm{})/|J_a|=1.05\pm 0.05`$. Note that finite size analyses, usually for sizes up to 80, for other quantities, like $`G_1`$, $`G_2`$, and $`n_d(l=9)`$, lead to similar estimates for the possible transition temperature. However, the close agreement may be fortuituous, depending on the type of the transition. For instance, for a Kosterlitz–Thouless transition, the peak in the specific heat does not occur exactly at the transition temperature, as $`L\mathrm{}`$. A detailed analysis of this subtle feature is, however, beyond the scope of the present study.

From simulations of the minimal model with $`\mathrm{\Theta }=0.05`$ and $`\mathrm{\Theta }=0.15`$, $`L=M=40`$, we infer that the characteristic temperature, at which the stripes become unstable, gets smaller when the concentration of defects, $`\mathrm{\Theta }`$, is increased. This observation may be explained by the fact that the effectively attractive interactions, caused, for instance, by the cage effect described above, may set in obviously at lower temperature when the average distance between the stripes decreases.

Applying an external field, $`H>0`$, the ground states change when $`H`$ exceeds specific critical values. For $`|2J_a|>H>|J_a|`$, the stripes are no longer straight, but they form a zig–zag structure. In that structure, supposing the field favours the ’+’ spins, each ’+’ cluster comprises two more spins than the ’$``$’ clusters directly below and above that cluster in the two adjacent chains. Defects bounding these ’$``$’ clusters are located exactly below and above the first and last spins of the ’+’ cluster. Spins and defects in each second chain are arranged identically. Obviously, the zig–zag structures carry a non–vanishing net magnetization. The degeneracy of the ground state is still high, albeit somewhat smaller than in the case of straight stripes at $`H<|J_a|`$, because the minimum length of ’+’ clusters is now three, instead of one. For larger fields, $`H>2|J_a|`$, the ’$``$’ clusters shrink drastically: ’$``$’ spins occur only in the pairs of next–nearest neighbouring defects; all other spins point in the direction of the external field.

Monitoring various quantities at fixed defect concentration, $`\mathrm{\Theta }=0.1`$, the stripes are observed to become unstable at lower temperatures for stronger fields, $`0<H<2|J_a|`$. Varying the field at fixed small temperature, for the same defect concentration, we found, that the destruction of the stripes seems to be accompanied by a fairly rapid increase in the magnetization, $`m`$, leading to an anomaly in the field derivative of the magnetization, similar to experimental findings on $`(Sr,Ca,La)_{14}Cu_{24}O_{41}`$ Buch . The change from straight to zig–zag stripes, on the other hand, leads to a jump in $`m(H)`$ at $`H=|J_a|`$ and $`T=0`$. The corresponding maximum in the field derivative of the magnetization is, however, extremely weak already at very low temperatures. More detailed investigations are needed to clarify the experimental relevance of these observations. They are beyond our present scope.

The high degeneracy of the ground states may be lifted by introducing a pinning potential. For example, in related simulations for $`H=0`$, we found that meandering and breaking up of the stripes seem to be qualitatively not affected by a weak one–dimensional, along the chain direction, harmonic regular pinning potential. More realistically, one may introduce a random pinning potential. If it is sufficiently strong, it is expected to destroy the defect lines even at zero temperature. In a weak random potential, the coherency of the defect lines gets lost or the lines collide only on a large spatial scale, the ’Larkin length’ Larkin . At some finite temperatures, the collision length between neighbouring lines due to thermal fluctuations becomes smaller than the Larkin length. Starting from this temperature, the random pinning potential can be neglected in thermodynamics, though, it can be dominant for dynamic phenomena. Thus, our model is thermodynamically robust with respect to weak random potential except of a very low temperature region, in which probably a glassy state occurs.

### III.2 Analytical results

Here we derive the asymptotics of the spin correlation functions in the chains, $`G_1(r)=S_{i,j}S_{i+r,j}`$, for sufficiently large distances $`r`$, first at zero temperature, also valid at infinite temperature, and then at non–vanishing, but small temperatures. We start with the obvious statement that

$$S_{i,j}S_{i+r,j}=(1)^{n_0(i,i+r)}$$
(2)
where $`n_0(i,i+r)`$ is the number of defects or ’zeros’, in the interval $`(i,i+r)`$. Eq. (2) can be rewritten as

$$S_{i,j}S_{i+r,j}=(e^{i\pi n_0(i,i+r)}+e^{i\pi n_0(i,i+r)})/2$$
(3)
Assuming that $`n_0(i,i+r)=\overline{n}_0(r)`$ is sufficiently large, we apply the Gaussian statistics to the deviation $`\delta n_0(r)=n_0(i,i+r)\overline{n}_0(r)`$. Thus, the asymptotic expression of the correlation function reads

$$G_1(r)\mathrm{cos}\pi \overline{n}_0(r)\mathrm{exp}\left(\frac{\pi ^2}{2}(\delta n_0(r))^2\right)$$
(4)
The average $`\overline{n}_0(r)`$ is related to the concentration of zeros $`c_0=\mathrm{\Theta }`$ by $`\overline{n}_0(r)=c_0r`$. The calculation of $`\left(\delta n_0(r)\right)^2`$ is more tricky. First we calculate, at zero temperature, the probability $`p(n_0,r)`$ of a fixed value of $`n_0`$ at fixed $`r`$. We find, disregarding here the value of the minimal distance between defects,

$`p(n_0,r)=\mathrm{exp}(r(c\mathrm{ln}{\displaystyle \frac{c_0}{c}}+(1c)\mathrm{ln}{\displaystyle \frac{1c_0}{1c}}))/`$
$`\sqrt{2\pi rc(1c)}`$ (5)
with $`c=n_0/r`$. For large numbers, the probability (5) has the Gaussian form near the maximum leading to the final result:

$$(n_0\overline{n}_0(r)^2=rc_0(1c_0)$$
(6)
Plugging this expression into Eq. (4), we obtain the asymptotics of $`G_1(r)=S_{i,j}S_{i+r,j}`$ at zero temperature

$$G_1(r,T=0)\mathrm{cos}(\pi c_0r)\mathrm{exp}\left(\frac{\pi ^2}{2}c_0(1c_0)r\right)$$
(7)
The asymptotics is valid for $`rr_c`$, with the correlation length $`r_c=2/(\pi ^2c_0(1c_0)`$. The condition of having a minimal distance of $`q+1`$ spacings between consecutive defects can be taken into account by replacing $`c_{(0)}`$ by $`c_{(0)}/(1c_{(0)}q)`$. For our model $`q=1`$.

Note that the above considerations hold also in the high–temperature limit, $`T\mathrm{}`$, as mentioned before.

At finite, but small temperatures, $`k_BT|J_a|`$, the long–distance asymptotics of the correlation function $`G_1(r)`$ changes dramatically due to the meandering of the defect stripes. This process may be described by the free–fermion approximation Pokro . In this approach the lines are represented as trajectories of free fermions. Their entropic repulsion is treated as statistical repulsion of the fermions. The meandering of the lines means that, going from one moment of discrete time to the next one, each fermion can move to the left or right by one site with the amplitude (or probability) $`z=\mathrm{exp}(2|J_a|/k_BT)`$. These processes are described by the Hamiltonian

$$=z\underset{i}{}(a_{i+1}^+a_i+a_i^+a_{i+1})$$
(8)
where the fermion operators $`a_i`$ and $`a_i^+`$ obey cyclic bondary conditions $`a_{i+L}^{(+)}`$= $`a_i^{(+)}`$. The Hamiltonian is diagonalized by Fourier transformation

$`=2z{\displaystyle \underset{p}{}}(\mathrm{cos}p)\alpha _p^+\alpha _p;`$
$`p=2\pi m/L(m=1,2,\mathrm{}L)`$ (9)
where

$$a_k=\underset{p}{}e^{ipk}\alpha _p/\sqrt{L}$$
(10)
The energy band in $`p`$–space extends from $`\pi `$ to $`\pi `$, but it is filled only partly, from $`p_F`$ to $`p_F`$, where $`p_F=\pi c_0`$ (to take into account the specific minimal distance between defects, $`c_0`$ has to be substituted here and in the following as before). It means that in the ground state $`\alpha _p^+\alpha _p=1`$ for $`|p|=p_F`$, and $`\alpha _p^+\alpha _p`$= 0 for $`p_F<|p|<\pi `$. In this approach the number of zeros between $`i`$ and $`i+r`$ is given in terms of the fermion operators by

$$n_0(i,i+r)=\underset{k=i}{\overset{i+r}{}}a_k^+a_k$$
(11)
The average of this quantity is obviously equal to $`rc_0`$, as at zero temperature, but its variance, $`\delta n^2=(n_0(r)\overline{n}_0(r))^2`$, is drastically different. Indeed,

$$\delta n^2=\underset{k,k^{}=i}{\overset{i+r}{}}:a_k^+a_k::a_k^{}^+a_k^{}:$$
(12)
where :$`XY`$: denotes the normal product of the operators $`X`$ and $`Y`$. Applying the Wick theorem, one obtains

$$\delta n^2=\underset{k,k^{}=i}{\overset{i+r}{}}a_k^+a_k^{}a_ka_k^{}^+$$
(13)
with the same summation limits. The diagonal term of this sum, $`k=k^{}`$, gives the contribution $`rc_0(1c_0)`$, as at zero temperature. However, it will be completely compensated by the non–diagonal terms. Indeed, the simultaneous correlation function for free fermions is known LL  to be equal to

$$a_ka_k^{}^+=\mathrm{sin}(p_F(kk^{}))/(\pi (kk^{}))$$
(14)
It is easy to check that $`a_ka_k^{}^+=a_k^+a_k^{}`$. In the limit of small $`p_F`$, the summation in Eq.(13) can be replaced by an integration, which can be explicitly performed without difficulty under the additonal condition $`p_Fr1`$, giving

$$\delta n^2_{nondiag}=p_Fr/\pi +(\mathrm{ln}r)/4\pi ^2$$
(15)
The first term compensates the diagonal contribution. Thus, at $`0<k_BT|J_a|`$, the correlation function $`G_1(r)=S_{i,j}S_{i+r,j}`$ is finally approximated as

$$G_1(r,T=0)\mathrm{cos}(\pi c_0r)/\sqrt{r}$$
(16)
Of course, the algebraic decay of the correlations is in accordance with previous findings on free fermions in two dimensions Pokro ; Villain ; Pers . The presence of a pinning potential is expected to establish long–range order in the correlations. Formally, there is no continuous change from the exponential decay of the correlations, at $`T=0`$, to the algebraic decay at non–vanishing temperatures. However, the final expression, Eq. (16), is valid only for a system whose size, $`L`$, perpendicular to the stripes exceeds the collision length of the fermions, $`l_{coll}=1/(zc_0^2)`$, which goes to infinity as $`T0`$. At a fixed value of $`L`$, there exists a crossover temperature $`T_{cr}|J_a|/\mathrm{ln}(Lc_0^2)`$ at which the exponential decay of the correlations goes over into an algebraic one. Note that the size, $`L1/(zc_0^2)`$, is rather large at the concentration we mostly considered in the simulations, $`c_0=\mathrm{\Theta }=0.1`$, and the crossover effect plays no role there. However, such sizes are not big in experimental systems.

At large temperatures, the correlations $`G_1`$ are believed to decay exponentially, see Eq. (7) which also holds at infinite temperature. Thence one expects a transition from algebraic decay at low temperatures to an exponential decay at high temperatures, in accordance with the simulational results.

## IV Beyond the minimal model

In the following, we shall present results on the full model, Eq. (1), with finite ferromagnetic couplings, $`J`$, between neighbouring spins in a chain and antiferromagnetic couplings, $`J_0`$, between next–nearest spins in a chain separated by a defect. In the simulations we choose $`J_0/J=6.25`$ and $`J_a/J`$ ranging from zero to minus one. This choice is, again, motivated by the experimental findings mentioned above.

At zero temperature and small fields, one obtains the same highly degenerate ground states of perfectly straight or zig–zag stripes as in the minimal model. Likewise, the behaviour at low temperatures and $`H<2|J_a|`$ is characterised by the meandering of the stripes as in the minimal model, followed by the stripe instability at higher temperatures. However, the instability may be masked, for instance, in the specific heat at $`J_a/J=0.3`$, $`\mathrm{\Theta }=0.1`$, and $`H=0`$, for systems of sizes up to $`L`$= 80, as depicted in Fig. 8 and to be discussed in the following.

In this case, in addition to the weak, almost size–independent maximum at low temperatures due to stripe meandering, the specific heat displays a rather pronounced peak at higher temperatures being much stronger than the one in the minimal model. The peak, however, gets smaller when the system size increases. Indeed, it is non–critical, stemming from energy fluctuations by breaking bonds, $`J`$, between spins in the chains. It persists when setting $`J_a=0`$, i. e. in the one–dimensional limit exhibiting, of course, no phase transition at all. In that limit, the lower maximum in $`C`$ disappears, because there are no stripes. Actually, the defects and their mobility play an important role in breaking the bonds betwen spins in the chains. For instance, when a defect moves next to a flipped spin, the spin on the other side of the defect will be flipped rather easily, costing an energy of 2$`J`$ in the ferromagnetic coupling. In contrast, in chain without defects, an energy of 4$`J`$ is needed to create the elementary excitation comprising a pair of neighbouring spins with opposite signs.

We also analysed, for the case depicted in Fig. 8, the simulational data for the spin correlations, $`G_1`$ and $`G_2`$, the cluster distribution, $`n_d(l)`$, and the average distance between defects in adjacent chains, $`d_m`$. The data provide strong evidence that the defect stripes become unstable at about the same temperature, measured in $`k_BT/|J_a|`$, as in the minimal model. Thence the impact of the spin flips on the location of the stripe instability seems to be rather small for this choice of parameters. In principle, the spin flips may lead to new mechanisms, different from the effectively attractive interactions between the defects discussed for the minimal model, to destruct the coherency of the stripes. To detect the stripe instability in the specific heat at $`J_a/J=0.3`$, presumably significantly larger systems have to be simulated. As shown in Fig. 8, e.g., for $`L=M`$= 80 merely a shoulder in $`C`$ starts to develop, at about the temperature where the minimal model shows a peak in $`C`$, compare with Fig. 6. Note that very long runs, especially for large systems, are needed to get sufficiently good statistics for the Monte Carlo data.

Applying a magnetic field, $`H>0`$, the breaking up of the stripes is found to shift to lower temperatures, as in the minimal model.

When choosing a smaller, but non–vanishing ratio of $`J_a/J`$, one gets closer to the minimal model. We did simulations for $`J_a/J=0.1`$. In fact, there the instability of the defect stripes is also indicated by a maximum in the specific heat, as in the minimal model, already for rather small systems, e.g., $`L=40`$, followed by the large non–critical peak due to the spin flips in the chains.

On the other hand, when weakening $`J`$ with respect to $`J_a`$, the stripe instability, as indicated by the rapid increase of the minimal distance $`d_m`$, may occur quite close to the pronounced, non–critical maximum in $`C`$ being due to the spin flips in the chains. In particular, for $`J_a/J=0.5`$ and, especially, $`1.0`$, one then observes another clearly visible maximum nearby the non–critical peak in $`C`$ already for small and moderate system sizes, e.g. for $`L`$= 40, due to the stripe instabilty. The location of the instability, measured in units of $`k_BT/J`$, increases with increasing ratio $`J_a/J`$. Interestingly enough, near the instability, the probability of finding pairs of next–nearest neighbouring defects, $`n_d(l=1)`$, now does not show any longer the overshooting phenomenon, compared to complete disorder at $`T\mathrm{}`$, in contrast to the situation in the minimal model and for small ratios $`J_a/J`$, see Fig. 9. The stripes become unstable at temperatures of the order of the ferromagnetic spin coupling $`J`$, and then the tendency to form pairs of next–nearest neighbouring defects is diminished by thermal disordering.

## V Summary

In this paper a two–dimensional Ising model with defects being mobile along the chain direction has been introduced. Albeit the model has been motivated by recent experiments on cuprates with low dimensional magnetic interactions, the model is believed to be of genuine theoretical interest as well.

In particular, based on analytical, asymptotical considerations at low and high temperatures as well as on Monte Carlo techniques, the model is found to desribe formation of defect stripes, their thermal meandering and, at higher temperatures, their destabilization.

The meandering and the instability of the stripes is discussed in the framework of a minimal model, assuming infinitely strong couplings between the spins along the chains. The instability is signalled by pronounced anomalies in spin correlation functions, in the spin cluster distribution along the chain, in the specific heat, and in the minimum distance between defects in neighbouring chains. The breaking up of the stripes seems to be caused by an effectively attractive interaction between the defects mediated by the spins.

The main features of the stripe instability persist when replacing the infinite couplings by, presumbaly, experimentally more realistic values. However, the anomaly in the specific may be masked for rather small systems, and thermal disordering and spin flips may also reduce the pairing tendency of the defects.

New experimental data on a stripe instability in $`(Sr,Ca,La)_{14}Cu_{24}O_{41}`$, together with a more detailed discussion on possible theoretical interpretations, will be presented elsewhere.

###### Acknowledgements.

We like to thank T. W. Burkhardt, T. L. Einstein, and D. Stauffer for useful discussions. One of us, V.L.P., acknowledges the financial support from the Humboldt foundation and from NSF under grant DMR 00721115 as well as the kind hospitality at the Institute of Theoretical Physics of Aachen University.