text:
# Massive and luminous YSO IRAS 05361+3539 and its environment

## 1 Introduction

Massive YSOs are protostars (Palla & Stahler, 1993) which are either surrounded by ultracompact HII regions or will eventually be hot enough to produce HII regions. They are mostly found buried deep inside clouds of gas and dust. The major difficulties in the study of massive YSOs are (a) there are fewer massive YSOs compared to low-mass YSOs and most of them are at a distance greater than 1 kpc from the Sun and b) most massive YSOs suffer large extinction (Av $``$ 10) and hence are difficult to study in the optical wavelengths. However, they can be studied in the radio and infrared wavelengths. Further, it is possible to study the environment of these objects with seeing limited spatial resolution in the near-infrared using array detectors like NICMOS.

Recent studies by Churchwell (1997) revealed that massive YSOs undergo similar bipolar outflows like the low-mass YSOs, but the rate of mass outflow is larger by orders of magnitude ($``$3). Similarly luminosities of massive YSOs are also higher. The total mass locked in the outflow is often found to be larger than the central star. Churchwell (1997) proposed that the larger outflow mass could be due to the in-falling matter directly diverted into the bipolar jets. It is clear from earlier works (Shepherd & Churchwell 1996; Churchwell 1997; Hartquist & Dyson 1997; Hunter et al. 1997) that the impact of outflows from massive YSOs shapes the history of star-formation in the parent molecular cloud.

We have started a NIR observational program on studies of regions of massive star formation. This is our first paper in the series. In this paper we study the massive and luminous YSO IRAS 05361+3539 and the star formation in the neighborhood using 2MASS (The Two Micron All Sky Survey) data and new NIR observations from Mt. Abu, India. The region so far is not very well studied. IRAS source 05361+3539 (G173.58 +2.45) was studied in the millimeter lines (12CO and 13CO) by Shepherd & Churchwell (1996). The source is situated at a kinematic distance of 1.8 kpc (Wouterloot & Brand 1989). Earlier, Wouterloot et al. (1988) detected H2O maser from the source. The source is embedded inside a large molecular cloud. The FIR fluxes meet the conditions of Wood & Churchwell (1989) for UCHII regions and total FIR fluxes correspond to a B2.5 central star. Shepherd & Churchwell (1996) found bipolar flows in the 12CO velocity map with high velocities (up to -27.1 km/s in the blue shifted lobe) with a total mass of 32M locked in the outflow. They estimated from the IRAS fluxes the central star mass to be 7M.

In section 2 we present the observations and data reduction procedures, section 3 deals with results and discussion and we summarize our conclusions in section 4.

## 2 Observations and data reduction

### 2.1 Observations from Mt. Abu

IRAS 05361+3539 (hereafter IRS1) was observed from the 1.2 meter Infrared Telescope, Mt. Abu, India, using a 256$`\times `$256 HgCdTe array (NICMOS-3, made by Infrared Laboratories, Arizona, USA). The telescope is located at an elevation of 1700 meters from the mean sea level and the instrument is described in Nandakumar (2000). Both the telescope and the NIR camera are owned by the Physical Research Laboratory, Ahmedabad, India.

The source was observed in two sessions: in J, H and K’ bands on 10th January 2000, and in narrow band filters centered on 2.12$`\mu `$m (H2 $`v=10S(1)`$), 2.16$`\mu `$m (Br$`\gamma `$), and 2.14$`\mu `$m (continuum) on 25th February 2000. The FWHM of the narrow-band filters was 0.042$`\mu `$m. The plate scale was 1′′ / pixel for the broad-band images and 0.5′′/ pixel for the narrow-band images. The nights were photometric during the observations. The seeing was 2′′ during the broad band observations and 1.5′′ during the narrow band observations. A large number of dithered sky frames were obtained ( by shifting the telescope 3 arcmin off the source in north-south-east-west directions) in all the filters for sky subtraction and for making flat frames. Although the 2.14$`\mu `$m filter is supposed to exclude the nearby emission lines, we found from filter transmission curves that there is serious contamination from the 2.12$`\mu `$m line. We therefore will not consider the 2.14$`\mu `$m line image in this paper and will present the narrow band images as emission line plus continuum. We observed three standard stars (AS11, AS19, and AS21 from Hunt et al. 1998) during the observations. We used the relation given by Wainscoat & Cowie (1992) to obtain K’ magnitudes from K magnitudes.

The data reduction was done using IRAF software tasks. All the NIR images went through standard pipeline procedures like sky-subtraction and flat-fielding. Individual object frames were of 30s of integration in J and H band and 3s in K’ band. The images were co-added to obtain a final image in each band (J, H and K’) of total integration time of 150s. Further, in each band two such images were obtained. The images were further filtered using a median filter of 3$`\times `$3 pixels for removing noise at sky-level, thus making them suitable for photometric analysis. We used DAOPHOT (Stetson 1987) task for deriving the photometry in the unvignetted field of view of 2 arcmins. The zero point was obtained using the three observed standard stars. From the J, H and K’ images containing the IRAS source, stars were identified down to 3$`\sigma `$ level (peak signal), subsequently each frame was visually inspected at different contrast levels to cross check the detection and remove false detections by DAOFIND. Aperture photometry was performed on the images using PHOT task with aperture radius of 4.5 pixels in the J, H and K’ images. The sky was sampled using 5-pixel wide annuli centered on each star with inner radius at 5 pixels. The aperture size was decided using the brightest and isolated star. From the two images of 150s in each band, we found that the overall photometric error is of $`\pm `$0.07 mag. This estimation includes uncertainties in the determination of zero-point ($`\pm `$ 0.03 mag) from the observed three standard stars. We found that in the K’ images we could detect 12 stars up to 14.2 mag in the field of view of 2 using the above mentioned procedure.

The narrow-band images also went through similar image processing. The total integration time was 250s in H2 and Br$`\gamma `$ filters.

### 2.2 2MASS data

We extracted stars from the 2MASS point source catalogue which were within 8 diameter of the IRAS source. The data were downloaded from the 2MASS Homepage available free for the Astronomical community. The 2MASS observations were carried out on 3rd February 1998111The date of observation was obtained from the header of the 2MASS fits images.. The 2MASS point source catalogue consists of J, H and Ks magnitudes of stars. The Ks band (bandpass =2.00$`\mu `$m to 2.32$`\mu `$m with center at 2.17$`\mu `$m) is very similar to K band. A histogram of the 2MASS sources in the Ks band within a field of view of 8 around the IRAS source was plotted. The histogram showed that the completeness limit is close to 15 mag.

Since the 2MASS catalogue goes deeper than the Mt. Abu images we decided to use the catalogue magnitudes to plot color-color and color-magnitude diagram. However, the magnitudes obtained from the Mt. Abu images will be used to find prospective variable stars and to study the morphology of the IRAS source and associations.

Table 1 shows a comparison of J,H,K’ and Ks magnitudes of stars detected within the unvignetted field of 2 in the K’ image with the 2MASS point source catalogue. We found that for standard stars the difference in K’ and Ks magnitudes is less than 0.05 which is less than the photometric error. Hence we will treat the K’ and Ks magnitude scales as similar in this paper.

## 3 Results and Discussion

### 3.1 Star formation in the cloud

Figures 1 and 2 show the J-H/H-Ks color-color diagram and the Ks/H-Ks color-magnitude diagram of the sources from 2MASS data. In figure 1 the solid curve is the locus of points corresponding to unreddened main sequence stars (Koornneef, 1983). The two dashed lines are parallel to the reddening vector with magnitude of Av=30. They form the reddening band (drawn from the base and tip of the unreddened main sequence) and bound the region in which stars with normal photosphere fall (also see Hunter et al. 1995). Similarly, in Figure 2, the solid line represents the locus of unreddened main sequence and the dashed arrow shows the direction and magnitude of the reddening vector. We did a search in the SIMBAD database to find out any known sources around IRS1 with a search radius of 4.5. Apart from the IRAS source and its H2O maser source, a star of A0V type (HD37317, 7.79 mag in Ks) was found in the SIMBAD search. This star’s co-ordinates are 05h 39m 19.9s +35 38 30′′, that is approximately 2.8 south-west from the IRAS source. We found from the 2MASS point source catalogue the colors of the star to be H-Ks=0.05 mag and J-H=0.031 mag. The colors of IRS1 are H-Ks=0.93 mag and J-H=1.45 mag.

Stars in figure 1 can be divided into three main groups a) those lying on the left side of the reddening band, b) those lying in the reddening band and c) those lying on the right side of the reddening band. The stars in the first group (a) plotted as open circles can be further divided into two sub-groups: one, those having J-H color less than 2 and second, those with J-H color more than 2. The former subgroup are mostly foreground stars as supported by their low values of Av (less than 5 mag), and the stars in the latter subgroup could be spurious detections since they have J and H magnitudes fainter than 17 and 16 mag respectively. The second group (b) mostly consist of normal stars with low values of Av ($``$10) and background stars with high values of Av ($``$10) and are also plotted as open circles. The third group (c) contains stars showing excess emission in H and Ks. Such sources are mostly YSOs (Lada & Adams 1992; Lada et al. 1993; Gomez, Kenyon & Hartmann 1994). YSOs can be further divided into Class I, Class II and Class III type sources based upon their Spectral Energy Distributions (SEDs) (Strom et al. 1989; Kenyon et al. 1993; Hartmann 1998). We have plotted sources redder than IRS1 and falling on the right hand side of the reddening band as ‘filled circles’ and are considered to be Class I type or protostars (since IRS1 is a known protostar; also see Lada & Adams 1992). Sources which show low H-Ks color ($``$1.5) and J-H color ($``$1.0) and also lying on the right hand side of the strip of reddening lines are also plotted as ‘open circles’. All these sources are faint in the Ks band (14 to 15.2 mag). Even though most of them are within the limit of completeness, we will need deep K band images to verify the existence of the sources. Their Av values range from 13 to 23 mag, suggesting that they could be Class II type sources (Kenyon et al. 1993).

It is clear from the J-H/H-Ks color-color diagram that the region around IRS1 is undergoing a phase of star formation. Three out of five Class I sources detected (excluding IRS1) are found within one arcmin radius of IRS1. Two of these sources are brighter than 14.2 mag in K’ and are detected in the Mt. Abu images. A comparison of magnitudes of these two sources #$`9`$ and #$`12`$ (see Table1) show that one of them #($`12`$) has varied over the time of observation between 2MASS (February 1998, as given in the header of the FITS images) and Mt. Abu (January 2000). The star has become fainter by 1.1 mag in K’ and H bands in about 23 months. In the Mt. Abu J band image the star is not detected. Since the difference in magnitude is much larger than the photometric errors and the decrease in the brightness is consistent in K’ and H bands, we believe that this is real. The star also shows extreme reddening in the J-H/H-Ks color diagram, a typical characteristic of Class I sources and presence of circumstellar material. Variability in low-mass protostars is known (Hartmann 1998). The total luminosity depends on the mass accretion rate, and low mass protostars like FU Orionis type of objects have shown variability of several magnitudes in the optical wavelength (Bell et al. 1995, and references therein). It is possible that we have witnessed a FU Orionis kind of behavior from the star. However, we need further NIR observations to verify this.

### 3.2 Photometry and Spectral Energy Distribution of IRAS 05361+3539

IRAS 05361+3539 (IRS1) is a massive luminous YSO and is associated with an ultra-compact HII region (Shepherd & Churchwell 1996). The low-resolution NVSS contour plot shows emission towards IRS1 at mJy level($``$4)(Condon et al. 1998). The FIR IRAS fluxes of IRS1 at 100$`\mu `$m, 60$`\mu `$m, 25$`\mu `$m, and 12$`\mu `$m are 1310Jy (upper limit), 29.15Jy, 6.72Jy, and 1.18Jy respectively. We estimated the spectral index of IRS1 to be 1.1 between 60$`\mu `$m and 1.12$`\mu `$m (Jband). The source therefore belongs to the group of Class I type sources (Strom et al. 1989). The SED is shown as a plot (crosses) of Log$`(\lambda F)`$/Log$`(\lambda )`$ in figure 3 where $`\lambda `$ is the wavelength in $`\mu `$m and F is the flux in Watts/meter2. We determined the temperature distribution of the circumstellar matter by fitting a model (Anandarao, Pottasch, & Vaidya 1993) to the observed SED. We have assumed photospheric temperature of 20,000K, stellar radius equal to seven solar radii, and plane parallel geometry for the dust shells. The dashed line in Fig 3 shows the model. The derived dust temperatures are 800K at 4AU and 80K at 400AU from the central source. The uncertainties in the dust parameters could be as large as 10-20 % due mainly to the inherent non-uniqueness of the model and to some extent to the uncertainties in the assumed stellar parameters. These results confirm that IRS1 is a Class I type source. The range of dust shell parameters derived from the model seem to support the accretion disk scenario (Hartmann 1998; Adams, Lada & Shu 1987) in which case, the grain heating is due to two processes: one due to the reprocessing of UV photons and the other due to viscous heating in the disk (Hillenbrand et al. 1992). Following Hillenbrand et al. (1992) we estimate the minimum radius of an accretion disk attributing the entire 25$`\mu `$m flux to the accretion heating of the grain (i.e. T=140K). We derive radius of the disk to be 160 AU if the inclination is 0 (face on).

### 3.3 Morphology of IRAS 05361+3539 and Detection of a NIR Jet

Figure 4 shows K’, H2 (emission line + continuum) and Br$`\gamma `$ (emission line + continuum) images of IRS1. Figure 5 gives the contour map of H2 superposed on the K’ image. The plate scale in Fig 4 is 1′′/pixel for the K’ image and 0.5′′/pixel for the narrow band images.

The NIR images reveal (figures 4, 5) that the IRS1 source is associated with a nebulosity extended up to 5′′ in the northern direction and a filamentary structure of length 6′′ in the east (seen more prominently in the H2+continuum image).

The eastern filamentary structure bends beyond 7′′ from IRS1 and continues another 20′′ in the southeast direction. Since we do not have a continuum free H2 emission line image we cannot quantify the amount of H2 emission from the region. However, we qualitatively argue from the brightness in the K’ and the H2+continuum image that at least 40% of the total brightness in the narrow band image is due to pure H2 emission. Further, the filamentary structure matches with the eastward jet found in the low resolution CO map of Shepherd & Churchwell (1996). Therefore, it is likely that the H2 emission is tracing the jet from the YSO. The bending of the jet could be because of physical obstruction of its flow due to the presence of putative dense matter. We do not see significant structure in the Br$`\gamma `$ \+ continuum image. If one assumes a particle density of 105 /cm3, then there can not be sufficient flux of energetic UV photons available from a B2 star to ionize the matter except in optically thin regions.

The structure in the northern direction is six times brighter than the jet. There are three distinct possibilities regarding the nature of this structure. This could be an unresolved star closely associated with the IRS1 source. However, in the 2MASS point source catalogue there is no entry corresponding to that position. We will need high-resolution images in K band to prove that the northern structure is an unresolved star.

A second alternative is that the extended nebulosity is due to the presence of the ultra-compact HII region. From the near-infrared morphology, it appears that this could be a cometary UCHII region (Wood & Churchwell 1989). Supposing that the central star has a relative motion with respect to the parent molecular cloud, say a few km/s, such a supersonic motion can create a region of low density behind the star (Hughes & Viner 1976; Weaver et al. 1977). If we consider a dynamical time of 5.3$`\times `$104 years (Shepherd & Churchwell 1996) and a speed of 1km/s to 2km/s (Jones & Walker 1988) then, in as many years the star would have traveled a distance of 0.1 pc (upper limit). The extent of the extended nebulosity in the H2 \+ continuum image corresponds to about 0.04 pc ($``$5′′, see figure 4). The density behind the central star should be much less (103 particles/cm3) for the UV photons to reach up to 0.04 pc.

NIR spectroscopy of UCHII regions have shown that they are bright in Br$`\gamma `$ emission line (Armand et al. 1996; Doherty et al. 1994). Aγ for the UCHII region can be calculated from Av using the relation Aγ/Av=0.125 (for Rv=5.0 extinction law of Cardelli, Clayton & Mathis 1989), and it turns out to be 1.87 mag. The density of matter can also be calculated from the relation between Av and column density, Av/NH=1$`\times `$ 10-21 mag/cm2 (Cardelli, Clayton & Mathis 1989). We determined the density of matter within 1 pc of the source assuming that most of the extinction is from within 1pc of the source. This gives a lower limit of density of 5x103particles/cm3. The upper limit can be 105 particles/cm3, if we consider that most of the extinction is within 0.05 pc. The value of 0.05 pc has been assumed based on the fact that very close to the star up to 5000AU the density can be as high as 108 particles/cm3 and up to 0.5 pc it can be 105 particles/cm3 (Churchwell 1997). A sizable Stromgren sphere with the lower limit of particle density can be obtained for a central B2 type star. But, we do not see any significant brightness in the Br$`\gamma `$ \+ continuum image compared to the K’ band image except at the core of the nebulosity. A deeper Br$`\gamma `$ image will be required for dertermining the extent of the nebulosity.

As a third possibility, the extended nebulosity can be considered as a dense clump of molecular hydrogen in the parent molecular cloud. We found that the region is bright in the H2+continuum image and only qualitatively we can say that pure H2 emission is present. H2 emission can arise due to either UV fluorescence, or collisionally excited by the impact of stellar winds from IRS1 (Genzel, 1992).

From the above arguments, it appears that it is difficult to determine from the present data the nature of the extended nebulosity. High resolution radio continuum images and medium resolution NIR spectroscopy are necessary to resolve this issue.

## 4 Conclusion

The parent molecular cloud seems to be an active star forming region. We have found a total of 6 prospective Class I type sources in the region including IRAS 05361 +3539. A number of faint Class II type sources are also detected in the region. One of the Class I sources (star no.12) detected close to IRS1 shows extreme reddening in the J-H/H-Ks color-color diagram. This star also appears to be fainter in H and K’ band by 1.1 magnitude in the Mt. Abu images when compared with the 2MASS magnitudes. The time difference between these observations was 23 months. Therefore, star # 12 could be a variable protostar of FU Orionis type. We need further observations to prove it.

The infrared spectral index of IRS1 is estimated to be 1.1 from IRAS and NIR fluxes. Our model predicts the possibility of an accretion disk with dust temperatures from 80K to 800K and with an extent of several hundreds of AU. Traces of molecular H2 emissions have been detected along the CO jet axis mapped by Shepherd & Churchwell (1996) in the eastward direction. A nebulosity is also detected towards north extending up to 5′′ from IRS1. This could be due to either an unresolved star, an UCHII region, or a dense clump of molecular hydrogen.

###### Acknowledgements.

The present observations were made under a collaborative project between TIFR and PRL. We thank the staff at Mt. Abu IR Telescope facility for support during observations. We especially thank the referee Dr. Todd R. Hunter for his very useful comments on the paper which significantly improved the quality of the paper. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.


format: text

text:
# Evolution of Differentially Rotating Supermassive Stars to the Onset of Bar Instability

## 1 Introduction

The existence of supermassive black holes (SMBHs) and the theory of their role as fuel sources for active galactic nuclei and quasars are supported by a growing body of solid observational evidence (see, e.g., the reviews of Rees 1998 and Macchetto 1999, and references therein). For example, VLBI observations of the Keplerian motion of masing knots orbiting the nucleus of NGC 4258 indicate the presence of a central object of mass $`3.6\times 10^7M_{\mathrm{}}`$ and radius $`<13`$ pc (Miyoshi et al., 1995). The most reasonable explanation is that this extremely compact object is a black hole. High resolution HST observations of matter orbiting the nuclei of M87 and NGC 4151 also show nearly Keplerian motion around SMBHs (Macchetto et al., 1997; Winge et al., 1999). In addition, reverberation mapping and photoionization methods indicate the presence of SMBHs in the centers of 17 Seyfert quasars (Wandel, Peterson, and Malkan, 1999).

However, the nature of the progenitors of SMBHs is rather uncertain (see Rees 1984, 1999 for overviews). SMBH formation scenarios that involve the stellar dynamics of dense clusters and the hydrodynamics of supermassive stars have both been proposed. The suggested stellar dynamical routes focus on the evolution of dense star clusters. These clusters are formed via a conglomeration of stars, produced by fragmentation of the primordial gas. In one such cluster scenario, massive stars form via stellar collisions and mergers in the cluster and then evolve into stellar-mass black holes. The merger of these holes, as they grow and settle to the cluster center, leads to the build up of one or more SMBHs (see, e.g., Quinlan and Shapiro 1990, Lee 1995, and references therein). Alternatively, dense clusters composed of compact stars may become dynamically unstable to relativistic core collapse and thus may give birth to SMBHs (Zel’dovich and Podurets, 1965; Shapiro and Teukolsky, 1985a, b; Quinlan and Shapiro, 1987).

The proposed hydrodynamical routes to SMBH formation center on the evolution of supermassive stars (SMSs). SMSs may contract directly out of the primordial gas or they may develop from a supermassive gas cloud that has been built up from the fragments of stellar collisions in clusters (Sanders 1970; Begelman and Rees 1978; Haehnelt, Natarajan, and Rees 1998; however, see also Loeb and Rasio 1994; Abel, Bryan, and Norman 2000; Bromm, Coppi, and Larson 1999). The evolution of SMSs will ultimately lead to relativistic dynamical collapse (Iben, 1963; Chandrasekhar, 1964a, b) and thus possibly to the formation of SMBHs.

In addition to their relevance to AGNs and quasars, SMBHs and their formation events are likely candidates for detection by proposed space-based gravitational wave detectors, like the Laser Interferometer Space Antenna (LISA) (Folkner, 1998). LISA would be very sensitive to the long wavelength, low frequency gravitational radiation that supermassive objects are expected to emit (Thorne and Braginsky, 1976; Thorne, 1995). For example, LISA may be able to detect the collapse of a SMS to a SMBH or the coalescence of two SMBHs (LISA Pre-Phase A report, 1995). The rates of such events are largely unknown due to the uncertainty in the formation mechanism of SMBHs.

Identifying the scenario by which SMBHs form is hence of fundamental importance to a number of areas of astrophysics. This work examines one of the possible progenitors of SMBHs discussed above, SMSs. The structure and evolution of SMSs have been explored by many authors (see, e.g., Zel’dovich and Novikov 1971 and Shapiro and Teukolsky 1983, and references therein). In this paper, we focus specifically on the quasistatic evolution of SMSs as they contract and cool due to thermal emission.

Two related factors that affect the cooling phase of a SMS are the strengths of its viscosity and magnetic fields and the nature of its angular momentum distribution. If strong enough, viscosity will enforce uniform rotation throughout the star as it contracts. However, the strength of molecular viscosity is unlikely to be large enough to accomplish this, so turbulent viscosity or magnetic fields must provide the necessary dissipation (Zel’dovich and Novikov, 1971; Spruit, 1999a, b; Baumgarte and Shapiro, 1999). However, it is not known whether these agents are sufficiently strong to maintain uniform rotation in SMSs as they evolve.

If a SMS is rotating prior to its contraction phase, conservation of angular momentum requires the star to spin up and therefore increase its ratio of kinetic energy $`T`$ to gravitational potential energy $`W`$, $`\beta =T/|W|`$. If viscosity or magnetic fields are able to maintain uniform rotation, the star will be spun up to its mass–shedding limit $`\beta _{shed}=8.99\times 10^3`$. The value is small because a SMS is dominated by radiation pressure and thus is well represented by an $`n`$=3 polytropic equation of state (see §2.1.1). As a result, it is quite centrally and cannot support much rotation when spinning uniformly before mass at the equator is driven into Keplerian orbit. Thereafter, it will evolve along a mass–shedding sequence, losing both mass and angular momentum, and will contract to the point of relativistic instability and collapse. Baumgarte and Shapiro (1999) have recently investigated this evolutionary scenario in detail and have described how it may lead to SMBH formation. Specifically, they identified the critical configuration at the onset of the relativistic radial instability. They also demonstrated that the nondimensional (in units where $`G`$=$`c`$=1) ratios $`T/|W|`$; $`R/M`$, where $`R`$ is the radius and $`M`$ is the mass; and $`J/M^2`$, where $`J`$ is the angular momentum, are universal numbers (i.e., they do not depend on $`M`$, $`R`$, or $`J`$ for the crititcal configuration or on the history of the star).

This work investigates the extreme opposite regime of SMS cooling evolution, the low-viscosity and low-magnetic field limits. In this case, as the star contracts and spins up, neither viscosity nor magnetic fields will be able to enforce uniform rotation (even if the star initially rotates uniformly). The outcome of the evolution will then depend on the star’s initial angular momentum distribution (or rotation law), as the star will conserve its angular momentum on cylindrical shells as it contracts. One possible outcome is that the star will again be spun-up to the mass-shedding limit (at a value of $`\beta _{shed}`$ higher than that for uniform rotation). The subsequent evolution is not clear. It is likely that mass and angular momentum will be ejected from the equator as the star undergoes further cooling. It is possible that it will continue to evolve in quasistationary equilibrium until reaching the onset of dynamical collapse. If it implodes to a SMBH, the matter is likely to emit an appreciable burst of gravitational radiation just prior to entering the event horizon.

An alternative outcome of low-viscosity, low-magnetic field contraction of a differentially rotating SMS is that it will encounter the dynamical “bar mode” instability at a critical value $`\beta _{bar}`$ prior to reaching $`\beta _{shed}`$. Such global rotational instabilities in fluids arise from nonaxisymmetric modes, which vary like $`e^{\pm im\varphi }`$, where $`m`$=2 is the bar mode. As the bar mode develops, angular momentum and mass will be transported outwards (Tohline, Durisen, and McCollough, 1985; Pickett, Durisen, and Davis, 1996; Imamura, Durisen, and Pickett, 2000; New, Centrella, and Tohline, 2000; Brown, 2000).

If enough angular momentum is removed, collapse to a SMBH could ensue. A SMS undergoing a bar mode distortion would be a strong source of quasi-periodic gravitational radiation because of its highly nonaxisymmetric structure. We note that a bar might form during catastrophic collapse in any scenario, if $`\beta (1/R)`$ grows above the critical value for the dynamical bar mode instability.

The outcome of the cooling and the resultant spin-up of a low-viscosity SMS depends critically on its initial angular momentum distribution. However, the nature of the angular momentum distribution of SMSs is largely unknown. The purpose of this paper is to examine the evolution of SMS models by assuming several different initial rotation laws to determine, for example, if any of them might cause the star to encounter the dynamical bar mode, prior to mass-shedding and collapse (i.e., to determine if $`\beta _{bar}<\beta _{shed}`$).

To make this determination, we analyze several Newtonian equilibrium sequences of SMS models. The models along each sequence are constrained to have the identical mass $`M`$, angular momentum $`J`$, and angular momentum distribution but decreasing entropy, due to cooling. Such a sequence mimicks the quasistatic evolutionary track of a SMS as it slowly cools and contracts. We find that along the sequence, the stars flatten as measured by the ratio of polar radius $`R_p`$ to equatorial radius $`R_e`$. We examine each sequence to determine whether the limits $`\beta _{shed}`$ and/or $`\beta _{bar}`$ are reached. Note that 3D hydrodynamical simulations are necessary to determine precisely $`\beta _{bar}`$ for each sequence. However, previous linear and nonlinear stability analyses indicate that $`0.14\beta _{bar}0.27`$ (Chandrasekhar, 1969; Shapiro and Teukolsky, 1983; Durisen and Tohline, 1985; Imamura, Friedman, and Durisen, 1985; Managan, 1985; Hachisu, Tohline, and Eriguchi, 1988; Tohline and Hachisu, 1990; Pickett, Durisen, and Davis, 1996; Toman et al., 1998; Imamura, Durisen, and Pickett, 2000). See §4 for further discussion of these stability limits.

We note that secular version of the bar instability, which sets in at a lower value of $`\beta `$, is not relevant for the scenario described here. The reason for this is that the viscosity that could drive the mode is assumed to be vanishingly small and, hence, the secular timescale is too long (Baumgarte and Shapiro, 1999). The timescale for gravitational radiation to drive the secular bar mode is also much longer than the evolution (cooling) timescale (Chandrasekhar, 1970).

In the low-viscosity and low-magnetic field limits, the angular momentum of mass shells on concentric cylinders in a contracting star will be preserved (Bodenheimer and Ostriker, 1973; Tassoul, 1978). Thus we confine our discussion to equilibrium sequences constructed with rotation laws that conserve the angular momentum distribution in this manner. The sequences examined here include one we construct below with the $`n^{}`$=3 rotation law (see §2.2). This sequence represents the cooling evolution of a SMS that is nearly spherical and rotating uniformly at formation prior to its contraction phase. It is the low-viscosity, low-field analog of the uniformly rotating sequence considered in Baumgarte & Shapiro (1999); in the present case, differential rotation will ensue once cooling gets underway. Because the initial rotation profiles of SMSs are unknown, we also examine alternative sequences constructed by Hachisu, Tohline, and Eriguchi (1988; hereafter HTE) that describe the evolution of SMSs with nonuniform initial rotation laws.

The remainder of this paper is organized as follows. In §2, we describe the numerical methods used to construct our $`n^{}`$=3 equilibrium sequence and the sequences constructed by HTE. We present our $`n^{}`$=3 sequence and summarize the sequences constructed by HTE in §3. In §4, we outline possible evolutionary scenarios for SMSs, based on the sequences described in §3. We discuss these results in §5.

## 2 Numerical Methods for SMS Model Construction

We have constructed an equilibrium sequence of SMS models, with a rotation law appropriate for SMSs that rotate uniformly prior to contraction. The individual models along the sequence have the same rotation law, $`M`$, and $`J`$, but decreasing entropy and axis ratios $`R_p/R_e`$. Each sequence is a quasistatic approximation to the evolution of a SMS that is contracting due to thermal emission. Modeling this phase of the evolution of SMSs with such equilibrium sequences is appropriate as the cooling timescale is longer than the hydrodynamic timescale (for $`M10^{13}M_{\mathrm{}}`$, Baumgarte and Shapiro 1999, and references therein)

### 2.1 The Self-Consistent Field Method

The structure of a fluid rotating about the $`z`$ axis with constant angular velocity $`\mathrm{\Omega }=\mathrm{\Omega }(\varpi )`$, where $`\varpi `$ is the distance from the rotation axis, is described by the following expression:

$$\frac{1}{\rho }P+\mathrm{\Phi }+h_0^2\mathrm{\Psi }(\varpi )=0,$$
(1)
where $`P`$ is the pressure, $`\mathrm{\Phi }`$ is the gravitational potential, $`\mathrm{\Psi }(\varpi )=(1/h_0^2)\mathrm{\Omega }^2(\varpi )\varpi 𝑑\varpi `$ is the centrifugal potential, and $`h_0`$ is a constant. Such a fluid is said to be in hydrostatic equilibrium because the forces due to its pressure and to its gravitational and centrifugal potentials are in balance.

The method we have used to construct SMSs in hydrostatic equilibrium is Hachisu’s grid based, iterative, axisymmetric Self-Consistent Field (HSCF) technique (Hachisu 1986, see also New 1996).

#### 2.1.1 Equation of state

The HSCF method requires the fluid to have a barotropic equation of state (EOS) $`P=P(\rho )`$. All of the SMS models we have constructed have a polytropic EOS for which

$$P=K\rho ^{1+1/n},$$
(2)
where $`n`$ is the polytropic index and $`K`$ is the entropy constant. The enthalpy of a polytrope is

$$H=(n+1)K\rho ^{1/n}.$$
(3)
The structure of a SMS is well represented by an $`n`$=3 polytrope (see §17.2 of Shapiro and Teukolsky 1983). This is because SMSs are dominated by thermal radiation pressure (for $`M10^6M_{\mathrm{}}`$, Zel’dovich and Novikov 1971; Fuller, Woosley, and Weaver 1986) and are convective with constant entropy per baryon throughout (see Loeb and Rasio 1994 for a simple proof of this property). The value of $`K`$ is determined by the entropy per baryon $`s`$ in the star and is given by

$$K=\frac{a}{3}\left(\frac{3s}{4m_Ha}\right)^{4/3},$$
(4)
where $`a`$ is the radiation density constant and $`m_H`$ is the mass of a hydrogen atom (see eq. 17.2.6 in Shapiro and Teukolsky 1983). The constant $`s`$ can be expressed in terms of the matter temperature $`T_m`$, baryon denisty $`n_b`$, and $`a`$:

$$s=\frac{4}{3}\frac{aT_m^3}{n_b}.$$
(5)
For polytropic EOSs, the sequences constructed with Hachisu’s SCF method are parameterized solely by the axis ratio $`R_p/R_e`$. The method requires that the grid cell locations of the two surface points, $`A=R_e`$ and $`B=R_p`$, be specified.

#### 2.1.2 Rotation Laws

In the low-viscosity limit, the angular momentum of a contracting star will be preserved on cylindrical mass shells (Bodenheimer and Ostriker, 1973; Tassoul, 1978). Thus the rotation laws used in sequences representative of such evolution must enforce this conservation. The $`n^{}`$=3 law and each of the laws employed in the sequences constructed by HTE are defined such that the specific angular momentum profile $`j(m)`$ is identical for each model constructed. Here $`m`$ is the mass interior to cylindrical radius $`\varpi `$. Because $`j`$ is the angular momentum per unit mass and is constant on cylindrical surfaces, a $`j(m)`$ law does ensure that the angular momentum distribution on cylindrical mass shells will be preserved from model to model along the sequence.

#### 2.1.3 Converged Configurations

The quality of the converged configurations can be estimated in terms of the virial error $`VE`$. The $`VE`$ provides a measure of how well the energy is balanced and is defined as

$$VE2T+W+3P𝑑V,$$
(6)
where V is the volume of the model. The quantity $`VE`$ is zero in strict equilibrium. The $`VE`$s for models on this sequence range from $`10^2`$ (for the very smallest axis ratios) to $`10^6`$, where the individual terms in equation (9) are of order unity (see Table 1).

Our axisymmetric equilibrium configurations were constructed on uniform, $`1024\times 1024`$ $`(\varpi ,z)`$ grids. Equatorial symmetry through the $`z`$=0 plane was assumed; this symmetry condition is valid for systems with $`\mathrm{\Omega }(\varpi )`$ and barotropic EOSs (Tassoul, 1978).

Note that initially we had difficulty obtaining converged models on the $`n`$=$`n^{}`$=3 sequence (see §2.2, 3.1) with axis ratios $`<0.400`$. For these models, the HSCF iterations oscillated between two or more states without converging. In order to obtain converged models with these axis ratios, we implemented a modification to the HSCF algorithm suggested in Pickett, Durisen, and Davis (1996). Instead of using $`\rho _i`$ from iteration $`i`$ as input for the next iteration step $`i+1`$, a combination of $`\rho _i`$ and the old density $`\rho _{i1}`$ is used (“under-relaxation”):

$$\rho _i^{^{}}=(1\delta )\rho _i+\delta \rho _{i1},$$
(7)
where $`\delta `$ is a constant (in this work, $`0.75\delta 0.99`$).

### 2.2 Scale of Models

As mentioned above, polytropic sequences constructed with the HSCF code are parameterized soley by the axis ratio $`R_p/R_e`$. That is, for a given EOS and rotation law, the choice of $`R_p/R_e`$ produces a unique model. A dimensional scale (i.e., mass, radius, etc.) for this model may be chosen after its construction. The individual models along a sequence representative of the cooling evolution of a SMS have the same $`M`$ and $`J`$. As will be demonstrated below, the choice of $`M`$ and $`J`$ for a sequence sets the dimensional scale for its individual models.

The HSCF computations are carried out in dimensionless form, with all quantities normalized in terms of $`G`$, $`R_e`$, and $`\rho _{max}`$. Here $`\rho _{max}`$ is the maximum density of the model. In HSCF units, $`M`$ and $`J`$ are normalized according to:

$`\widehat{M}`$ $`=`$ $`{\displaystyle \frac{M}{R_e^3\rho _{max}}},`$ (8a)
$`\widehat{J}`$ $`=`$ $`{\displaystyle \frac{J}{G^{\frac{1}{2}}R_e^5\rho _{max}^{\frac{3}{2}}}}.`$ (8b)
Here, and in what follows, carets denote normalized quantities.

The HSCF method produces unique values of the normalized quantites (e.g., $`\widehat{M}`$, $`\widehat{J}`$, etc.) for each input axis ratio $`R_p/R_e`$. After a converged model has been constructed, the normalized quantities are known. If physical values of $`M`$ and $`J`$ are then chosen to specify the sequence, equations (8) can be used to compute $`R_e`$ and $`\rho _{max}`$ for each model along the sequence:

$$R_e=G^1\left(\frac{J}{\widehat{J}}\right)^2\left(\frac{\widehat{M}}{M}\right)^3,$$
(9)
$$\rho _{max}=\frac{M}{\widehat{M}R_e^3}.$$
(10)
Once $`R_e`$ and $`\rho _{max}`$ are computed, they can be used to convert normalized quantities to physical units. For example, the definition of $`H`$ from equation (3) and its normalization, $`\widehat{H}=H/(GR_e^2\rho _{max})`$, can be used to determine a physical value for $`K`$:

$$K=\frac{1}{1+n}G\widehat{H}_{max}R_e^2\rho _{max}^{11/n}.$$
(11)
Note that when $`n`$=3, equation (10) leads to the following simplification of equation (11):

$$K=\frac{1}{4}G\widehat{H}_{max}\left(\frac{M}{\widehat{M}}\right)^{\frac{2}{3}}.$$
(12)
This is the well known result that the scale of $`K`$ depends solely on the choice of $`M`$, for $`n`$=3 polytropes.

## 3 Equilibrium Sequences

In this section, we summarize the properties of $`n`$=3 polytropic equilibrium sequences. The $`n^{}`$=3 sequence constructed by the present authors is discussed in §3.1 and the sequences constructed by HTE are discussed in §3.2.

### 3.1 $`n^{}`$=3

As mentioned in §1, the $`n^{}`$=3 law corresponds to the rotation profile of a uniformly rotating, spherical $`n`$=3 polytrope. In the low-viscosity limit, a sequence constructed with this law is representative of the evolution of a SMS that rotates uniformly prior to its cooling/contraction phase. Thus this sequence begins with the same starting model as that used in Baumgarte and Shapiro’s study of SMS evolution in the high-viscosity limit (Baumgarte and Shapiro, 1999). The $`n^{}`$=3 specific angular momentum distribution is given by (Bodenheimer and Ostriker, 1973):

$$j(m)=a_1+a_2\left(1\frac{m(\varpi )}{M}\right)^{\alpha _2}+a_3\left(1\frac{m(\varpi )}{M}\right)^{\alpha _3},$$
(13)
where $`M`$ is the total mass of the system, and the numerically determined constants are $`a_1=13.27`$, $`a_2=163.3`$, $`a_3=176.5`$, $`\alpha _2=0.2353`$, and $`\alpha _3=0.2222`$.

In Table 1, we list the properties of selected configurations constructed for the $`n^{}`$=3 rotation law. The data given (in HSCF normalized units, see §2.2) include the axis ratio $`R_p/R_e`$, polytropic entropy constant $`K`$ (normalized to the entropy constant of the initial spherical model $`K_0`$), mass $`\widehat{M}`$, angular momentum $`\widehat{J}`$, maximum enthalpy $`\widehat{H}_{max}`$, total energy $`\widehat{E}=E/GR_{e}^{}{}_{}{}^{5}\rho _{max}^{}{}_{}{}^{2}`$, ratio of kinetic to potential energy $`\beta `$, ratio of the centrifugal force $`F_c`$ to the gravitational force $`F_g`$ on the equator, and the virial error $`VE`$. When the quantity $`(F_c/F_g)_{eq}`$, equals unity, the mass shedding limit $`\beta _{shed}`$ has been reached.

Figure 1 displays $`R_p/R_e`$, $`\beta `$, $`(F_c/F_g)_{eq}`$, and $`E`$ as a function of $`K/K_0`$ for the individual models on this sequence. Note that $`K`$, which is a measure of the specific entropy, decreases monotonically with $`R_p/R_e`$ along the sequence. (The ratio $`K/K_0`$ can be related to time once the cooling law is specified; this has been treated by Baumgarte and Shapiro 1999 for the case of uniform rotation.) The energy $`E`$ also decreases monotonically. Density contour plots of selected sequence models are shown in Figure 2.

This sequence terminates due to mass shedding. The last spheroidal model we were able to construct had an axis ratio $`R_p/R_e=0.002`$ and $`(F_c/F_g)_{eq}`$=0.99. We were unable to construct a converged model with $`R_p/R_e`$=0.001, as the iteration process oscillated between models with $`(F_c/F_g)`$ less than and greater than one. Thus we conclude that this sequence terminates due to mass shedding at or very near $`R_p/R_e`$=0.001, with $`\beta _{shed}`$ 0.30. Note that there is no physical effect that forces the mass shedding limit to occur at an axis ratio of 0. In fact, these two points do not coincide on equilibrium sequences constructed with other equations of state and/or rotation laws (see, e.g, Hachisu 1986).

Note that no toroidal branch of this sequence exists. The HSCF code can construct tori if the inner radius of the torus $`R_i`$ is specified instead of $`R_p`$. All of the torioidal configurations we attempted to construct did not converge for this sequence (and had $`(F_c/F_g)_{eq}>1`$ during the iteration process).

Recall that the dynamical bar mode instability is likely to set in at $`\beta _{bar}0.27`$. The model with $`\beta \beta _{bar}`$=$`0.27`$ on the $`n^{}`$=$`3`$ sequence has $`R_p/R_e`$=$`0.004`$ and $`K/K_0`$=$`0.381`$.

In order to estimate the possible physical scales of SMS models on this sequence, we must choose the constant values $`M`$ and $`J`$ (see §2.2). For example, say the initial nearly spherical star has $`M`$=$`10^6M_{\mathrm{}}`$, $`(R_e)_0`$=$`10^{17}\mathrm{cm}`$, and $`\beta _0=(T/|W|)_0`$=$`10^5`$, where a subscript 0 denotes the value is for the (nearly) spherical star. The value of $`(R_e)_0`$ is chosen such that the compaction, $`(G/c^2)M/R`$, for the model that reaches $`\beta \beta _{bar}`$=$`0.27`$ is of the same order of magnitude as the uniformly rotating critical configuration of Baumgarte and Shapiro (1999) for gravitational collapse due to relativistic gravitation. For all larger values of the product $`\beta _0(R_e)_0`$ , an evolved configuration with $`M`$=$`10^6M_{\mathrm{}}`$ will experience a dynamical bar instability prior to reaching the critical compaction for radial collapse. The value of $`\beta _0`$ is quite arbitrary except for the slow rotation constraint, $`\beta _01`$, necessary to ensure that the initial cloud is spherical. Cosmological considerations and calculations of early structure formation will ultimately determine the likely distribution of parameters for primordial SMSs.

Having chosen $`M`$, $`(R_e)_0`$ and $`\beta _0`$, one can exhibit the scaling of all of the remaining quantities derived below in terms of these choices, i.e. in terms of ($`M_6M/10^6M_{\mathrm{}}`$), ($`R_{17}(R_e)_0/10^{17}\mathrm{cm}`$), and ($`\beta _5\beta _0/10^5`$).

The choice of $`M`$, $`(R_e)_0`$ and $`\beta _0`$ determines $`J`$. To see this, recall that for a uniformly rotating star $`J=(2TI)^{1/2}`$, where $`I`$ is the moment of inertia. On this sequence, $`|\widehat{W}|_0`$=$`8.96\times 10^3`$ and $`\widehat{I}_0`$=$`5.82\times 10^3`$. Thus, $`\widehat{J}_0`$= $`3.23\times 10^5\beta _5^{1/2}`$. Equations (8b) and (10) can then be used to convert $`\widehat{J}_0`$ to $`J`$=$`1.09\times 10^{61}[\beta _5R_{17}M_6^3]^{1/2}\mathrm{g}\mathrm{cm}^2\mathrm{s}^1`$. These choices of $`M`$ and $`J`$ for the sequence can be used to dimensionalize its individual models according to equations (9-12). For instance, the nearly spherical initial model has $`(\rho _{max})_0`$= $`2.57\times 10^{11}M_6R_{17}^3\mathrm{g}\mathrm{cm}^3`$ and $`K_0`$=$`3.84\times 10^{18}M_6^{2/3}\mathrm{g}^{1/3}\mathrm{cm}^3\mathrm{s}^2`$. The model with $`\beta \beta _{bar}=0.27`$ and axis ratio 0.004 has $`R_e`$=$`4.14\times 10^{13}\beta _5R_{17}\mathrm{cm}`$, $`\rho _{max}`$=$`1.27\times 10^5M_6[\beta _5R_{17}]^3\mathrm{g}\mathrm{cm}^3`$, and $`K`$=$`1.46\times 10^{18}M_6^{2/3}\mathrm{g}^{1/3}\mathrm{cm}^3\mathrm{s}^2`$. The model with axis ratio 0.002 has $`R_e`$=$`4.04\times 10^{13}\beta _5R_{17}\mathrm{cm}`$, $`\rho _{max}`$=$`4.18\times 10^5M_6[\beta _5R_{17}]^3\mathrm{g}\mathrm{cm}^3`$, and $`K`$=$`1.20\times 10^{18}M_6^{2/3}\mathrm{g}^{1/3}\mathrm{cm}^3\mathrm{s}^2`$.

### 3.2 Hachisu, Tohline, and Eriguchi 1988

Hachisu, Tohline, and Eriguchi (1988), hereafter HTE, used the SCF method to construct $`n`$=3 polytropic sequences with the following parameterized angular momentum distribution:

$$j(m)=(1+q)(J/M)\left[1(1m(\varpi )/M)^{1/q}\right].$$
(14)
Here, the index $`q`$ specifies the rotation law. Note that the limiting case of $`q`$=0 corresponds to the j-constant rotation law.

The j-constant ($`q`$=0) sequence HTE constructed is entirely toroidal. That is, no spheroidal models exist with this rotation law (for any of the polytropic indices $`0n3`$ they investigated). Thus the j-constant sequence is not relevant to the discussion of the cooling phase of SMS evolution, which begins with a spheroidal configuration.

HTE also constructed sequences with $`q`$=0.5, 0.7, 1.0, 1.3, and 1.5. All of these sequences make continuous transitions between spheroidal and toroidal configurations at relatively high values of $`\beta _{trans}>`$ 0.33, without encountering mass shedding. Note that if $`M`$ and $`J`$ are held constant, the entropy $`K`$ and energy $`E`$ decrease monotonically along the spheriodal branch of each of these sequences. This confirms that they are suitable representations of the cooling evolution of SMSs.

From HTE’s Figure 7f, one can estimate that $`K_{bar}/K_00.4`$ for each of these sequences (where $`K_{bar}`$ is the value of K for a model with $`\beta `$=0.27). Thus the onset of the bar mode instability for HTE’s sequences occurs at approximately the same value of $`K/K_0`$ as on the $`n^{}`$=$`3`$ sequence.

## 4 Evolutionary Scenarios

In the low-viscosity limit, the outcome of the cooling evolution and resultant spin up of a SMS depends critically on its initial angular momentum distribution. Because the nature of the angular momentum distribution in SMSs is largely unknown, we have presented equilibrium sequences of SMS models with several different rotation laws in §3. Here we outline possible scenarios for the outcome of the cooling evolution of SMSs with low viscosity, based on these sequences.

As mentioned in §1, one possible outcome is that a SMS will spin up to a critical value $`\beta _{bar}`$, at which it will encounter the dynamical bar mode instability. As the bar mode develops, angular momentum and mass will be transported outwards. The angular momentum removed will hasten, to some degree, the eventual collapse of the SMS (and the subsequent possible formation of a SMBH). A SMS undergoing the bar mode is potentially a strong source of gravitational radiation because of its highly nonaxisymmetric structure and hence might be detectable by LISA.

The nonlinear nature of the bar mode instability necessitates the use of 3D hydrodynamical simulations to precisely determine $`\beta _{bar}`$ for each equilibrium sequence. However, previous linear and nonlinear stability analyses indicate that $`0.14\beta _{bar}0.27`$ for a wide range of polytropic indices and rotation laws (Shapiro and Teukolsky, 1983; Durisen and Tohline, 1985; Imamura, Friedman, and Durisen, 1985; Managan, 1985; Hachisu, Tohline, and Eriguchi, 1988; Tohline and Hachisu, 1990; Pickett, Durisen, and Davis, 1996; Toman et al., 1998; Imamura, Durisen, and Pickett, 2000). The lower limit, $`0.14`$, may be most appropriate for spheroids with off-center density maxima and tori (Hachisu, Tohline, and Eriguchi, 1988; Tohline and Hachisu, 1990).

The nonlinear hydrodynamics study of Pickett, Durisen, and Davis (1996) indicates that $`\beta _{bar}`$ may be less than $`0.27`$ for rotation laws that place a significant amount of angular momentum in equatorial mass elements. Their results predict that the $`m`$=2 stability limit may be less than $`0.20`$ for models with the $`n^{}`$=3 rotation law. Their simulations, of $`n`$=1.5 polytropes, also suggest that a one-armed spiral, $`m`$=1 mode may become increasingly dominant over the $`m`$=2 mode as the equatorial concentration of angular momentum increases. Note that the grid resolutions used in Pickett, Durisen, and Davis (1996) were likely not sufficient to accurately model the development of instabilities in models with these extreme differential rotation laws (Toman et al., 1998). However, the results of Pickett, Durisen, and Davis (1996) and the linear and nonlinear stability analyses of Toman et al. (1998) confirm that $`0.27`$ is likely to be an approximate upper limit to $`\beta _{bar}`$, for a variety of polytropic indices and rotation laws. The analysis presented in the remainder of this manuscript assumes that the $`m`$=2 bar mode is the dominant mode and that $`\beta _{bar}0.27`$.

The new $`n^{}`$=3 sequence constructed by the present authors and the sequences constructed by HTE that are considered in this work indicate that a bar mode phase is probable for SMSs with a wide range of differential rotation laws. Indeed, HTE’s $`q`$=0.1-1.5 (see equation 14) sequences do not terminate due to mass shedding. Each of these sequences makes a continuous transition from spheroidal to toroidal configurations at values of $`\beta _{trans}>0.33>0.27\beta _{bar}`$. Thus, $`n`$=3 models with these $`q`$-indexed laws would likely encounter the bar mode as spheroids.

The spheroidal $`n^{}`$=3 sequence terminates due to mass shedding at $`\beta _{shed}0.30`$. Because $`\beta _{shed}>0.27`$, we expect that models with this rotation law become dynamically unstable to the bar mode near $`R_p/R_e0.004`$ (this model is marked with a solid dot in Figure 1), before reaching the mass shedding limit.

Note that even if the actual value of $`\beta _{bar}`$ is less than $`0.27`$, the qualitative nature of our results would not change. That is, the sequences examined here would still have models that are unstable to the bar mode. In addition, the quantitative estimates of the characteristics of the gravitational radiation emission presented in §5 would only change by a numerical factor of order 1-10, even if $`\beta _{bar}`$ were as low as $`0.22`$.

3D hydrodynamical simulations are needed to properly study the evolution of rapidly rotating stars that encounter the bar mode instability. To date, we are unaware of any published hydrodynamical studies of the bar mode in $`n`$=3.0 polytropes. However, many authors have simulated this instability in polytropes with lower $`n`$ and various rotation laws (e.g., Tohline, Durisen, and McCollough 1985; Pickett, Durisen, and Davis 1996; Imamura, Durisen, and Pickett 2000; New, Centrella, and Tohline 2000; Brown 2000). All of these simulations agree on the initial phase of the instability, during which the initially axisymmetric object deforms into a bar shape through which spiral arms are shed. Recent Newtonian simulations indicate that the long-term outcome of the bar instability in an $`n`$=1.5 polytrope (with an initial Maclaurin rotation law) is a persistent bar-like structure that emits gravitational radiation over many rotation periods (New, Centrella, and Tohline, 2000; Brown, 2000) (the persistent nature of the bar structure was recently confirmed with the post-Newtonian simulations of Saijo et al. 2000). Simulations relevant to the dynamical bar mode instability have recently been performed for $`n`$=1 polytropes in full 3D general relativity by Shibata, Baumgarte, and Shapiro (2000).

## 5 Conclusions and Future Work

The thermal emission of a rotating supermassive star causes it to contract and spin up. If the viscosity and internal magnetic fields are weak, the SMS will rotate differentially during its cooling and contraction phase. We have examined models of the quasistatic cooling evolution of SMSs, represented by equilibrium sequences with six differential rotation laws.

The evolution of a cooling SMS that is initially nearly spherical, with very low, uniform rotation can be represented by the $`n^{}`$=3 sequence we constructed. This sequence terminates due to mass–shedding at $`\beta _{shed}0.30`$. This $`\beta _{shed}`$ exceeds the value $`0.27`$, which is a likely upper limit for the onset of the rotationally induced, dynamical bar instability. Thus it appears likely that a star with this rotation law would become unstable to the bar mode prior to reaching the mass shedding limit.

We have also considered sequences constructed by HTE, in order to follow the evolution of SMSs with five different initial differential rotation laws. No mass–shedding limits exist on these sequences. Thus, a cooling SMS characterized by one of these rotation profiles will also encounter the dynamical bar mode.

This mode will transport mass and angular momentum outward and thus may hasten the onset of collapse and possible SMBH formation (Tohline, Durisen, and McCollough, 1985; Pickett, Durisen, and Davis, 1996; Imamura, Durisen, and Pickett, 2000; New, Centrella, and Tohline, 2000; Brown, 2000). Because of their nonaxisymmetric structure, SMSs undergoing the bar mode could be strong sources of quasi-periodic gravitational radiation for detectors like LISA, even before collapse.

The frequency of these quasiperiodic gravitational waves can be estimated from the expected bar rotation rate $`\mathrm{\Omega }_{bar}`$. The model on our $`n^{}`$=$`3`$ sequence with $`\beta \beta _{bar}`$=$`0.27`$ has a central rotation rate $`\mathrm{\Omega }_c`$=$`1.02\times 10^1\mathrm{Hz}[M_6\beta _5^3R_{17}^3]^{1/2}`$. In previous hydrodynamics simulations of the bar mode instability (New, Centrella, and Tohline, 2000), $`\mathrm{\Omega }_{bar}`$ was $`0.4\mathrm{\Omega }_c`$. With this relation between $`\mathrm{\Omega }_{bar}`$ and $`\mathrm{\Omega }_c`$, the gravitational wave frequency $`f_{GW}`$ can be estimated to be

$`f_{GW}`$ $`=`$ $`2f_{bar}=2{\displaystyle \frac{\mathrm{\Omega }_{bar}}{2\pi }}{\displaystyle \frac{0.4}{\pi }}\mathrm{\Omega }_c`$ (15)
$``$ $`1\times 10^2\mathrm{Hz}[M_6\beta _5^3R_{17}^3]^{\frac{1}{2}}.`$
For a SMS of $`10^6M_{\mathrm{}}`$, which begins as a slowly rotating star of radius $`10^{17}\mathrm{cm}`$ with $`\beta _0`$=$`10^5`$, this yields a frequency of $`1\times 10^2\mathrm{Hz}`$. This frequency is in the range in which LISA is expected to be most sensitive, $`10^4`$-$`1\mathrm{H}\mathrm{z}`$ (see, e.g., Thorne 1995, Folkner 1998). The choice $`\beta _0`$=$`2\times 10^4/R_{17}`$ is the maximum value of $`\beta _0`$ for which $`f_{GW}`$ (=$`10^4\mathrm{Hz}`$) is still in LISA’s range of sensitivity.<sup>1</sup><sup>1</sup>1We note that Peebles (1969) has discussed the spin-up of primordial masses by tidal torquing and finds that $`JMR^2`$. Naively scaling his results (derived for $`M10^{11}M_{\mathrm{}}`$ configurations) to our regime, we find that $`\beta _010^4`$, so that the resulting frequency of gravitational radiation would lie near the lower limit of LISA’s range.

The strength of the gravitational wave signal can be estimated roughly to be

$`h`$ $``$ $`{\displaystyle \frac{G}{c^4}}{\displaystyle \frac{\ddot{Q}}{d}}{\displaystyle \frac{G}{c^4}}{\displaystyle \frac{MR_{bar}^2f_{bar}^2}{d}}`$ (16)
$``$ $`4\times 10^{15}\left({\displaystyle \frac{d}{1\mathrm{Gpc}}}\right)^1M_6^2\beta _5^1R_{17}^1,`$
where $`\ddot{Q}`$ is the second time derivative of the star’s quadrupole moment and $`d`$ is the distance, which we scale to 1 Gpc (the Hubble distance is $``$ 3 Gpc). Here we have used $`R_{bar}R_e`$=$`4.14\times 10^{13}\beta _5R_{17}\mathrm{cm}`$ for the $`n^{}`$=$`3`$ model that reaches the point $`\beta \beta _{bar}`$.

As mentioned previously, recent simulations of the bar instability indicate that the outcome is a long-lived bar-like structure (New, Centrella, and Tohline, 2000; Brown, 2000; Saijo et al., 2000). The bar will decay on a secular timescale due to dissipative effects. For differentially rotating SMSs, the largest source of dissipation will be gravitational radiation. The gravitational radiation damping timescale $`\tau _{GW}`$ is approximately

$$\tau _{GW}\frac{T}{\left(\frac{dE}{dt}\right)_{GW}},$$
(17)
where $`T`$ is the rotational kinetic energy and $`(dE/dt)_{GW}`$ is the rate at which gravitational radiation carries energy away from the system. Recall that $`T=\beta |W|`$. The gravitational potential energy $`|W|GM^2/R_{bar}`$. Thus,

$$T=\beta \frac{GM^2}{R}.$$
(18)
The radiation rate can be estimated as (Shapiro and Teukolsky, 1983)

$$\left(\frac{dE}{dt}\right)_{GW}\frac{G}{c^5}\frac{M}{R^2}v^6.$$
(19)
In this case the characteristic velocity of the system is $`v=(\beta GM/R)^{1/2}.`$ Substitution of equations (18) and (19) into equation (17) yields

$`\tau _{GW}`$ $``$ $`{\displaystyle \frac{c^5}{G^3}}{\displaystyle \frac{R_{bar}^4}{\beta ^2M^3}}`$ (20)
$``$ $`1\times 10^4yrsM_6^3[\beta _5R_{17}]^4.`$
The number of cycles $`𝒩`$ for which the signal will persist is

$`𝒩`$ $`=`$ $`\tau _{GW}f_{GW}`$ (21)
$`=`$ $`4\times 10^9[M_6\beta _5R_{17}]^{5/2}.`$
The quasiperiodicity of such a signal will assist in its detection (Schutz, 1997).

The fraction of the mass $`(dM/M)_{GW}`$ radiated via gravitational radiation over the interval $`\tau _{GW}`$ can be estimated as

$`\left({\displaystyle \frac{dM}{M}}\right)_{GW}`$ $``$ $`\tau _{GW}\left({\displaystyle \frac{dE}{dt}}\right)_{GW}{\displaystyle \frac{1}{Mc^2}}`$ (22)
$``$ $`{\displaystyle \frac{G}{c^2}}{\displaystyle \frac{\beta M}{R}}`$
$``$ $`1\times 10^3M_6[\beta _5R_{17}]^1.`$
It appears that SMSs with a wide range of differential rotation laws may become unstable to the development of the dynamical bar mode. However, linear and nonlinear stability analyses are needed to precisely determine the onset of the $`m`$=2 dynamical instability for the equilibrium sequences considered here. Such analyses are also needed to determine the relative importance of various unstable nonaxisymmetric modes in SMS models (as a dominant $`m`$=1 mode would change the characteristics of the gravitational radiation emission). In addition, three-dimensional hydrodynamical simulations are necessary to study the mass and angular momentum redistribution induced by these nonaxisymmetric instabilities, to compute the gravitational waveforms emitted, and to determine the final fate of the star.

Future investigations involving hydrodynamical simulations of SMS models would benefit from improved knowledge of initial conditions, such as an appropriate value for $`\beta _0`$. These appropriate initial conditions could be determined from studies of large-scale structure and cosmology.

Most importantly, studies are required to assess the roles of viscosity and magnetic fields in rotating SMSs to judge whether they are sufficient to drive these configurations to uniform rotation prior to bar instability.

We thank Joel Tohline for helpful conversations. This work has been supported in part by NSF Grants AST 96-18524 and PHY 99-02833 and NASA Grants NAG5-7152 and NAG5-8418 to the University of Illinois at Urbana-Champaign. A portion of this work was performed under auspices of the U.S. Department of Energy by Los Alamos National Laboratory under contract W-7405-ENG-36.


text:
# Magnetic and electronic structures of superconducting RuSr2GdCu2O8

## I Introduction

The reported discovery of the coexistence of ferromagnetism and superconductivity by Tallon et. al.  has attracted a great deal of interest, since ferromagnetism breaks the degeneracy of spin-up and spin-down partners of Cooper pairs. The experiments demonstrated that RuSr2GdCu2O8 exhibits ferromagnetic ordering of Ru moments below a Curie temperature, $`T_C`$ 133 K, and becomes superconducting at a lower temperature, $`T_c`$ up to 40 K. From a theoretical point of view, Pickett et. al. calculated the electronic structure of RuSr2GdCu2O8 and discussed how the superconductivity, which may behave as a Fulde-Ferrell-Larkin-Ovchinniko (FFLO) type superconductor, can coexist with the ferromagnetism. However, the X-ray determined crystal structure used in their calculations was different from that recently determined by neutron diffraction \- and so they did not consider the structural distortions (i.e., oxygen rotations) in the RuO6 octahedra.

The crystal structure of RuSr2GdCu2O8 with a P4/mbm tetragonal space group, determined by the neutron experiment, is shown in Fig. 1 (a). This is a structure similar to that of YBa2Cu3O7, where Y, Ba and Cu (chain atom) are replaced by Gd, Sr and Ru, respectively. Ru lies at a six-coordinated position in the octahedron composed of six neighboring oxygens (four ORu and two Oapical), while Cu lies in a five-coordinated position (four OCu and one Oapical). Since the interatomic distance of Cu-Oapical (2.190 Å) is much larger than that of Ru-Oapical (1.912 Å), the Cu has a rather weak interaction with Oapical. Thus, this may yield a two-dimensional (2D) electronic structure of the Cu-OCu layer. The RuO6 octahedra are rotated by about 14 around the c-axis to fill the space more efficiently, as shown in Fig. 1 (b). These rotations may lead to a significant reduction in the $`d`$-electron bandwidth due to a deviation of the angle ∠Ru-ORu-Ru from 180, and so the magnetism of Ru may be strongly perturbed by the structural distortion. Hence, it is necessary to again investigate the stability of the magnetic structure in RuSr2GdCu2O8 by first-principles calculations with the precise crystal structure parameters.

With the experimentally determined structure, our highly precise first-principles full-potential linearized augmented plane wave(FLAPW) calculations on three different magnetic orderings demonstrate that antiferromagnetic ordering of the Ru moments is energetically favorable over the previously proposed ferromagnetic ordering. Furthermore, in agreement with our result, recent neutron diffraction experiments by Lynn et. al. observed antiferromagnetic Ru ordering. In the present paper, we report the first-principles results on RuSr2GdCu2O8, and discuss the magnetic and electronic structures and their implications for the coexistence of superconductivity.

## II Model and calculation method

In our FLAPW calculations, we employed the experimentally determined room temperature crystal structure. For the magnetic structures, we assumed three kinds of magnetic orderings for Ru and Gd moments as depicted in Fig.2: (a) a FM structure with ferromagnetic orderings for both Ru and Gd moments; (b) an AFM-I structure with a C-type antiferromagnetic ordering of Ru moments and a ferromagnetic ordering of Gd moments (more precisely, a ferrimagnetic structure); and (c) an AFM-II structure with C-type antiferromagnetic orderings for both Ru and Gd moments. (Note that the symmetry of both the AFM-I and AFM-II structures makes two inequivalent sublattice sites for Ru, Oapical and Cu.) While the recent experiment revealed antiferromagnetic alignments for both Ru and Gd moments with G-type orderings, in which nearest neighbors of Ru (Gd) moments along all three crystallographic axes are coupled antiferromagnetically, we employed the C-type antiferromagnetic ordering in order to reduce the large computational effort this would entail. (The G-type ordering requires doubling the unit cell of the C-type ordering.) This is justified because the moment alignment along the c-axis may be less important than that along the a-axis since the distance between the neighboring Ru (Gd) atoms in the c-axis (11.56 Å) is significantly greater than that in the a-axis (3.84 Å). Thus, the assumption of C-type magnetic structure would not alter the results or the insight obtained.

Calculations were performed based on the local spin density approximation (LSDA) with the Hedin-Lundquist exchange correlation, in which the core states are treated fully relativistically and the valence states are treated semi-relativistically. Although the effects of electronic correlation in a strongly correlated system may be taken into account within a scheme such as LDA+U, RuSr2GdCu2O8 shows metallic character even in the RuO layer, as will be presented in Sec. IV, which causes that effect to be weak. Hence, while our results are restricted within LDA, they may be sufficient to discuss qualitative conclusions. The LAPW basis functions were used with a cut-off, $`|k+G|<3.6`$ a. u., corresponding to about 1800 plane waves. Muffin-tin (MT) sphere radii in a. u. were chosen as: 2.4 (Ru), 2.5 (Sr), 2.5 (Gd), 2.2 (Cu) and 1.2 (O). Inside the spheres, the angular momentum expansion was truncated at $`l=8`$ for wave functions, charge density and potential. Integrations were performed over a first Brilliouin zone using 126 special $`k`$-points, corresponding to 16, 16 and 32 $`k`$-points in the irreducible Brilliouin zone for the FM, AFM-I and AFM-II structures, respectively.

## III Results

### A Magnetic Moments

The calculated magnetic moment of each atom inside the MT sphere for the FM, AFM-I and AFM-II structures is given in Table I. The moment of Gd is close to 7 $`\mu _B`$ \- as expected. In the AFM-I structure, the magnitude of the Ru, Oapical and Cu magnetic moments with opposite spin directions is slightly different since the ferromagnetically ordered Gd moments break the symmetry of the two sublattice sites. Among all three structures, the Ru magnetic moments are almost of the same magnitude, about 1.55 $`\mu _B`$. Note that since the magnitude depends on the MT sphere radius, the rather large MT sphere radius (2.4 a.u.) we used gives rather large values. However, the experimental value (1.18 $`\mu _B`$) of the Ru magnetic moment determined by neutron diffraction is smaller than our results. This may be attributed to possible canted antiferromagnetic ordering of the Ru magnetic moments, for which further investigations are necessary. The induced magnetic moments of Oapical and ORu are very sensitive to the Ru magnetic structures, while the Cu and OCu moments are negligibly small.

### B Total energy

Total energies for all three FM, AFM-I and AFM-II structures are also given in Table I. From a comparison between the FM and AFM-I structures, we found that an antiferomagnetic ordering of Ru moments is favored over the ferromagnetic ordering. The total energy difference, 22.5 meV/Ru-atom, between both FM and AFM-I corresponds to the Neél temperature ($`T_N`$) of 73 K, crudely estimated by mean field theory with $`S=3/2`$. Although the calculated $`T_N`$ is lower than the experimentally determined $`T_N`$, 136 K, our prediction for the magnetic structure is consistent with the experiment by Lynn et. al. Further, we found that an antiferromagnetic ordering of the Gd moments (AFM-II) reduces the total energy, but only by 2.3 meV/Gd-atom, from that of the AFM-I structure, corresponding $`T_N7`$ K (assumed $`S=7/2`$), which is approximately in agreement with the observed low Neél temperature, 2.5 K, for the Gd moment alignment.

## IV Discussion

### A Density of states

Figure 3 shows the total density of states (DOS) and partial DOS of each atom for the FM, AFM-I and AFM-II structures. The peaks at 2 eV above the Fermi level ($`E_F`$) and 3 eV below $`E_F`$ in the total DOS for all three structures correspond to Gd 4$`f`$ states which are strongly localized, and their states do not significantly affect the electronic structure around $`E_F`$. An exchange splitting in the total DOS in the AFM-I structure is due to the ferromagnetism of Gd, but the difference between majority and minority spin DOS is very small except in the region of the localized Gd 4$`f`$ states. Note that the total DOS for all three structures shows metallic character.

From the partial DOS in Fig. 3, we can see that the magnetism of Ru is dominated by antibonding $`t_{2g}`$ states, and the major ORu and Oapical peaks correlate strongly with those of Ru. In contrast, the partial DOS for Cu $`d`$ and OCu $`p`$ states, which are strongly hybridized, are insensitive to the magnetism of Ru. The exchange splitting of the Ru $`d`$ states for all three structures is about 1 eV, while the induced exchange splitting in the Cu $`d`$ states is smaller by two orders of magnitude - 0.01 to zero eV. Our results for the exchange splitting and magnetic moments of the Cu $`d`$ states is smaller, even in the FM case, than those in the undistorted FM case reported by Weht et. al. This may be due to the differernt structure parameters employed in the two calculations. The fact that the Ru-Cu distance, 4.10 Å, in our calculations is larger than in their case, 3.58 Å, results in a rather weak magnetic interaction and leads to smaller exchange splittings and magnetic moments in the Cu $`d`$ states. Futher, the relative atomic position of Oapical between the Ru and Cu may be sensitive to the magnetism. However, the exchange splitting and the magnetic moments in the Cu $`d`$ are quite small. The small exchange coupling is due to the unique electronic structure of the layered Ru $`t_{2g}`$ and Cu $`e_g`$ states separated by $`p`$ orbitals of Oapical, as discussed by Pickett et. al., and is roughly valid regardless of the magnetic ordering of Ru moments: Ru $`t_{2g}`$ states couple only to the $`p_x`$ and $`p_y`$ orbitals of Oapical, but do not couple to the Cu-O $`d_{x^2y^2}`$ states.

### B Magnetism of Ru

We focus here on the magnetism of Ru. Figure 4 shows the projected DOS in real space for Ru $`d`$, ORu $`p`$ and Oapical $`p`$ states for the FM, AFM-I and AFM-II structures. The $`x`$ and $`z`$-axes for the Ru (ORu) DOS are chosen as directions to the neighboring ORu (Ru) site and c-axis, respectively. Compared to the DOS of the FM structure presented by Weht et. al., which did not include a structural distortion due to RuO6 rotations, our calculated majority $`d_{xy}`$ states (cf., Fig. 4 (a)) are more localized below $`E_F`$, due to deviation of the angle ∠Ru-ORu-Ru from 180 and the elongation of Ru-ORu bonds by the RuO6 rotations. This results in the $`d_{xy}`$ states being almost fully occupied, and in a depletion of the majority spin hybridization channel at $`E_F`$. In the AFM-I and AFM-II structures, the majority $`d_{xy}`$ and $`d_{xz(yz)}`$ states on a Ru site can hybridize with the minority spin states on neighboring Ru sites through ORu $`p_y`$ and ORu $`p_z`$ orbitals, respectively, by the superexchange mechanism. The $`d_{xz(yz)}`$ states also couple strongly to the Oapical $`p_x(y)`$ orbitals. The charge configuration of Ru seems to be $`t_{2g}^3`$ with a high spin state, namely close to Ru5+. Generally, the superexchange interaction between fully occupied $`t_{2g}^3`$ states in nearest neighbor magnetic ions tends to be stabilized with an antiferromagneticic ordering of their moments. On the other hand, the double exchange mechanism through itinerant electrons may not be favorable, since the majority $`t_{2g}`$ states are strongly localized and almost fully occupied - as seen in Fig. 4 (b) and (c). Therefore, the magnetism of Ru in RuSr2GdCu2O8 is determined by the superexchange mechanism.

### C Band structure and Fermi surface

The calculated band structures for RuSr2GdCu2O8 are shown along high symmetry directions in Fig. 5. The majority and minority spin states in AFM-II are degenerate due to the symmetry of the antiferromagnetic Ru and Gd moment alignments, as shown in Fig. 2. In the AFM-I structure, a small difference between majority and minority spin states was observed due to the assumed ferromagnetic Gd moments, but its difference is negligibly small. The band structure of AFM-I is practically the same as that observed in the AFM-II structure. The bands around $`E_F`$ for all three structures arise from the Ru-ORu-Oapical bands (solid light circles) and Cu-OCu bands (solid dark circles). In the AFM-I and AFM-II structures, there are three bands crossing $`E_F`$ \- as indicated in Fig. 5 (c). They arise from a Ru-ORu band with antibonding $`d_{xy}p_y`$ ($`dp\pi `$) orbitals, and two Cu-OCu bands composed of antibonding $`d_{x^2y^2}p_x`$ ($`dp\sigma `$) obitals with even and odd symmetry derived by a mirror operation in the Cu-OCu bilayer. The two Cu-OCu bands are almost the same as those in the FM structure. As expected from the structural similarity between RuSr2GdCu2O8 and YBa2Cu3O7, their band structures arising from Cu-OCu layers are close to each other, showing a strong two-dimensionality. The energy states are little changed upon going from $`\mathrm{\Gamma }`$ to Z along the $`z`$ direction. Also, the strong $`dp\sigma `$ hybridization leads to a wide bandwidth of 10.1, 9.7 and 9.7 eV for FM, AFM-I and AFM-II, respectively. These values are almost the same as observed in the 2D $`dp\sigma `$ bands in YBa2Cu3O7, 9.0 eV. The slightly larger bandwidth may be due to the smaller interatomic distance (1.93 Å) of Cu-OCu bonds and the smaller buckling angle (5.7) in the CuO2 plane of RuSr2GdCu2O8, compared with those of YBa2Cu3O7 (1.95 Å and 7.8).

From a comparison of band structures between the FM and AFM-I structures, we found that the antibonding $`dp\sigma `$ bands are not influenced significantly by the magnetism of Ru except around the M point in the zone. Around M, where they possess small components of Cu $`d_{3z^2r^2}`$ character, the Cu $`e_g`$ states couple to the Ru $`e_g`$ states through the Oapical $`p_z`$ orbitals, but do not couple to Ru $`t_{2g}`$ states. Therefore, the quite small exchange splitting in the $`dp\sigma `$ states may be mediated through a more indirect coupling path with Oapical-ORu. Note that the Cu $`t_{2g}`$ states are coupled to the Ru $`t_{2g}`$ states through the Oapical $`p_{x(y)}`$ orbitals, but the Cu $`t_{2g}`$ states are almost fully occupied, so the hybridization channel at $`E_F`$ is depleted.

The calculated Fermi surfaces (FS) at $`k_z=0`$ for the AFM-II structure are shown in Fig. 6 (a). For convenience, a sketch of the FS is given in Fig. 6 (b), which makes for an easier comparison with the high $`T_c`$ results. Two closed FS centered at the zone corner arising from the Cu-OCu bands and a closed FS from the Ru-ORu band, are found. The inner and outer FS (1 and 2 in Fig. 6 (b)) correspond to antibonding $`dp\sigma `$ states with even and odd symmetry under $`z`$ reflection, respectively. Note that the FS has no exchange splitting, due to the antiferromagnetic ordering of the Ru and Gd moments, implies that the superconductivity need no longer be of the FFLO type \- in contrast to the FM case. A folding in $`k`$-space due to the antiferromagnetic ordering of Ru makes intersections between the bands (a, b and c in Fig. 6 (b)), showing the band-band interaction. As observed in the undistorted FM, the FS clearly demonstrates similar nesting structures with those in the other high $`T_c`$ cuprate superconductors; these will give rise to singularities in the generalized susceptibility, and lead to possible anomalous behavior of the electronic properties. Further, the FS (3 in Fig. 6 (b)), derived from Ru $`d_{xy}`$ states, shows a pocket of electrons centered at $`\mathrm{\Gamma }`$, which may lead to the superconductivity in the Cu-OCu layer with hole-doped character, as observed in experiments

In contrast, the FS in the FM structure is complicated due to the Ru-ORu bands. The majority spin Ru $`d_{xz(yz)}`$ bands (4 in Fig. 5 (a)) become more dispersive and cross at $`E_F`$ compared to the AFM-I structure, which makes an additional hole pocket centered at $`\mathrm{\Gamma }`$, leading to be a less hole-doped character in Cu-OCu layers. Hence, the FM structure may show low $`T_c`$ or non-superconductivity.

## V summary

We have performed first-principles FLAPW calculations to investigate magnetic structures in superconducting RuSr2GdCu2O8. Contrary to the previously ferromagnetic RuSr2GdCu2O8 structure proposed from susceptibility experiments and first-principles calculations, we found that antiferromagnetic ordering of the Ru moments is energetically favored over ferromagnetic ordering. Our results are consistent with those of a recent neutron diffraction experiment, in which the antiferromagnetic ordering of Ru moments was predominantly observed. The magnetism arises from Ru $`t_{2g}`$ states through the superexchange mechanism. The Ru $`t_{2g}`$ states have a very small interaction with the Cu $`e_g`$ states \- which results in quite a small exchange splitting in the antibonding $`dp\sigma `$ states. The Fermi surfaces derived from the antibonding $`dp\sigma `$ states have similar nesting structures with those in the 2D high $`T_c`$ cuprate superconductors, which may give a basis for superconductivity in this material. Further experimental and theoretical investigations are necessary to confirm this conclusion. In any case, the coexistence of antiferromagnetism and superconductivity in RuSr2GdCu2O8 is clearly allowed.

## Acknowledgements

Work supported by the U. S. Department of Energy, Office of Science (at Northwestern University under Grant No. DE-F602-88ER45372 and at Argonne National Laboratory under Contract No. W-31-109-ENG-38) and by a grant of computer time at NERSC.


text:
# Hadron Structure Studied with the Electromagnetic Probe – from Giant Resonances to Meson Production

## 1 Introduction

In this talk, which I dedicate to the memory of Michael Danos, I would like to give a brief overview on a number of important issues in the field of photonuclear physics which have deepend considerably our understanding of nuclear structure over the past 40 years. It is not possible to give here a comprehensive review, and I had to make a selection which is governed by the scientific work of Danos, who has made many important contributions to nuclear photoreactions, and also by my own research interests. A very good represention of benchmark papers in this field may be found in . Let me remind you that in nuclear physics the electromagnetic probe has always been a very important tool. This is illustrated by a few highlights from the early development of nuclear physics:

* 1909: Coulomb scattering of $`\alpha `$-particles from a gold foil by Geiger and Marsden  from which Rutherford concluded in 1911 that most of the mass of an atom is concentrated in a tiny, almost pointlike nucleus in its center .
* 1934: Photodisintegration of the deuteron by Chadwick and Goldhaber  marks the beginning of photonuclear physics.
* 1947: Discovery of the giant dipole resonance by Baldwin and Klaiber  as a collective phenomenon almost exhausting the Thomas-Reiche-Kuhn dipole sum rule completely.
* 1951: First electron scattering experiment by Lyman, Hanson and Scott  establishing electron scattering as a very important tool for nuclear structure studies.

The important role of Danos’s research in the field of photonuclear physics is best illustrated by his influential 1961-Lectures at the University of Maryland . In these lectures Danos covered quite a variety of different topics: properties of the e.m. field, multipole decomposition, quantization, gauge invariance, Siegert’s theorem, dispersion relations, TRK and other sum rules, to name a few.

## 2 The dynamic collective model of the giant resonances

The collective phenomenon of the giant dipole resonance (GDR) can well be explained in the framework of the hydrodynamical model of Steinwedel and Jensen as an oscillation of a proton fluid against a neutron fluid. In particular, the dependence of the position of the GDR on the mass number is governed by the fact that for a spherical nucleus the eigenfrequency is proportional to the nuclear radius. It was the important observation of Danos (and independently of Okamoto) that for an axially symmetric deformed nucleus the GDR will be split into two peaks corresponding to oscillations along and perpendicular to the symmetry axis .

In view of the additional collective surface degrees of freedom, Danos and Greiner  developed in 1964 a unified dynamic collective model of the giant resonances (DCM) which includes in particular the coupling between the rotation-vibration and the giant resonance d.o.f. leading to additional dynamic effects. A comparison between the predictions for the total $`\gamma `$-absorption cross section and experimental data is shown in Fig. 1 for 165Ho including the giant quadrupole resonance contributions. The very good agreement is evident.

An important consequence of this dynamic coupling is the considerably strong dipole transition strength from the GDR states to the low lying rotational and vibrational states leading to sizeable Raman scalar and tensor scattering into these states (see Fig. 2, left panel) as experimentally observed by Fuller and Hayward . Another interesting feature of the DCM is the fact that a deformed nucleus with a nonvanishing ground state spin becomes optically anisotropic (nonvanishing tensor polarization) and thus its absorption and scattering cross sections depend on the nuclear orientation. This is shown in the right panel of Fig. 2 for 165Ho with spin $`I=7/2`$. A review may be found in .

Subsequently, the DCM was further developed in order to describe a much more general class of potential energy surfaces for the low energy collective d.o.f. allowing a unified description of the GDR for nuclei with quite different collective characteristics . An example is shown in Fig. 3.

## 3 Subnuclear degrees of freedom

In his 1961-Lectures Danos had already alluded to the possibility that a nucleon, in view of its extended internal structure, could be slightly deformed in a nuclear medium. Indeed, at the end of the 60s several groups, including Danos, H.T. Williams and myself, were developing a phenomenological model of such modifications by admixing into the nuclear wave function configurations, where one or several nucleons are internally excited, say as a $`\mathrm{\Delta }(1232)`$-resonance . Such wave function components have been coined “nuclear isobar configurations” (IC). Reviews may be found in . The characteristic features of these IC is that they possess small probabilities and a short range structure. Examples are shown in Fig. 4.

These subnuclear d.o.f., as manifest in isobar admixtures in nuclear wave functions, contribute also to the e.m. interaction. Alternatively, their contribution may be described by effective two-body exchange currents involving IC as shown diagrammatically in the left panel of Fig. 5 together with similar current contributions (MEC) from meson d.o.f. mediating the strong interaction, shown in the right panel for $`\pi `$-exchange. The largest MEC effects are found in $`E1`$-transitions which, however, are largely covered by the Siegert operator . Particularly strong MEC contributions to $`M1`$ are found in deuteron electrodisintegration near threshold, for which I show the inclusive cross section in Fig. 6. One readily notices the sizeable increase of the cross section by MEC and IC, particularly dominant at higher momentum transfers (see also ).

As last example for the importance of subnuclear d.o.f. I show in Fig. 7 differential cross sections of deuteron photodisintegration in the $`\mathrm{\Delta }`$-resonance region, where $`\mathrm{\Delta }`$-d.o.f. dominate. For a long time, the most sophisticated theoretical approaches were unable to describe properly the experimental date in this energy region. Only recently, a careful analysis has shown, that a proper treatment of meson retardation is absolutely necessary for a satisfactory description of experimental data .

## 4 Relativistic effects in the GDH sum rule

The Gerasimov-Drell-Hearn sum rule (1965/66) links a ground state property, the anomalous magnetic moment, to the energy weighted integral from threshold up to infinity over the spin asymmetry $`\sigma ^P(k)\sigma ^A(k)`$, the difference of the total photoabsorption cross sections for circularly polarized photons on a target with spin parallel and antiparallel to the spin of the photon, i.e.

$`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dk}{k}}\left(\sigma ^P(k)\sigma ^A(k)\right)`$ $`=`$ $`\mathrm{\hspace{0.17em}\hspace{0.17em}4}\pi ^2\kappa ^2{\displaystyle \frac{e^2}{M_t^2}}I,`$
where $`I`$ denotes the spin of the particle, $`M_t`$ its mass, and $`\kappa `$ its anomalous magnetic moment, defined by the total magnetic moment $`\stackrel{}{M}=(Q+\kappa )\frac{e}{M_t}\stackrel{}{S}`$, with $`eQ`$ as its charge and $`\stackrel{}{S}`$ its spin operator. An important consequence is, that a particle has to possess an internal structure if $`\kappa 0`$.

Here I will consider the GDH sum rule and in particular the spin asymmetry for the deuteron  in order to illustrate the fact, that even at low energies relativistic effects can become quite important in certain polarization observables. Application of the GDH sum rule to the deuteron reveals a very interesting feature. On the one hand, one expects a very small anomalous magnetic moment for the deuteron, because the deuteron has isospin zero ruling out the contribution of the large nucleon isovector anomalous magnetic moments to the deuteron magnetic moment. In fact, the experimental value is $`\kappa _d=.143`$ resulting in a GDH prediction of $`I_d^{GDH}=0.65\mu `$b, which is more than two orders of magnitude smaller than the nucleon values. On the other hand, one has the following absorptive processes: (i) photodisintegration $`\gamma +dn+p`$, (ii) single pion production (coherent: $`\gamma +dd+\pi ^0`$, and incoherent: $`\gamma +dN+N+\pi `$), (iii) two pion production etc. The processes (ii) and (iii) are dominated by quasifree production and, therefore, one may estimate from them a positive GDH contribution of the order of the sum of proton and neutron, namely $`I_p^{GDH}(\mathrm{})+I_n^{GDH}(\mathrm{})=438\mu \text{b}`$. In order to obtain the small total deuteron GDH value one, therefore, needs a large negative contribution of about the same size for cancellation.

Indeed, from the photodisintegration channel, which is the only photoabsorption process below the pion production threshold, a sizeable negative contribution arises at very low energies near threshold from the $`M1`$-transition to the resonant $`{}_{}{}^{1}S_{0}^{}`$ state, because this state can only be reached if the spins of photon and deuteron are antiparallel, and is forbidden for the parallel situation. The photodisintegration channel has been evaluated within a nonrelativistic framework but with inclusion of the most important relativistic contributions. The results are summarized in Fig. 8, where the spin asymmetry and the GDH integral as function of the upper integration limit is shown. One readily notices a very drastic decrease of the spin asymmetry at quite low energies by relativistic contributions, here mostly by the leading order relativistic spin-orbit current. At about 500 MeV convergence of the GDH-integral is achieved for the photodisintegration channel yielding a sume rule contribution of $`413`$ $`\mu `$b cancelling almost completely the estimated contribution from meson production. Without the relativistic contributions, one would have found a GDH contribution from photodisintegration of $`619`$ $`\mu `$b, which would lead to a negative total GDH-value for the deuteron. This demonstrates very convincingly the importance of relativistic effects in the spin asymmetry.

## 5 Eta photoproduction on the deuteron

Finally, I will briefly consider photoproductions of mesons on the deuteron. Meson photoproduction is the primary absorptive process on the nucleon and the interest of meson production on nuclei arises from the possibilities (i) to study the elementary neutron amplitude, (ii) to study the meson-nucleon interaction, (iii) to study medium effects, and (iv) to study nuclear structure.

As an example of recent interest, I will present some very new results on incoherent $`\eta `$-meson production on deuterium near threshold. It turns out that the simple impulse approximation fails badly to describe the experimental data. This is seen very clearly in Fig. 9 where the total cross section is shown. The strong enhancement by final state interaction (FSI) is most evident. In fact, it is so strong that the first order rescattering contribution underestimates FSI considerably close to threshold, thus making a three-body calculation necessary which leads to a much improved description, although not completely satisfactory. But in this calculation, FSI has been included in $`S`$-waves only. A more detailed test would be provided by a measurement of the meson spectrum shown in Fig. 10 which exhibits very clearly the signature of the virtual $`{}_{}{}^{1}S_{0}^{}`$-resonance of $`NN`$-scattering near threshold.

## 6 Concluding remarks

These examples should suffice for illustrating that the field of photo- and electronuclear physics has seen quite an impressive development both in breadth and depth over the past 40 years. This was made largely possible by the advent of new experimental tools and techniques and in theory by new ideas and the extraordinary progress in computing power. As a result, we have gained much deeper insights into the structure of hadrons, but still quite a few problems remain unresolved, in particular in the high energy regime, where a transition to quark-gluon dynamics is expected.

At present, the main research in this field is focused on the following topics:

* short range structure of nuclei, in particular, the nature of 2-body correlations,
* the size of relativistic effects in nuclear wave functions and operators,
* the role of subnuclear d.o.f. in terms of nuclear IC and meson exchange operators and their relation to the underlying quark-gluon d.o.f. of QCD,
* study of polarization observables,
* meson production on light nuclei.