text:
# Grains in Photo-Ionized Environments

## 1. Introduction

Grains are ubiquitous in the interstellar medium (ISM), and they can be detected either directly through their far-infrared emission or indirectly through extinction or polarization studies. Despite the vast number of observations, many questions regarding grain composition and grain physics remain unanswered. Further study is therefore required, and detailed models are needed to interpret the results. Ever since the pioneering study of Spitzer (1948), it has been widely recognized that grains play an important role in the heating and cooling of the diffuse ISM (see also the more recent studies by Bakes & Tielens 1994, and Weingartner & Draine 2001a, hereafter WD). Grains also play an important role in the physics of H ii regions and planetary nebulae (e.g., Maciel & Pottasch 1982, Baldwin et al. 1991, hereafter BFM, Borkowski & Harrington 1991, and Volk, these proceedings) and photo-dissociation regions (PDR’s, e.g., Tielens & Hollenbach 1985).

The interactions of grains with their environment include a host of microphysical processes, and their importance and effects can only be judged by including all of these processes. This can, in turn, only be done with a complete simulation of the environment. In this paper we consider in detail the physics of grains in both ionized and neutral environments, and model these with the spectral synthesis code Cloudy. We will present a calculation of photo-electric heating and collisional cooling rates for a range of physical conditions and grain materials and for a range of grain sizes (including a realistic size distribution). A comparison of our results with benchmark calculations using the WD code will also be presented.

## 2. The Photo-Ionization Code Cloudy

Cloudy is a well known and widely used photo-ionization code which is publicly available at http://www.pa.uky.edu/gary/cloudy. This code is not only useful for modeling fully ionized regions, but calculations can also be continued into the PDR. In order to make such a calculation realistic, the presence of a detailed grain model is required. The first grain model was introduced to Cloudy in 1990 to facilitate more accurate modeling of the Orion nebula (for a detailed description see BFM). In subsequent years, this model has undergone some revisions and extensions, but remained largely the same.

Recently, Cloudy has undergone several major upgrades, described in Ferland (2000a), Ferland (2000b), and van Hoof, Martin, & Ferland (2000). Some of these recent upgrades were aimed at improving the accuracy and versatility of Cloudy as a PDR code. These comprise an upgrade in the collision strengths of the \[C i\] and \[O i\] infrared fine-structure lines, the inclusion of a full CO model, and a comprehensive upgrade of the grain model. The latter was necessary for two reasons. First, the discovery of crystalline silicates in stellar outflows (e.g., Waters et al. 1996), and other detailed observations of grain emission features by the Infrared Space Observatory (ISO), meant that the code had to become much more flexible to allow such materials to be included in the modeling. Second, even before the ISO mission it had already become clear that the photo-electric heating and collisional cooling of the gas surrounding the grains is dominated by very small grains (possibly consisting of polycyclic aromatic hydrocarbons or PAH’s). The charging of very small grains could not be modeled very accurately with the original grain model (grain physics becomes increasingly non-linear as a function of charge for smaller grains and the average potential approach of the original grain model eventually breaks down for molecule-sized grains). In view of these facts we have undertaken a comprehensive upgrade of the grain model in Cloudy. The two main aims were to make the code more flexible and versatile, and to make the modeling results more realistic. The new grain model is included in version 96 of Cloudy.111A beta release is currently available at http://nimbus.pa.uky.edu/cloudy/cloudy\_ 96.htm

## 3. Resolving the Grain Size Distribution

In the original grain model, opacities for a handful of grain species were hard-wired in the code. Furthermore, only a single size bin would be used for the entire size distribution. This is not a very good approximation since most grain properties depend strongly on size, e.g., the grain opacities (see the right panel of Figure 1). This approach was nevertheless adhered to in BFM because of computational restrictions.

Resolving the size distribution into many small bins improves the modeling in several ways. First, an equilibrium temperature can be calculated for each bin separately (see the left panel of Figure 1). This is important since grain emissions are a strongly non-linear function of temperature.

More importantly, resolving the size distribution also enables other improvements of the grain treatment: stochastic heating can be treated correctly for the smallest grains in the size distribution (see § 4.), and the calculations yield much more accurate results for the total photo-electric heating and collisional cooling rates of the gas by the grains (see § 5.).

To improve the model, we have implemented the following changes:

1 – We have included a Mie code for spherical particles in Cloudy. Assuming that the grains are homogeneous spheres with a given complex refractive index (optical constant) one can use Mie theory (Mie 1908) to calculate the absorption and scattering opacity. This has to be done separately for every wavelength since the refractive index depends on wavelength. Good overviews of Mie theory can be found in van de Hulst (1957), and Bohren & Huffman (1983). Our Mie code is based on the program outlined in Hansen & Travis (1974) and references therein. The optical constants needed to run the code are read from a separate file. This allows greater freedom in the choice of grain species. Files with optical constants for a range of materials are included in the Cloudy distribution. However, the user can also supply optical constants for a completely different grain type.

2 – It is possible to use arbitrary grain size distributions. The user can either choose one of a range of preset functions (with numerous free parameters), or supply the size distribution in the form of a table. Single-sized grains can also be treated.

3 – The size distribution can be resolved in an arbitrary number of size bins (set by the user), and the absorption and scattering opacities and all the physical parameters (charge, temperature, etc.) are calculated for each bin separately.

### 3.1. The Dust Emission Spectrum in Orion

In Figure 2 we display the dust emission from the face of the Orion blister illuminated by $`\theta ^1`$ Ori C, using the BFM model with resolved and unresolved grain size distributions. These calculations include grains in the ionized region, the PDR, and fully molecular regions. Stochastic heating effects were turned off for these calculations to highlight the effect of resolving the size distribution. These effects would have been small anyway since this model is based on the BFM Orion size distribution which does not contain grains smaller than 30 nm. One can see that resolving the size distribution has a noticeable effect on the emitted spectrum since the grain emissivity depends strongly on size.

## 4. Stochastic Heating of Small Grains

It is well known that in conditions where the cooling time of the grain is shorter than or comparable to the average time between two significant heating events, a stochastic treatment of the grain temperature is necessary. This effect is important for grains smaller than roughly 20 nm (most notably for PAH’s), and/or in regions where the photon density is very low (e.g., the diffuse ISM). The importance of this effect was first suggested by Greenberg (1968), and the topic has evolved considerably since (e.g., Desert el al. 1986, Dwek 1986), mainly focusing on efficient numerical techniques for calculating the temperature distribution of the grains. A brief history can be found in Guhathakurta & Draine (1989, hereafter GD).

Code implementing stochastic treatment of PAH’s was already included in a revision to the old grain model. However, that code only worked with the two PAH species that were hardwired in Cloudy, and therefore was only infrequently used. In order to obtain accurate modeling results, the code should work on very small grains in any size distribution (which is only possible now that the size distribution has been resolved), and should include the effects automatically when they have a noticeable effect on the emitted spectrum222In Cloudy version 96beta2 stochastic heating effects still need to be turned on explicitly by the user. Once the numerical stability of the new algorithm has been sufficiently validated, the situation will be reversed and stochastic heating effects will be included automatically unless explicitly turned off by the user.. To achieve these goals, the stochastic heating code has been extensively rewritten for the new grain model. It now works efficiently with all grain types and sizes, and under all conditions.

The original code in Cloudy was based on the algorithm of GD. The aim of this algorithm is to calculate the probability distribution of grain temperatures (or equivalently: grain enthalpies). When the code was upgraded for the current release of Cloudy, certain improvements were added to this algorithm for better performance.

The first difference with the work of GD is that we assumed $`\lambda _{\mathrm{cutoff}}=\mathrm{}`$ and dropped the term containing the integrated incident flux for wavelengths longer than $`\lambda _{\mathrm{cutoff}}`$ from the energy balance. The quantity $`\lambda _{\mathrm{cutoff}}`$ has no physical meaning, and was merely a numerical invention introduced by GD to avoid dealing with zero temperature grains. Our algorithm incorporates a different and more natural solution for that problem. This improved convergence of the probability distribution dramatically and helped us to avoid calculating grain enthalpy bins with extremely low probabilities. The algorithm of GD was rewritten in such a way that it is no longer necessary to set up the grain enthalpy grid a priori, but use an adaptive stepsize algorithm instead. This makes the algorithm much more flexible and efficient, while still guaranteeing proper convergence. The code can automatically detect when stochastic heating might be important, or when it is safe to simply use the equilibrium temperature for calculating the emitted spectrum. To highlight the importance of stochastic heating, we show in Figure 3 the spectrum emitted by a 5 nm silicate grain in typical diffuse ISM conditions with stochastic heating effects either included or not.

## 5. Changes to the Grain Physics

We have modified certain aspects of the grain physics following the discussion in WD. Below we highlight certain aspects of these changes. A detailed description will be presented in a forthcoming paper (van Hoof et al., in preparation).

1 – We include the bandgap between the valence and conduction bands in our potential well model for silicates. This change only affects the results for negatively charged grains ($`Z_\mathrm{g}1`$).

2 – Reduction of the potential barrier for negatively charged grains is included using an analytic fit to numerical calculations. Two effects are important here: quantum tunneling and the Schottky effect. Quantum theory predicts that an electron with insufficient energy to overcome a barrier still has a finite chance of tunneling through. This effect has been modeled using the WKB approximation which gives a simple analytic expression for the tunneling probability for a barrier of given width and height. Quantum tunneling is only important for small grains. For large grains the Schottky effect will dominate, which describes the lowering of the potential barrier by an image potential in the grain. This effect has been accurately modeled by Draine & Sutin (1987).

We will approximate both effects by assuming that the barrier is effectively reduced in height from $`(Z_\mathrm{g}+1)e^2/(4\pi ϵ_0a)`$ to $`E_{\mathrm{min}}`$. The magnitude of the combined tunneling/Schottky effect was calculated by WD. However, the fitting function they used has the wrong limiting behavior for large grains where it should asymptotically approach the classical Schottky expression. We therefore repeated these calculations using the same assumptions, but adopted a different fitting function that does exhibit the correct limiting behavior:

$$E_{\mathrm{min}}=\theta _\nu \frac{e^2}{4\pi ϵ_0a}\left[1\frac{0.3}{(a/\mathrm{nm})^{0.45}\nu ^{0.26}}\right],$$
(1)
where $`\theta _\nu [\nu =(Z_\mathrm{g}+1)]`$ describes the Schottky effect and is defined in Draine & Sutin (1987). The term in square brackets describes the quantum-mechanical correction. This change only affects the results for grains with $`Z_\mathrm{g}<1`$.

3 – The treatment of the photo-electric effect has been improved, following the discussion in WD. This includes new expressions for the photo-electric yield and the energy distribution of ejected electrons.

4 – The treatment of electron sticking probabilities has been updated, again following WD. Especially for very small grains the sticking efficiency has been substantially lowered to obtain better agreement with laboratory studies of molecules. This has an important impact on the photo-electric heating rate of the gas since the electron recombination rate has to be matched by electron loss processes to preserve the charge balance. The loss processes are usually dominated by the photo-electric effect.

5 – Certain physical constants have been updated. Most notably the work function for graphite has been lowered. This results in an increased photo-electric heating rate of the gas since less of a potential barrier needs to be overcome to ionize graphite grains.

6 – The treatment of collisional heating or cooling of the grains and the gas has been improved.

Our treatment deviates from the WD code in two ways. Most importantly, we use a different grain charge model, which will be discussed in more detail in § 5.1. Secondly, we use slightly different physics for collisions between ions and grains. This only gives rise to very small differences at the 1 – 2% level when compared to the WD treatment.

### 5.1. The New Hybrid Grain Charge Model

The original grain model in Cloudy (which we will call the average grain potential model) is described in BFM. In that model an average grain potential is calculated by finding the potential for which the charge gain rate exactly matches the loss rate. This method was first proposed by Spitzer (1948), and is an excellent approximation for large grains. However, it is now clear that photo-electric heating and collisional cooling of the gas are dominated by very small grains. For such grains the average grain potential approximation does not work very well because grain physics becomes increasingly non-linear as a function of charge for smaller grain sizes. This fact, combined with the fact that grain charges are quantized, has led to a new approach where the charge distribution is fully resolved, and heating and cooling rates are calculated for each charge state separately (see e.g., WD). This ensures accurate results, but leads to an appreciable increase in computational overhead. This is especially the case for large grains since the width of the charge distribution increases with grain size. Hence the paradoxical situation arises that most of the computing time is spent on grains which contribute least to heating and cooling, and which are also the grains for which the average grain potential model works best!

In this paper we present a hybrid grain potential model which saves most of the computational speed of the original average grain potential model, but nevertheless gives sufficient accuracy when compared to fully resolved charge distribution calculations. The basic philosophy is that for very small grains ($`a<1`$ nm) only a few charge states have a significant population. Hence we adopt the $`n`$-charge state approximation, in which all grains are treated by using exactly $`n`$ contiguous charge states, independent of size. The higher $`n`$ is, the more accurate the results will be (exactly how accurate will be discussed in § 5.2.). The default for Cloudy calculations is $`n=2`$, but the user can request a larger number if higher precision is desired.

Since the $`n`$-charge state model does not fully resolve the charge distribution, a different algorithm from WD is needed to calculate the fractional populations $`f_i`$ of each of the $`n`$ charge states $`Z_iZ_1+i1`$. These populations must first of all obey the following normalization:

$$\underset{i=1}{\overset{n}{}}f_i=1.$$
(2)
Secondly, we require that the electron gain rates $`J_i^{}`$ and electron loss rates $`J_i^+`$ summed over all charge levels match exactly:

$$\underset{i=1}{\overset{n}{}}f_i(J_i^+J_i^{})=0,$$
(3)
similar to the average grain potential model. Eqs. 2 and 3 are sufficient to determine the charge state populations if $`n=2`$, but for $`n>2`$ we need additional equations. These equations need to satisfy the following constraints. First, the resulting level populations should always be greater or equal zero. Second, the level populations should change continuously when the electron gain and loss rates change continuously. Third, the level populations should asymptotically approach the results from fully resolved calculations for increasing values of $`n`$. We have adopted the following algorithm:

1 – The $`n`$ charge states are split up in two groups of $`n1`$ contiguous charge states. The first group contains charge states $`[Z_1,Z_{n1}]`$, and the second $`[Z_2,Z_n]`$. The value for $`Z_1`$ is determined iteratively (see step 4).

2 – The relative level populations in the first group $`f_i^1`$ are determined using an algorithm very similar to the one used in fully resolved calculations, i.e., assume $`f_1^1=1`$, calculate $`f_2^1=f_1^1J_1^+/J_2^{}`$, $`f_3^1=f_2^1J_2^+/J_3^{}`$, etc.333Note that this procedure is not correct for charge transfer with multiply charged ions. In order to avoid having to solve a full set of linear equations, we will approximate this process as multiple single-charge-transfer events. The resulting errors are expected to be small as collision rates for multiply charged ions are usually quite low because their velocities are small compared to electrons. The collision rates are normally suppressed even further by the positive grain charge., and then re-normalize to $`_{i=1}^{n1}f_i^1=1`$. We use an analogous procedure for the populations $`f_i^2`$ of the second group for $`i[2,n]`$.

3 – Determine for both groups the net charging rate

$$J_k=\underset{i=k}{\overset{n2+k}{}}f_i^k(J_i^+J_i^{})(k=1,2).$$
(4)
4 – Iterate $`Z_1`$ and repeat steps 1 – 3 until $`J_1\times J_20`$. Then find $`0\alpha 1`$ such that

$$\alpha J_1+(1\alpha )J_2=0.$$
(5)
5 – Determine the final charge state populations as follows

$$f_i=\alpha f_i^1+(1\alpha )f_i^2(\text{with }f_n^10,f_1^20).$$
(6)
One can verify that this algorithm satisfies all constraints.

The hybrid grain potential model is efficient because it avoids the overhead for large grains, while still giving accurate results for both small and large grains. An added bonus is that most of the time an excellent initial estimate for $`Z_1`$ can be derived from the previous zone, reducing the overhead even further. The model works for very small grains because only few charge states are populated and it can reconstruct the actual charge distribution. It works for large grains because the grain potential distribution approaches a delta function for increasing grain size (as opposed to the charge distribution which becomes ever wider). Our method therefore asymptotically approaches the average grain potential model, which we already know is very accurate for large grains.

### 5.2. Validating the Hybrid Grain Charge Model

In order to validate the new grain charge model, we calculated the photo-electric heating and collisional cooling rates for a range of physical conditions, two grain species, and a wide range of grain sizes (including a realistic size distribution). We then compared these calculations with benchmark results from the WD code, which fully resolves the charge distribution.

We model the physical conditions with simple assumptions: the plasma only contains hydrogen, the electron temperature and density are fixed at prescribed values, the incident spectrum is assumed to be a blackbody (either full in the warm ISM and ionized cases, or cut off at 13.6 eV in the cold ISM and PDR cases). To allow for an accurate comparison, both codes used the same optical constants for classical graphite and astronomical silicate (Draine & Lee 1984). In Cloudy the photo-electric heating by local diffuse emission was switched off, as well as thermionic emissions and collisional cooling by neutral hydrogen since these are not treated in the WD code. The WD code was modified according to Eq. 1. Except for the grain charge model, the treatment of the grain physics in Cloudy and the WD code are very similar.

We present results for selected single sized grains, as well as a realistic reconstruction of the ISM grain size distribution (labeled A6) taken from Weingartner & Draine (2001b). In the ISM and H ii region cases we used a standard mix of graphite and silicate as prescribed by Weingartner & Draine (2001b)444Note that the A6 size distribution was derived by matching the extinction curve using grain opacities defined in Li & Draine (2001). In this study we use the same size distribution, but use opacities from Draine & Lee (1984) instead, which is not consistent. This inconsistency is irrelevant for our purposes.. In the planetary nebula cases, we made separate models for either graphite or silicate as these two species are not expected to coexist spatially in the nebular material on theoretical grounds555Note that spectroscopy of planetary nebulae by ISO has revealed a surprisingly large number of cases that show both silicate and graphitic dust features. One famous example is NGC 6302 (Molster et al. 2001). It is however usually assumed that the two species reside in different parts of the nebula..

The physical parameters for the models are listed in Table 1; the results of the calculation are presented in Tables 2, 3, and 4.

### 5.3. Discussion

A summary of the comparison of the photo-electric heating and collisional cooling rates from Tables 2, 3, and 4 is given in Table 5. In general the results are in excellent agreement, with only a few outliers for single sized grains in the $`n=2`$ and $`n=3`$ cases. The results for the size distribution cases always agree to better than 25%, even for $`n=2`$. This is well within the accuracy with which we know grain physics to date. There are still major uncertainties in the photo-electric yields and the sticking efficiency for electrons, both of which have a strong effect on the photo-electric heating rate. Also the work function and bandgap for astrophysical grain materials are poorly known and can have a strong effect as well. This is unfortunate since photo-electric heating and collisional cooling in photo-ionized environments are important effects. All these uncertainties mainly stem from the fact that the composition of interstellar grains is still poorly known.

From Table 5 one can see that the collisional cooling rates usually are in better agreement than the photo-electric heating rates for a given set of physical parameters. It is furthermore clear that the accuracy of the $`n`$-charge state approximation increases as $`n`$ increases, as should be expected. Closer inspection of Tables 2, 3, and 4 for single sized grains reveals that the largest errors in the $`n=2`$ and $`n=3`$ cases are for the 0.5 nm grains, for $`n=4`$ for either 2 nm or 10 nm grains, and for $`n=5`$ for 10 nm grains, i.e., the grain size for which the errors are largest shifts upwards for higher values of $`n`$. This is expected since the $`n`$-charge state approximation will fully resolve the charge distribution of the smallest grains for $`n>3`$. Note that the results for the 100 nm grains are always in excellent agreement, even when $`n=2`$, despite the fact that the actual charge distribution is much wider than that.

The agreement between the Cloudy and the WD results is very satisfactory for realistic size distributions, and should be sufficient for all realistic astrophysical applications. Therefore the hybrid grain charge model presented above (with $`n=2`$) will be the default for Cloudy modeling.

## 6. Outlook

The rapid increase in computing power over the past decades has made it possible now to build detailed grain models that can be run in a modest amount of time. Computing power is therefore not a limiting factor anymore, and the gaps in our knowledge of grain physics are almost exclusively what is hampering progress at the moment. These gaps need to be filled in, both through laboratory work, and observations of interstellar and circumstellar grains. Some of the key areas have already been indicated above. They include a better understanding of the exact chemical composition of the grains, and the size distribution as a function of environment, as well as the corresponding physical properties (electron and ion sticking efficiencies, photo-electric yields, work functions, and bandgaps) for realistic astronomical materials. Depletion of trace elements onto grains should also be understood better since that has an important effect on the abundances of the gas in which the grains have been formed and alters the emitted spectrum and cooling rate of the gas. An improved model of the life cycle of dust (including grain formation, coagulation, and destruction in the ISM, as well as the formation of crystalline materials in circumstellar envelopes, which was one of the great surprises from the ISO mission) should go hand in hand with the above improvements. In the long run theory combined with laboratory experiment should move towards building full quantum-mechanical models of very small grains, which are now treated by heuristic, semi-classical methods. Currently the complexity of such models is still computationally prohibitive.

We are currently using Cloudy to model the Ney-Allen nebula (Ney & Allen 1969), which is situated close to $`\theta ^1`$ Ori D. The silicate emission from this nebula was used (combined with data from long period variable stars) as a template to define astronomical silicate (Draine & Lee 1984). Now this effort will come full circle by using the derived astronomical silicate data to model the Ney-Allen nebula. This will be an important test of our understanding of interstellar grains. The improvements presented above are important for this effort since the silicate 10 $`\mu `$m feature is situated in the Wien tail of the emission, and effects from resolving the size distribution (see Fig. 2) as well as temperature spiking of very small grains are important. This work will be presented in a forthcoming paper (van Hoof et al., in preparation).

#### Acknowledgments.

We wish to acknowledge financial support by the National Science Foundation through grant no. AST–0071180, and NASA through its LTSA program, NAG 5–3223. This research was also supported in part by the Natural Sciences and Engineering Research Council of Canada. J.C.W. acknowledges support from an NSF International Research Fellowship.

## References

Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822

Baldwin, J. A., Ferland, G. J., Martin, P. G., Corbin, M. R., Cota, S. A., Peterson, B. M., & Slettebak, A. 1991, ApJ, 374, 580 (BFM)

Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (Ney York: Wiley)

Borkowski, K. J., & Harrington, J. P. 1991, ApJ, 379, 168

Desert, F. X., Boulanger, F., & Shore, S. N. 1986, A&A, 160, 295

Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803

Dwek, E. 1986, ApJ, 302, 363

Ferland, G. J. 2000a, in Astronomical Data Analysis Software and Systems IX, eds. Manset, N., Veillet, C., & Crabtree, D., ASP Conference Series, Vol. 216, (San Francisco: Astronomical Society of the Pacific), p. 32.

Ferland, G. J. 2000b, Rev. Mex. Astr. Astrof. (Serie de Conferencias), 9, 153

Greenberg, J. M., 1968, in Stars and Stellar Systems, Vol. 7, Nebulae and Interstellar Matter, eds. Middlehurst, B. M., & Aller, L. H., (Chicago: Chicago University Press), p. 221

Guhathakurta, P., & Draine, B. T. 1989, ApJ, 345, 230 (GD)

Hansen, J. E., & Travis, L. D. 1974, Space Sci. Rev., 16, 527

Li, A., & Draine, B. T. 2001, ApJ, in press (astro-ph/0011319)

Maciel, W. J., & Pottasch, S. R. 1982, A&A, 86, 380

Martin, P. G., & Rouleau, F. 1991, in Extreme Ultraviolet Astronomy, eds. Malina, R. F., & Bowyer, S., (New York: Pergamon Press), p. 341

Mie, G. 1908, Ann. Phys., 25, 377

Molster, F. J., et al. 2001, A&A, 372, 165

Ney, E. P., & Allen, D. A. 1969, ApJ, 155, L193

Spitzer, L. 1948, ApJ, 107, 6

Tielens, A. G. G. M., & Hollenbach, D. J. 1985, ApJ, 291, 722

van de Hulst, H. C. 1957, Light Scattering by Small Particles, (Ney York: Wiley)

van Hoof, P. A. M., Martin, P. G., & Ferland, G. J. 2000, in Cosmic Evolution and Galaxy Formation: Structure, Interactions, and Feedback, eds. Franco, J., Terlevich, E., López-Cruz, O., & Aretxaga, I., ASP Conference Series, Vol. 215, (San Francisco: Astronomical Society of the Pacific), p. 220 (extended version: astro-ph/0001196)

Waters, L. B. F. M., et al. 1996, A&A, 315, L361

Weingartner, J. C., & Draine, B. T. 2001a, ApJS, 134, 263 (WD)

Weingartner, J. C., & Draine, B. T. 2001b, ApJ, 548, 296


format: text

text:
# Standing Shocks around Black Holes: An Analytical Study

## 1 Introduction

In recent years, study of standing and oscillating shocks in accretion flows has become very important since it is recognized that the spectral states of black holes as well as Quasi-Periodic Oscillations (QPOs) observed in light curves of black hole candidates are directly related to the radiative transfer properties of a compact Comptonizing region close to a black hole (e.g., Chakrabarti and Titarchuk 1995, CT95; Ling et al. 1997; Chakrabarti & Manickam, 2000, CM00; Muno et al., 2000; Feroci et al. 2000; Homan et al., 2001). These shocks have been studied most extensively in the literature (Chakrabarti 1989, C89; Chakrabarti 1990, C90) and their properties have been verified by several independent groups (Yang and Kafatos, 1995; Nobuta and Hanawa, 1994; Lu and Feng, 1997) of workers. If for a given set of initial parameters the standing shocks exists analytically, numerical simulations also would find them (Chakrabarti and Molteni, 1993; Molteni, Lanzafame and Chakrabarti, 1994; Ryu et al., 1995; Molteni, Ryu and Chakrabarti, 1996) otherwise, the shock would be oscillating (Molteni, Sponholz and Chakrabarti, 1996; Ryu, Chakrabarti and Molteni, 1997) causing QPOs in the X-rays emitted from the post-shock region. Furthermore, it has been observed that at least in some of the black hole candidates, such as GRS1915+105, outflows are produced from the region which also emits the Comptonized photons (Fender et al, 2000; Dhawan et al., 2000) i.e., the post-shock region or the centrifugal barrier dominated region according to our present understanding of the flow solutions.

It is therefore pertinent to ask if the properties of the shocks, such as location, strength, compression ratio etc. could be understood solely analytically. The solutions obtained so far by Chakrabarti and his collaborators and other groups always resort to numerical means. Roughly, the method was the following: for a given set of parameters (such as the specific energy and specific angular momentum), it is first determined if the flow allows more than one X-type sonic points (C89). This is because, at the horizon, matter must have velocity same as that of the velocity of light and for causality arguments, the flow must be supersonic. Thus, once the flow becomes supersonic through the outermost sonic point, and forms a shock (i.e., jumps to a subsonic branch) it must pass through the inner sonic point to become supersonic at the horizon. Next question is whether the specific entropy at the inner sonic point is higher compared to that at the outer sonic point. This is because, at the shock, entropy must be generated and post-shock flow which passes through the inner sonic must have higher entropy. The final and the most important question is whether all the Rankine-Hugoniot conditions (also known as the shock conditions) are satisfied somewhere between the two X-type sonic points. The three Rankine-Hugoniot conditions could be combined to obtain a combination of Mach numbers (Shock Invariant). C89 obtained this expression which is continuous across the shock and used this to obtain the shock location by iterative technique.

In the present paper, we follow the same philosophy, but obtain the shock locations analytically. There were two motivations to do this. First, from the theoretical point of view, it is challenging to find solutions of a large number of non-linear equations which must satisfy a number of conditions mentioned above. Second, from the observers’ point of view, any observation which could require standing shock waves, could be explained using more fundamental parameters, such as specific energy and angular momentum or even better, using accretion rates of Keplerian and sub-Keplerian flows as in a two-component flow solution of CT95. Analytical work also gives insight into why the shocks form in the first place. The boundary of the parameter space for shock formation is also obtained by analytical means. These findings are important as they would tell us when QPOs may or may not be seen. These details would be discussed in future.

In the next Section, we present the model equations and shock conditions. In §3, we present the sonic point analysis and write down the expression for Shock Invariant (C89) for the shake of completeness. In §4, we present the analytical expressions of the sonic points and discuss how the parameter space is divided into regions of one or three sonic points. In §5, we present expression of shock locations and again divide the parameter space into regions which may or may not have shocks. We compare our solutions with numerical work. In §6, we discuss some of the astrophysical implications of our solutions and finally, in §7, we make concluding remarks.

## 2 Model Equations and Shock Conditions

We start with a thin, axisymmetric, inviscid, steady flow which is in equilibrium in a direction transverse to the flow. The model equations which govern the motion of matter accreting into black hole are given as follows (C89):

A. Radial Momentum Equation:

$$\vartheta \frac{d\vartheta }{dx}+\frac{1}{\rho }\frac{dP}{dx}+\frac{\lambda ^2}{x^3}+\frac{1}{2(x1)^2}=0$$
(1)
In a non-viscous flow, integration of this equation leads to the Energy Conservation Equation :

$$=\frac{\vartheta _e^2}{2}+\frac{a_e^2}{\gamma 1}+\frac{\lambda ^2}{2x^2}+g(x)$$
(2)
where, $`g(x)`$ is the pseudo-Newtonian potential introduced by Paczyński & Wiita (1980) which is given by, $`g(x)=\frac{1}{2(x1)}`$. Here, $`\rho `$ is the mass density, $`P`$ is the isotropic pressure, $`\vartheta _e`$ and $`a_e`$ are the non-dimensional radial and the sound velocities measured in units of velocity of light $`c`$, $`x`$ is the non-dimensional radial distance measured in units of the Schwarzschild radius $`r_g=2GM/c^2`$, $`M`$ being the mass of the black hole and $`G`$ being the gravitational constant. $`\gamma `$ is the adiabatic index of the flow and $`P=K\rho ^\gamma `$ is assumed to the equation of state. $`K`$ is the measure of the specific entropy which is constant except at the shock location where local turbulence generate some entropy enabling the flow to pass through the inner sonic point. The subscript $`e`$ refers to the quantities measured on the equatorial plane.

B. Continuity Equation:

$$\frac{d}{dx}(\vartheta _e\rho xh)=0$$
(3)
which is integrated to obtain the Mass Conservation Equation :

$$\dot{M}=\vartheta _e\rho xh,$$
(4)
where, $`h`$ is the half-thickness of the flow at radial coordinate $`x`$. Here $`\dot{M}`$ is the mass accretion rate apart from a geometric constant. If we assume that the flow is in hydrostatic equilibrium in the transverse direction, the vertical component of gravitational force balances the pressure gradient force. Hence the expression for the half thickness of the disc is given by (C89),

$$h(x)=a_ex^{1/2}(x1).$$
(5)
We write the mass flux conservation equation in terms of $`\vartheta _e`$ and the sound speed on the equatorial plane $`a_e=\sqrt{\gamma P/\rho }`$ in the following way,

$$\dot{}=\vartheta _ea_e^qx^{3/2}(x1)=\vartheta _ea_e^qf(x)$$
$`(4a)`$
where, $`q=\frac{\gamma +1}{\gamma 1}`$ and $`f(x)=x^{3/2}(x1)`$.

The shock conditions which we employ here are the following (C89):

(a) the energy flux is continuous across the shock —

$$_+=_{},$$
(6)
(b) the mass flux is continuous across the shock —

$$\dot{M}_+=\dot{M}_{}$$
(7)
and finally, (c) the momentum balance condition —

$$W_++\mathrm{\Sigma }_+\vartheta _{e}^{2}{}_{+}{}^{}=W_{}+\mathrm{\Sigma }_{}\vartheta _{e}^{2}{}_{}{}^{}$$
(8)
where subscripts “$``$” and “$`+`$” refer, respectively, to quantities before and after the shock. A shock satisfying these conditions is termed as a Rankine-Hugoniot shock (Landau and Lifshitz, 1959). Here, $`W`$ and $`\mathrm{\Sigma }`$ denote the pressure and the density, integrated in the vertical direction \[see, e.g., Matsumoto et al. (1984)\], $`i.e.`$,

$$\mathrm{\Sigma }=\underset{h}{\overset{h}{}}\rho 𝑑z=2\rho _eI_nh,$$
(9)
and

$$W=\underset{h}{\overset{h}{}}P𝑑z=2P_eI_{n+1}h,$$
(10)
where, $`I_n=\frac{(2^nn!)^2}{(2n+1)!}`$, $`n`$ being the polytropic index as defined previously. In the subsequent analysis we drop the subscript $`e`$ if no confusion arises in doing so.

## 3 Sonic Point Analysis and Shock Invariants

In order to have a shock, the flow must be supersonic i.e., the stationary flow must pass through the sub-sonic flow to a supersonic flow. Discussions in this Section are based on early works \[see, Chakrabarti (1989, 1990)\].

### 3.1 Sonic Point Conditions

From the radial momentum equation and the continuity equation (Eqs. 1 and 3) we derive the sonic point condition (or critical point condition) in the usual way (C89). The first derivative of the radial velocity w.r.t radial distance is given by,

$$\frac{d\vartheta }{dx}=\frac{\frac{2a^2}{\gamma +1}\frac{dln(f)}{dx}\frac{dG}{dx}}{\left[\vartheta \frac{2a^2}{(\gamma +1)\vartheta }\right]}.$$
(11)
Here, $`G(x)=\frac{\lambda ^2}{2x^2}\frac{1}{2(x1)}`$ is the effective potential. Since the denominator must vanish at the sonic points, if the flow is assumed to be smooth everywhere, the numerator must also vanish simultaneously. The vanishing of the denominator gives,

$$\vartheta _c^2(x_c)=\frac{2}{(\gamma +1)}a_c^2(x_c).$$
(12)
The factor (which is unity only in isothermal flows) in front of $`a_c^2(x_c)`$ arises because the flow is assumed to be in vertical equilibrium. The vanishing of the numerator gives,

$$a_c^2(x_c)=\frac{(\gamma +1)(x_c1)}{x_c^2}\frac{[\lambda _K^2(x_c)\lambda ^2]}{(5x_c3)}.$$
(13)
The subscript $`c`$ denotes quantities at the critical points. Here $`\lambda _K`$ is the Keplerian angular momentum defined as $`\lambda _K^2=x_c^3/[2(x_c1)^2]`$. It is to be noted that since the left hand side is always positive, angular momentum at the sonic point must be sub-Keplerian, i.e., $`\lambda (x_c)<\lambda _K`$ (e.g., C90 and references therein).

### 3.2 Mach Number Relation at the Shock

From the equations given in §2, we now seek an invariant relation which must be satisfied at the shock (C89). We rewrite the condition of energy flux continuity (Eq. 6), and the pressure balance condition (Eq. 8) in terms of the Mach number $`M=\vartheta /a`$ of the flow,

$$\frac{1}{2}M_+^2a_+^2+\frac{a_+^2}{\gamma 1}=\frac{1}{2}M_{}^2a_{}^2+\frac{a_{}^2}{\gamma 1},$$
(14)
$$\dot{_+}=M_+a_+^\nu ^{^{}}f(x_s),$$
(15)
$$\dot{_{}}=M_{}a_{}^\nu ^{^{}}f(x_s),$$
(16)
where $`\nu ^{^{}}=\frac{2\gamma }{\gamma 1}`$, and

$$\frac{a_+^\nu }{\dot{}_+}\left(\frac{2}{3\gamma 1}+M_+^2\right)=\frac{a_{}^\nu }{\dot{}_{}}\left(\frac{2}{3\gamma 1}+M_{}^2\right),$$
(17)
where, $`\nu =\frac{3\gamma 1}{\gamma 1}`$ and $`x_s`$ is the location of the shock. From Eqs. (14-17) one obtains the following equation relating the pre- and post-shock Mach numbers of the flow at the shock (C89),

$$C=\frac{\left[M_+(3\gamma 1)+(2/M_+)\right]^2}{2+(\gamma 1)M_+^2}=\frac{\left[M_{}(3\gamma 1)+(2/M_{})\right]^2}{2+(\gamma 1)M_{}^2}.$$
(18)
The constant $`C`$ is invariant across the shock. The Mach number of the flow just before and after the shock can be written down in terms of $`C`$ as,

$$M_{}^2=\frac{2(3\gamma 1)C\pm \sqrt{C^28C\gamma }}{(\gamma 1)C(3\gamma 1)^2}.$$
(19)
The product of the Mach number is given by,

$$M_+M_{}=\frac{2}{\left[(3\gamma 1)^2(\gamma 1)C\right]^{1/2}}.$$
(20)
## 4 Analytical Expression of the Sonic Points and Behavior in Parameter Space

To obtain shock locations, we need to first obtain the locations of the sonic points, and ensure that at least two of them are X-type (C90). In §3.1, we presented the sonic point conditions. Using the definition $`M=\vartheta /a`$ of the Mach number, and substituting $`a_c^2(x_c)`$ from Eq. (13) into Eq. (2), we obtain following algebraic equation for $`x_c`$ given by,

$$𝒩x_c^4𝒪x_c^3+𝒫x_c^2𝒬x_c+=0.$$

where, $`𝒩=10`$; $`𝒪=16+2n3`$; $`𝒫=6+\lambda ^2(4n1)3`$; $`𝒬=8n\lambda ^2`$; $`=(1+4n)\lambda ^2`$ and $`n`$ $`(=\frac{1}{\gamma 1})`$ is the polytropic index.

We solve this equation analytically (Abramowitz and Stegun, 1970) and get the location of the sonic points. Details are given in Appendix A. For the purpose of critical or sonic points, $`D`$ of the Appendix would be denoted as $`D_c`$. The equation has four roots and Eq. A.8 can be used to check whether all of them are complex (pair of complex conjugates) or at least two are real and two are complex or all four are real. A necessary condition to form a shock wave is to have all four real roots. Of course, only one would be inside the black hole, and the other three would be outside, and out of these, due to topological reasons, only two would be X-type or saddle type and the one in between must be ‘O’-type or Center type. These are determined by computing derivative $`d\vartheta /dx`$ at the sonic point by using L’Hospital’s rule and checking if they are real. For our purpose, two derivatives at each sonic point must be real and of opposite signs in order that the sonic point be of X-type or saddle type.

Figure 1 shows the division of the parameter space. Denoting the discriminant $`D`$ (of Appendix A) by $`D_c`$ we find that the condition $`D_c<0`$ is the necessary condition to have three sonic points. The boundary $`D_c=0`$ separates this region on which two sonic points merge and the third one remains separate. Outside of this region $`D_c>0`$, and only one sonic point is possible and the other two roots are complex conjugate of each other. The dotted curve in the middle represents the condition:

$$\dot{}_o=\dot{}_i,$$

where $`\dot{}_i`$ and $`\dot{}_o`$ are the entropy accretion rates at the inner and outer sonic points respectively. The region above it contains parameters with $`\dot{}_i<\dot{}_o`$ and the region below it contains parameters with $`\dot{}_i>\dot{}_o`$. This latter region is suitable for shock formation in accretion flows.

In passing we wish to point out that from Eq. (11) one can easily show that locations where $`dM/dx=0`$ exactly coincide with the sonic points of the flow. Thus the number of extrema of $`M=M(x)`$ is the same as that of the sonic points and hence the division in Fig. 1 could give an idea about the behavior of $`M=M(x)`$ as well.

## 5 Analytical Expression for the Shock Location and Behavior in Parameter Space

A black hole accretion flow being transonic, it must satisfy two sonic point conditions at the cost of one extra unknown, namely, the the sonic point. Because of this extra condition, out of the three constants of motion, namely, $``$, $`\dot{M}`$ and $`\lambda `$, only two are to be supplied as free parameters. While computing shock locations, we used past experience derived from numerical methods, that only one of the shocks, namely $`x_{s3}`$ (C89 notation) is stable, and accordingly our procedure as delineated below attempts to compute only this location.

The flow will have a shock only at the point where the shock invariant condition is satisfied. Simplifying the shock invariant relation (Eq. 18), we obtain,

$$2(\gamma 1)(M_+^2+M_{}^2)[(3\gamma 1)^22(3\gamma 1)(\gamma 1)]M_+^2M_{}^2+4=0.$$
(21)
We consider a relativistic flow with $`\gamma =4/3`$ so that the polytropic index $`n=3`$. Then from Eq. (22), we obtain,

$$2(M_+^2+M_{}^2)21M_+^2M_{}^2+12=0.$$
(22)
We now expand the post-shock Mach number $`M_+^2`$ by a polynomial which must satisfy the following conditions:

(a) $`\frac{dM}{dx}`$ is zero at the central ‘O’-type sonic point. This is a general property of the flow (vide discussion at the end of last Section).

(b) Mach no. ($`M_+`$) at the location of the middle sonic point must match with that derived from approximate analytical solution obtained using the energy equation (Eq. 2 and 4a).

(c) Solution must pass through the position where $`\frac{dM}{dx}`$ is $`\mathrm{}`$ and a good guess of this the location (say, from the location of the sonic points) is known.

Similarly, we expand $`M_{}^2`$ by a polynomial which must satisfy the following conditions:

(a) $`\frac{dM}{dx}`$ is zero at the outer sonic point (location of which is already determined above). This is a general condition (vide discussion at the end of the previous Section).

(b) Mach No. at the outer sonic point must match with the analytical value obtained from the Sonic point condition (Eq. 12).

(c) Mach No. ($`M_{}`$) at the location of the middle sonic point must match with the analytical value derived from approximate analytical solution obtained using the energy equation (Eq. 2 and 4a).

Keeping in mind that an algebraic equation which is beyond quartic cannot be solved analytically \[see, Abramowitz and Stegun (1970)\], we expand $`M_\pm ^2`$ as quadratic equations so that Eq. (23) may become quartic. We shall show a posteriori that such an assumption introduces a very small and tolerable error in our computation.

If $`x_s`$ denotes the shock location, we assume,

$$M_\pm ^2=\underset{q=0}{\overset{2}{}}A_{[q,\pm ]}x_s^q$$
(23)
where, $`A_{[q,\pm ]}`$ are constant co-efficients to be determined from the conditions mentioned above. We find them to be,

$`A_{[2,+]}=\frac{1\left(M_+^2\right)_{mid}}{(x_{inf}x_{mid})^2}`$,

$`A_{[1,+]}=2x_{mid}A_{[2,+]}`$,

$`A_{[0,+]}=1+\left(2x_{mid}x_{inf}x_{inf}^2\right)A_{[2,+]}`$,

$`A_{[2,]}=\frac{\left(M_{}^2\right)_{mid}1}{(x_{out}x_{mid})^2}`$,

$`A_{[1,]}=2x_{mid}A_{[2,]}`$,

$`A_{[0,]}=1+x_{out}^2A_{[2,]}`$,

where, $`x_{mid}`$, $`x_{inf}`$ and $`x_{out}`$ are the middle (O-type) sonic point, the position where first derivative of Mach No. is infinity and the outer sonic point respectively.

We now substitute the above expression (Eq. 24) for the Mach number in the Mach invariant relation (Eq. 23) to obtain the following algebraic equation:

$`𝒜x_s^4+x_s^3+𝒞x_s^2+𝒟x_s+=0,`$ (24)
where,

$`𝒴=\left(3\gamma 1\right)\left(\gamma +1\right)`$,

$`𝒵=2(\gamma 1)`$,

$`𝒜=𝒴A_{[2,+]}A_{[2,]}`$,

$`=𝒴\left(A_{[1,+]}A_{[2,]}+A_{[2,+]}A_{[1,]}\right)`$,

$`𝒞=𝒴\left(A_{[0,+]}A_{[2,]}+A_{[1,+]}A_{[1,]}+A_{[2,+]}A_{[0,]}\right)𝒵\left(A_{[2,+]}+A_{[2,]}\right)`$,

$`𝒟=𝒴\left(A_{[0,+]}A_{[1,]}+A_{[1,+]}A_{[0,]}\right)𝒵\left(A_{[1,+]}+A_{[1,]}\right)`$,

$`=𝒴A_{[0,+]}A_{[0,]}𝒵\left(A_{[0,+]}+A_{[0,]}\right)4`$.

We solve this equation for $`X_s`$ analytically using the same procedure as in §4 (details in Appendix A). We denote the discriminant $`D`$ by $`D_s`$, and $`Q`$ and $`R`$ values as $`Q_s`$ and $`R_s`$ respectively for our discussion on the parameter space behavior of shocks.

In Figure 2, we redraw the parameter space as in Figure 1, but consider the formation of shocks alone. We find that $`Q_s>0`$ produces no shock from above, and $`Q_s=0`$ with $`R_s0`$ gives the boundary of the weakest shock (shocks with unit compression ratio). This boundary, though obtained using our approximate analytical method, generally coincides with the dotted curve of Figure 1. The edge of the boundary is obtained with an extra condition $`R_s=0`$. Thus, $`R_s`$ progressively decreases towards the edge along the dashed curve. This edge (denoted by $`D_s=R_s=0`$) ought to have coincided with the cusp of the $`D_c=0`$ curve drawn for the sonic point (see, Fig. 1 also), had the analytical method been exact. A small shift is the evidence that a small error is present at this corner of the parameter space. We also provide the region of the oscillating shocks ($`Q_s<0`$ and $`D_s<0`$). Here, the shock location is imaginary and therefore shock continuously oscillates back and forth causing a very interesting astrophysical effect known as Quasi-Periodic Oscillations (QPOs). This would be discussed in the next Section. The boundary between the shock and no-shock region from below is denoted by the dashed curve marked by $`D_s=0`$. Below the no-shock region where the energy and angular momentum of the flow are very low, the flow has only one sonic point.

In Figure 3, we draw shock location $`x_s`$ (along Y-axis) as a function of energy (along X-axis). Different curves are drawn for different specific angular momentum of the flow. The rightmost one is for $`\lambda =1.51`$ and the leftmost one is for $`\lambda =1.84`$, interval of $`\lambda `$ is $`0.01`$. As angular momentum is increased, the shock is located farther from the black hole. Comparing the locations from those obtained analytically, one notices that the same location is obtained for a specific angular momentum slightly more ($`3`$%) than that used in the numerical method. We therefore believe that the results obtained are very much reliable.

In Figure 4, we present a comparison of the boundary of the parameter space for which standing shocks may form as obtained by our analytical solution (shaded region) and by the numerical means (solid curve) existing in the literature (C89). The agreement is very good except in a region near the cusp (as also noted while discussing Fig. 2). Since very little parameter space is involved at this edge, we think that this small mismatch is tolerable.

## 6 Astrophysical Applications

Even though a black hole has no hard surface, it is remarkable that matter forms standing shocks around it in the same way as a shock is formed when a supersonic flow encounters a hard boundary. Shock waves heat up a gas and puff it up. This post-shock region intercepts soft photons from the pre-shock matter, particularly from the Keplerian disk which is located on the equatorial plane (CT95). In this scenario, the nature of the Comptonized radiation depends on the amount of matter in the sub-Keplerian and in the Keplerian flow: if the intensity of soft photon is very low, they cannot cool down the post-shock region by inverse Compton process and the spectrum remains very hot. On the other hand, if the intensity of soft photon is very high (i.e., the Keplerian rate is large), they cool down the post-shock region to the extent that the shock could not be sustained (pressure balance condition breaks down). This produces a soft state spectrum with a hard tail due to bulk motion Comptonization (CT95). There are several models in the literature which perhaps explain the soft/hard states. However, no model other than CT95 explains the power-law hard-tail in the soft state. Similarly, when the question of Quasi-Periodic Oscillations (QPOs) come, shock-oscillation model turns out to be a sufficiently satisfactory one (CM00).

When the parameters fall in the ‘no shock’ region of Fig. 2, shock location becomes imaginary. However, three sonic points are still present and the entropy of the flow at the inner sonic point continues to be higher compared to that at the outer sonic point. In this case, shock starts oscillating with a time period ($`T_s`$) comparable to the infall time from the post-shock region (Ryu, Chakrabarti and Molteni, 1997). Even when shocks form, if the infall time-scale turns out to be comparable to the cooling time, then the resonance condition is satisfied (Molteni, Sponholz and Chakrabarti, 1996; Chakrabarti and Manickam, 2000) and shocks oscillate in time scales of:

$$T_sx_s/v_s,$$

where, $`v_sR^1x_s^{1/2}`$ is the infall velocity and $`R`$ is the compression ratio at the shock (easily obtained analytically from our equations). Observed QPO frequencies are comparable to $`1/T_s`$.

When a mixture of Keplerian and sub-Keplerian matter is accreted, it is easy to obtain the parameters $`\lambda `$ and $``$ in terms of the Keplerian ($`\dot{M}_d`$) and sub-Keplerian ($`\dot{M}_h`$) accretion rates. Suppose, the viscosity parameter is such that the flow is deviated from a Keplerian disk at $`x=x_K`$ where its energy and angular momentum were $`_d`$ and $`\lambda _K`$ respectively. Further suppose that the sub-Keplerian halo has a constant energy $`0_h`$ and constant angular momentum $`\lambda _h`$, then the average angular momentum and energy of the transonic flow would be,

$$<>=\frac{\dot{M}_d_d+\dot{M}_h_h}{\dot{M}_d+\dot{M}_h}$$

$$<\lambda >=\frac{\dot{M}_d\lambda _d+\dot{M}_h\lambda _h}{\dot{M}_d+\dot{M}_h}$$

It is easy to compute the shock location of the resultant flow using our formalism given above.

It is to be noted from Fig. 3 that shock solutions are allowed only if the specific energy is positive. In other words, if a flow deviates from a cool Keplerian disk on the equatorial plane, the flow cannot have shocks as the specific energy in such a flow would be negative unless this flow is mixed with a substantial amount of sub-Keplerian matter with a positive energy. Typically, $`\dot{M}_h>>\dot{M}_d`$ and even with a small energy, the specific energy of the mixture becomes positive, giving rise to shocks and (unbound) winds. In the case when magnetic dissipation is present, flow energy could increase to a positive value and a solution with a shock would be allowed. The prospect of magnetic energy dissipation has been discussed by several workers in the literature (Shapiro, 1973; Bisnovatyi-Kogan and Blinnikov, 1976; Bisnovatyi-Kogan, 1998). Briefly, since the magnetic field rises as $`B_rr^2`$ and therefore magnetic pressure rises as $`P_{mag}r^4`$ while the gas pressure in the sub-Keplerian matter goes as $`P_{gas}r^{5/2}`$, magnetic field in excess of the equipartition value would escape from the disk buoyantly and may dissipate at the atmosphere as in the case of the Sun. If the flow has specific energy $`_h`$ at, say $`r=100r_g`$ where the flow was in equipartition, then at the shock, the energy would be at least $`2^4=16`$ times larger if all the magnetic energy is dissipated into the flow. Thus, a basically free-fall matter of $`10^4`$ would have an energy $`10^3`$ and thus a shock at a few tens of Schwarzschild radii would be expected.

An accreting flow can intercept hard X-rays emitted at the inner edge and this pre-heating effect need not be negligible. For instance, a flow emitting isotropically with $`6\%`$ efficiency will definitely intercept $`\theta \mathrm{\Theta }_{in}/4\pi `$ fraction of radiation in between the shocked region and the Keplerian disk. Assuming $`\theta 0.1`$, the energy deposition due to pre-heating is $`0.006`$ which is significant. This would energize Keplerian matter as well and shocks in the sub-Keplerian flow would be expected.

## 7 Concluding Remarks

So far in the astrophysical literature, existence of shocks in accretion flows has been indicated by steady and time-dependent numerical simulations. Study of these standing and oscillatory shocks in an accretion flows has been shown to be of great importance. Presently we show that the shocks could be studied completely analytically, at least in the case of thin, axisymmetric, inviscid flows with positive energy. We note that shock locations vary with flow parameters in a simple way — they form farther from a black hole when angular momentum is increased. This proves that they are mainly centrifugal pressure supported.

Given that the shocks, especially the standing shocks are the ideal locations at which a flow is heated up, hard X-rays are produced from the post-shock region after the soft-photons are processed by the flow by inverse-Comptonization process. Thus, the spectral states and the time-dependent behavior of the hard X-rays are directly related to the behavior of this region. For instance, CT95 computes steady state spectra with post-shock region as the source of hot, Comptonizing electrons. CM00 establishes that QPOs are due to oscillations of this region, since only hard X-rays are seen to exhibit QPOs. We therefore believe that shocks should be an important ingredient of an accreting system. However, if the disk is cool and Keplerian very far away, the specific energy must be negative. So the problem is not whether shocks should exist, the problem seems to be how to energize Keplerian matter as it becomes sub-Keplerian, by, say, magnetic energy dissipation, pre-heating etc. These work is in progress and would be reported elsewhere.

From the observers point of view, our work could also be useful, since the steady spectra, QPO frequencies, etc. are, in principle, determined analytically from a few free physical parameters. In future, we shall apply our understanding to obtain spectral properties and QPO behaviors more quantitatively.

This work is partly supported by a project (Grant No. SP/S2/K-14/98) funded by Department of Science and Technology (DST).

## Appendix A Method of the Analytical Solution

The procedure of obtaining an analytical solution of a quartic equation,

$$q^4+b_1q^3+b_2q^2+b_3q+b_4=0$$
$`A.1`$
is to first obtain a solution of the following cubic equation:

$$p^3+a_1p^2+a_2p+a_3=0,$$
$`A.2`$
where,
$`a_1=b_2`$, $`a_2=b_1b_34b_4`$, and $`a_3=4b_2b_4b_3^2b_1^2b_4`$ .

Let,

$`Q=\frac{3a_2a_1^2}{9}`$, $`R=\frac{9a_1a_227a_32a_1^3}{54}`$,

$`S=\sqrt[3]{R+\sqrt{Q^3+R^2}}`$, $`T=\sqrt[3]{R\sqrt{Q^3+R^2}}`$.

The discriminant is defined as

$$D=Q^3+R^2.$$
$`A.3`$
If $`D>0`$, one root is real and two roots are complex conjugate. In this case, the real solution is,

$$p_1=S+T\frac{1}{3}a_1.$$
$`A.4`$
If $`D=0`$, all roots are real and at least two are equal.

If $`D<0`$, all roots are real and unequal. They are:

$$p_1=2\sqrt{Q}\mathrm{cos}(\frac{1}{3}\theta )\frac{1}{3}a_1,$$
$`A.5`$
$$p_2=2\sqrt{Q}\mathrm{cos}(\frac{1}{3}\theta +120^o)\frac{1}{3}a_1,$$
$`A.6`$
and,

$$p_3=2\sqrt{Q}\mathrm{cos}(\frac{1}{3}\theta +240^o)\frac{1}{3}a_1,$$
$`A.7`$
where, $`\mathrm{cos}\theta =R/\sqrt{Q}`$.

One can now write a quadratic equation using any one of the real solutions of the cubic equation (see, Spiegel, 1968) as follows:

$$z^2+\frac{1}{2}\{b_1\pm \sqrt{b_1^24b_2+4p_1}\}z+\frac{1}{2}\{p_1\sqrt{p_1^24b_4}\}=0.$$
$`A.8`$
This is a quadratic equation which can be solved easily. Since we applied this procedure both of the sonic points and shocks, we denoted quantities like, $`D`$,$`Q`$ and $`R`$ by $`D_c`$,$`Q_c`$ and $`R_c`$ for sonic (critical) points and $`D_s`$,$`Q_s`$ and $`R_s`$ for shocks respectively.


text:
# The Phase Diagram of the Shastry-Sutherland Antiferromagnet

## I INTRODUCTION

The recent discovery of the two-dimensional spin gap system SrCu2(BO3)2 , and its representation as an antiferromagnetic $`S=\frac{1}{/}2`$ Heisenberg spin model equivalent to a model previously introduced by Shastry and Sutherland, has led to many recent studies of this model.

The Shastry-Sutherland model is a nearest-neighbor square lattice antiferromagnet, with additional diagonal interactions, in a staggered pattern, on alternate squares. The Hamiltonian is written as

$$H=J\underset{\mathrm{diag}}{}𝐒_i𝐒_j+J^{}\underset{\mathrm{axial}}{}𝐒_i𝐒_j,$$
(1)
This convention is adopted because in SrCu2(BO3)2 the diagonal bonds are, in fact, the shortest and the exchange constants are estimated to be $`J100K`$, $`J^{}68K`$. The model is illustrated in Figure 1.

It is clear that the system will exhibit Néel order for $`J/J^{}`$ small, and will form a gapped spin-singlet state for $`J/J^{}`$ large, with $`S=0`$ dimers on the diagonal bonds. Shastry and Sutherland showed that this dimer state is an exact eigenstate of the Hamiltonian for any $`J,J^{}`$ and that it is rigorously the ground state for $`J/J^{}>2`$.

The main controversial and challenging question is whether the system has an intermediate phase, and what kind of phase it is, if it exists. In the classical limit, when $`S\mathrm{}`$, the ground state is Néel ordered if $`J/J^{}1`$ and is helically ordered (see Fig. 1) otherwise, where the twist between one spin and its nearest neighbor is given by $`q=\mathrm{arccos}(J^{}/J)`$. For the $`S=\frac{1}{/}2`$ system, Albrecht and Mila used a Schwinger boson mean-field theory (SBMFT) and found an intermediate phase with helical long-range order (LRO) for $`1.1<J/J^{}<1.65`$, but with $`q`$ differing from its classical value. In our previous series study, we were unaware of this work, and did not consider the possibility of a helical LRO intermediate phase. We computed series for the Néel-ordered phase and for the dimer phase at $`T=0`$, and located a transition point at $`J/J^{}=1.45(1)`$ which appeared to be first-order.

Subsequently Koga and Kawakami calculated series expansions about disconnected plaquettes, and claimed to have identified two transition points: a second-order transition from Néel order to a plaquette singlet phase at $`J/J^{}1.16`$, followed by a first-order transition to the dimer phase at $`J/J^{}1.48`$. We find this interpretation unconvincing, partly because of uncertainty in the analysis of rather short series, and partly because such an intermediate phase is totally different from the helical-ordered phase presented in Ref. . Another recent development is the field-theoretical study of a generalized model with $`Sp(2N)`$ symmetry by Chung, Marston and Sachdev . For $`S=\frac{1}{/}2`$, they suggest a Néel to helical LRO phase transition at $`J/J^{}1.02`$, and a helical LRO to dimer short-range ordered (SRO) phase transition at $`J/J^{}2.7`$. This is similar to the scenario in Ref. , but with the helical phase extending over a much larger range. The quantitative accuracy of this approach for $`S=\frac{1}{/}2`$ is problematic. Both our previous work and the plaquette series of Koga and Kawakami provide strong evidence for the dimer phase setting in at around $`J/J^{}1.451.48`$. It is also interesting to note that for $`1/S>5`$ (S is a continuous variable in this theory) a phase with plaquette SRO is predicted. If finite $`N`$ fluctuations change the phase diagram significantly then it is conceivable, according to this theory, that there are four stable phases.

Knetter et al. subsequently performed a new series analysis of the dimer phase. They find that the lowest single-triplet excitation vanishes at $`J^{}/J=0.69`$ (or $`J/J^{}=1.45`$), very close to the transition point mentioned above; but that an $`S=1`$ two-triplet excitation energy vanishes even earlier, at $`J^{}/J=0.63`$ (or $`J/J^{}=1.59`$). This would seem to indicate important binding between the triplets, possibly leading to a condensate of triplets at or even below $`J^{}/J=0.63`$. Even more recently, Totsuka, Miyahara, and Ueda have discussed two-triplet binding effects in this model using both perturbation theory and exact diagonalization, and have shown that binding occurs even in the quintet channel of two triplets.

All the above studies show that the quantum system is Néel ordered up to and beyond $`J/J^{}=1`$, i.e. beyond the regime for the classical system: this is also consistent with the general phenomenon of quantum order by disorder. The transition to dimer order is located at $`J/J^{}1.451.58`$, so if there is an intermediate phase, it must exist within the region $`1<J/J^{}1.5`$. Our aim, in the present paper, is to explore these issues through new extended series expansions. We have obtained extended series expansions about isolated plaquettes, both with and without diagonal bonds, following Koga and Kawakami. Secondly, we have extended our previous series about the Néel ordered state and the dimer singlet state. Finally, we present new series expansions about states with helical LRO and columnar dimer order. The next Section presents technical details of the calculations, with analysis of series and discussion of results given in Section III. In the last section we present our conclusions.

We conclude that an intermediate phase between the Néel phase and the dimer phase very likely exists, in the region $`1.2J/J^{}1.5`$. The nature of this phase is much less certain, however. It does not appear to be plaquette ordered or helically ordered. It could perhaps posses weak columnar dimer order; or it could posses some other type of order, not considered here. In any case, this intermediate phase appears to be gapless, or nearly so.

## II Series Expansions

We use the linked-cluster expansion method to derive perturbation expansions for various choices of the unperturbed reference Hamiltonian. The technical details are discussed in many articles, and we refer the reader who is unfamiliar with these to a recent review. The following is a brief summary,

### A Dimer Expansions

The first term in the Hamiltonian (1), which represents disconnected dimers, is taken as the unperturbed Hamiltonian. The unperturbed ground state is then a product state of $`S=0`$ singlets on these dimers. As mentioned above, this state remains an exact eigenstate of the system for all $`J,J^{}`$, but is not the true ground state for $`J^{}/J>(J^{}/J)_c`$. The perturbation, being the second term in (1), mixes this state with states in which triplet excitations occur on the dimers.

Dimer expansions can be developed for both ground state properties and for excitations. Here because the dimer state is the exact ground state, we focus on the triplet excitation spectrum $`\mathrm{\Delta }(𝐤)`$. We have computed expansions to order $`(J^{}/J)^{21}`$, extending our previous calculation by 6 terms. This calculation involves 185332 linked clusters up to 11 sites. The resulting series coefficients are available on request. The minimum gap is at $`𝐤=0`$, where the series is:

$`\mathrm{\Delta }(𝐤=0)/J`$ $`=1x^2{\displaystyle \frac{x^3}{2}}{\displaystyle \frac{x^4}{8}}+{\displaystyle \frac{5x^5}{32}}{\displaystyle \frac{7x^6}{384}}{\displaystyle \frac{2051x^7}{4608}}{\displaystyle \frac{39091x^8}{55296}}`$ (5)
$`{\displaystyle \frac{268849x^9}{663552}}+{\displaystyle \frac{964411x^{10}}{6635520}}+{\displaystyle \frac{6597973x^{11}}{58982400}}{\displaystyle \frac{183919894867x^{12}}{191102976000}}2.180975296x^{13}`$
$`1.938901500x^{14}0.071485040x^{15}+0.820962395x^{16}2.236236507x^{17}`$
$`7.855506946x^{18}9.452785582x^{19}2.833667089x^{20}+4.387303805x^{21}+O(x^{22})`$
where $`x=J^{}/J`$. A standard Dlog Padé approximant analysis gives the gap vanishing at $`J^{}/J=0.6965(15)`$, or $`J/J^{}=1.436(3)`$ with the estimated critical index $`\nu 0.45(2)`$. Different Dlog Padé approximants show remarkable consistency with this exponent value, suggesting that this transition may belong to a new universality class. These results are in agreement with Müller-Hartmann et al., who used shorter series. But we recall that Knetter et al. found that a two-triplet $`S=1`$ bound-state engery vanishes even before $`J^{}/J=0.70`$, at around 0.63, which would entirely alter the position and the nature of this transition. We have not computed series for the two-particle states in this model. It would be very interesting to explore this transition further, using series or other methods.

We also computed a new series, to order $`(J^{}/J)^{17}`$, for the susceptibility $`\chi `$ corresponding to this momentum $`𝐤=0`$. The resulting series is

$`\chi (𝐤=0)`$ $`=1+{\displaystyle \frac{x^2}{2}}+{\displaystyle \frac{x^3}{8}}+{\displaystyle \frac{41x^4}{96}}+{\displaystyle \frac{281x^5}{1152}}+{\displaystyle \frac{6461x^6}{13824}}+{\displaystyle \frac{65953x^7}{165888}}+{\displaystyle \frac{186863x^8}{311040}}`$ (8)
$`+{\displaystyle \frac{374542669x^9}{597196800}}+{\displaystyle \frac{24573240371x^{10}}{28665446400}}+{\displaystyle \frac{2138459511091x^{11}}{2149908480000}}+1.308076648x^{12}`$
$`+1.604414343x^{13}+2.086275141x^{14}+2.638378225x^{15}+3.431256777x^{16}+4.419063264x^{17}+O(x^{18})`$
A standard Dlog Padé analysis of this series indicates a preferred critical point around $`J^{}/J=0.71`$. Biasing $`x_c`$ at 0.6965, the critical exponent is found to be very small, around 0.03, which is completely different from that of the classical 3D Heisenberg model. It might even be compatible with a logarithmic divergence.

### B Plaquette Expansions

Instead of perturbing about isolated dimers, any set of disconnected units can be used. Using a plaquette basis allows the investigation of plaquette type order. Because each plaquette has 16 rather than 4 states it is not possible to derive series of the same length as with dimer expansions.

Following Koga and Kawakami, we have computed two kinds of plaquette expansions (PE1 and PE2), which respectively take the set of plaquettes without the diagonal $`J`$ bonds and the set with diagonal $`J`$ bonds as the unperturbed Hamiltonian. To make the expansion possible it is necessary to introduce an expansion parameter $`\lambda 1`$ which modifies the interactions not included in $`H_0`$. This is illustrated in Figure 2. The series are then computed in powers of $`\lambda `$. The analysis evaluates these at $`\lambda =1`$, corresponding to the original Hamiltonian.

For PE1 we have computed series to order $`\lambda ^8`$, $`\lambda ^7`$, $`\lambda ^6`$ for the ground state energy, triplet excitation energies $`\mathrm{\Delta }(𝐤)`$ and staggered susceptibility $`\chi _{\mathrm{AF}}`$ respectively. This is one additional term for $`E_0`$ and two additional terms for $`\mathrm{\Delta }(𝐤)`$ and $`\chi _{\mathrm{AF}}`$ over Ref. . The calculation is computationally demanding, even with our efficient program. For example, the computation of $`\mathrm{\Delta }(𝐤)`$ to order $`\lambda ^7`$ took about 10 days and required 1.8GB memory on an SGI Origin 2400 system with a 400MHz R12000 CPU. The effort for each additional term requires factors of approximate 50 and 15 increase in CPU time and memory, respectively. For the second expansion (PE2) the series have been computed to order $`\lambda ^7`$ for the ground state energy and the singlet excitation spectrum, and to order $`\lambda ^6`$ for the triplet excitation $`\mathrm{\Delta }(𝐤)`$ and staggered susceptibility $`\chi _{\mathrm{AF}}`$. In Tables I and II we present the series for the ground state energy, the triplet gap $`\mathrm{\Delta }`$ at $`𝐤=0`$ and staggered susceptibility $`\chi _{AF}`$ for $`J/J^{}=1.25`$ and 1.4. Other series are available on request. The analysis is left for the next Section.

### C Columnar Dimer Expansions

A plausible candidate for the intermediate phase is the columnar dimer phase, as shown in Fig. 3. For both the spin-$`\frac{1}{2}`$ $`J_1J_2`$ Heisenberg model and the spin-$`\frac{1}{2}`$ Heisenberg model on an anisotropic triangular lattice, there appears to be a columnar dimer ordered phase adjacent to the Néel ordered phase. To make the expansion possible it is necessary to introduce an expansion parameter $`\lambda 1`$ which modifies the interactions not included in $`H_0`$. This is illustrated in Figure 3. The series are then computed in powers of $`\lambda `$. The analysis evaluates these at $`\lambda =1`$, corresponding to the original Hamiltonian.

We have computed series to order $`\lambda ^{10}`$, $`\lambda ^9`$, $`\lambda ^9`$ for the ground state energy, the antiferromagnetic susceptibility $`\chi _{\mathrm{AF}}`$, and the triplet excitation energies $`\mathrm{\Delta }(𝐤)`$ respectively. This calculation involves 290215 linked clusters up to 10 sites. In Table I we present the series for the ground state energy $`\chi _{AF}`$ for $`J/J^{}=1.25`$ and 1.4. Other series are available on request. The analysis is left for the next Section.

### D Ising Expansions for Néel order

For small $`J/J^{}`$, the model will have a Néel ordered phase with long range antiferromagnetic order (in the $`z`$ direction). To construct an Ising expansion, we introduce an exchange anisotropy parameter $`\lambda `$, and write the Hamiltonian as

$$H=H_0+\lambda V,$$
(9)
where

$`H_0`$ $`=`$ $`J{\displaystyle \underset{\mathrm{diag}}{}}S_i^zS_j^z+J^{}{\displaystyle \underset{\mathrm{axial}}{}}S_i^zS_j^z,`$ (10)
$`V`$ $`=J`$ $`{\displaystyle \underset{\mathrm{diag}}{}}(S_i^xS_j^x+S_i^yS_j^y)+J^{}{\displaystyle \underset{\mathrm{axial}}{}}(S_i^xS_j^x+S_i^yS_j^y),`$ (11)
The unperturbed ground state is the classical Néel state with energy $`E_0/N=J^{}/2+J/8`$. The perturbation $`V`$ flips pairs of neighbouring spins. Series are developed in powers of $`\lambda `$ and evaluated at $`\lambda =1`$, which recovers the original Hamiltonian. We have obtained series for the ground state energy per site, $`E_0/N`$, and the staggered magnetization (order parameter) $`M`$ to order $`\lambda ^{12}`$, extending our previous calculation by 3 terms. We have also computed a new series, to order $`\lambda ^{12}`$ for the perpendicular susceptibility $`\chi _{}`$. The resulting series for the ground state energy and the staggered magnetization (order parameter) $`M`$ for $`J/J^{}=1.25`$ and 1.4 are listed in Tables I and III; the series for other values of $`J/J^{}`$ are available on request. For details of the analysis we refer to our previous paper. Results of the analysis are presented in the next section.

### E Ising Expansions for Helical Order

As discussed by Albrecht and Mila, the classical system has planar helical order for $`J/J^{}>1`$. Starting from a reference spin in the $`z`$ direction, each neighboring spin is rotated by an angle $`q`$, as shown in Figure 1. The twist is determined by minimization of the energy and yields

$$q=\{\begin{array}{cc}\pi ,\hfill & JJ^{}\text{;}\hfill \\ \mathrm{arccos}(J^{}/J),\hfill & J>J^{}\text{.}\hfill \end{array}$$
(12)
To develop an Ising expansion about such a helically ordered state we transform the Hamiltonian by rotating the spin axes at each site. The transformed Hamiltonian is

$$H=H_1+\lambda (H_2+H_3),$$
(13)
where

$`H_1`$ $`=`$ $`J\mathrm{cos}(2q){\displaystyle \underset{\mathrm{diag}}{}}S_i^zS_n^z+J^{}\mathrm{cos}(q){\displaystyle \underset{\mathrm{axial}}{}}S_i^zS_j^z,`$ (14)
$`H_2`$ $`=`$ $`J{\displaystyle \underset{\mathrm{diag}}{}}\left[S_i^yS_n^y+\mathrm{cos}(2q)S_i^xS_n^x+\mathrm{sin}(2q)(S_i^zS_n^xS_i^xS_n^z)\right],`$ (15)
$`H_3`$ $`=`$ $`J^{}{\displaystyle \underset{\mathrm{axial}}{}}\left[S_i^yS_j^y+\mathrm{cos}(q)S_i^xS_j^x+\mathrm{sin}(q)(S_i^zS_j^xS_i^xS_j^z)\right].`$ (16)
and where exchange anisotropy is introduced through the perturbation parameter $`\lambda `$.

We have computed series for the ground state energy and the order parameter to order $`\lambda ^{12}`$, for various choices of $`q`$ and the ratio $`J/J^{}`$. The resulting series for the ground state energy and the staggered magnetization (order parameter) $`M`$ for $`J/J^{}=1.25`$ and 1.4 are listed in Tables I and III; other series can be supplied on request. We are not aware of previous series of this type for this model.

## III Results and Discussion

In this Section we present a variety of results obtained from analysis of the various series expansions. The analysis has been carried out using integrated first-order inhomogeneous differential approximants and Padé approximants to extrapolate each series to the physical value $`\lambda =1`$.

### A Ground State Energy

The nature of the ground state for any particular $`J/J^{}`$ can be determined by comparing the energies of various candidate states. Our analysis shows that, among all the series expansions we have calculated, the ground state energy from the Ising expansion about Néel order has the lowest energy for coupling $`J/J^{}<1.45`$, while the ground state energies obtained from the columnar dimer expansion and from the Néel order Ising expansion remain nearly equal through $`J/J^{}1.2`$. For example, we estimate the ground state energy at $`J/J^{}=1.25`$ as

$$E_0/NJ^{}=\{\begin{array}{cc}0.5628(15),\hfill & \text{Néel order;}\hfill \\ 0.5586(10),\hfill & \text{helical order;}\hfill \\ 0.5599(3),\hfill & \text{PE1;}\hfill \\ 0.5584(8),\hfill & \text{PE2;}\hfill \\ 0.563(1),\hfill & \text{columnar dimer.}\hfill \end{array}$$
(17)
and at $`J/J^{}=1.40`$ as

$$E_0/NJ^{}=\{\begin{array}{cc}0.561(3),\hfill & \text{Néel order;}\hfill \\ 0.554(1),\hfill & \text{helical order;}\hfill \\ 0.5548(5),\hfill & \text{PE1;}\hfill \\ 0.5525(10),\hfill & \text{PE2;}\hfill \\ 0.557(2),\hfill & \text{columnar dimer.}\hfill \end{array}$$
(18)
where we have used the classical-$`q`$ value in Eq. (12) in the calculation of helical order. In Table IV we show details of the analysis from the latter case $`J/J^{}=1.4`$ for both the Ising expansion about Néel order and the plaquette expansion (PE1), from which we estimate the values in (18). The ground state energies obtained from various expansions versus $`J/J^{}`$ are given in Fig. 4. Since the energy from the second plaquette expansion is slightly higher than that from the first plaquette expansion, we do not show the results of the second plaquette expansion. Our results are thus in disagreement with ref. , which claimed that the plaquette phase has the lowest energy. Of course these energies are all very close to each other and the estimated errors are subjective confidence limits. However inspection of the data in Table IV shows that the Ising expansion about Néel order consistently gives lower energy estimates.

Our results also argue against helical order being a stable ground state. In Figure 5 we show estimates of the ground state energy from Ising expansions about helical order, as a function of angle $`q`$, for the two coupling ratios $`J/J^{}=1.25`$, 1.40. In each case the curve shows no indication of a minimum at some $`q_c`$, but rather decreases monotonically towards $`q=\pi `$, corresponding to Néel order.

Over the region $`1.2J/J^{}1.45`$, the ground state energies obtained from the Néel ordered Ising expansion and from the columnar dimer expansion have some overlap after considering the error bars, so they are both good candidates for the ground state of the Shastry-Sutherland model.

### B Staggered Magnetization and Perpendicular Susceptibility

The staggered magnetization and perpendicular susceptibility will be non-zero in a phase with long range antiferromagnetic order, and are expected to vanish at a transition point to a magnetically disordered or spin-liquid phase. Effective Lagrangian theory predicts a relationship $`\rho _s=v^2\chi _{}`$ where $`\rho _s`$ is the spin stiffness and $`v`$ the spin-wave velocity. In Figures 6 we show estimates of $`M`$ and $`\chi _{}`$ versus $`J/J^{}`$, from the Ising expansion about Néel order. It can be seen that both quantities behave similarly, decreasing from their values at $`J=0`$ and vanishing at around $`J/J^{}=1.2\pm 0.1`$. The error bars are large and the vanishing point can not be obtained to high precision. However our new extended series shows a much sharper drop-off in $`M`$ than given in our previous work. Thus it is possible that the true transition point is at or below 1.20.

For the Ising expansion about helical order the magnetization $`M`$ at $`\lambda =1`$ (for the classical $`q`$ value) is zero, within error bars, for all $`J/J^{}`$. This strengthens our conclusion, from the ground state energy results, that a helical phase is not present.

### C Energy gap and susceptibility in the plaquette expansion

The main argument of Koga and Kawakami for the existence of a plaquette intermediate phase in the system is based on analysis of the triplet gap at $`𝐤=0`$ and the staggered susceptibility series. In their series analysis, they assumed that the minimum triplet gap is at $`𝐤=0`$, and the transition lies in the same universality class as the classical 3D Heisenberg model (i.e. critical exponents $`\gamma 1.4`$, $`\nu 0.71`$). Using standard Dlog Padé approximants, they found an apparent critical singularity at $`\lambda >1`$ for $`J/J^{}1.16`$. This would imply a non-zero spin gap for $`J/J^{}1.16`$ for the original Hamiltonian (1). The above assumption is valid for the transition to Néel order, but if there was an intermediate phase, the transition from the plaquette phase to this intermediate phase might not lie in the same universality class as the classical 3D Heisenberg model, and the minimum triplet gap might not be at $`𝐤=0`$. In the analysis with our longer series, we find indeed that the minimum triplet gap is not at $`𝐤=0`$ for $`J/J^{}1.25`$, and also that the transition does not lie in the classical 3D Heisenberg universality class.

Firstly, let us discuss the first plaquette expansion (PE1). Here for $`J/J^{}1.25`$, or for larger values of $`J/J^{}`$ with small values of $`\lambda `$, we find the minimum triplet gap is at $`𝐤=0`$; but for $`J/J^{}1.25`$ and $`\lambda 1`$, it is no longer at $`𝐤=0`$. For example for $`J/J^{}=1.40`$, the dispersion for various $`\lambda `$ along the line connecting $`𝐤=(0,0)`$ and $`𝐤=(\pi /4,3\pi /16)`$ is shown in Fig. 7, where we can see that for $`\lambda 0.75`$, the minimum gap is located at $`𝐤=0`$, while for larger $`\lambda `$, the minimum gap is located at about $`𝐤=(\pi /4,3\pi /16)`$, and the dispersion near the minimum is quite flat. The series for the triplet gap at $`𝐤=(\pi /4,3\pi /16)`$ is

$`\mathrm{\Delta }(\pi /4,3\pi /16)/J^{}`$ $`=10.3077152787\lambda 0.2293662299\lambda ^20.0346064996\lambda ^3`$ (20)
$`0.0741227651\lambda ^40.0176489659\lambda ^50.0400694385\lambda ^60.0030161628\lambda ^7+O(\lambda ^8)`$
Table V shows estimates of the critical point $`\lambda _c`$ and critical index from both the spin gap and susceptibility series for the coupling ratios $`J/J^{}=1.25`$ and 1.40, using unbiased Dlog Padé approximants, where for $`J/J^{}=1.40`$, we show the results for the spin gap at both $`𝐤=0`$ and $`𝐤=(\pi /4,3\pi /16)`$. Apart from the spin gap $`\mathrm{\Delta }(\pi /4,3\pi /16)`$ at $`J/J^{}=1.40`$, there is considerable scatter in the results, but several features are apparent. The spin gap series show decreasing estimates of $`\lambda _c`$ with increasing order for both $`J/J^{}`$ ratios. The data are more consistent with the conclusion $`\lambda _c1`$, rather than $`\lambda _c>1`$. For the spin gap $`\mathrm{\Delta }(\pi /4,3\pi /16)`$ at $`J/J^{}=1.40`$, the unbiased Dlog Padé approximants give very convergent results, with the critical point $`\lambda _c=0.988`$ and critical index $`\nu =0.43`$. This indicates instability of the plaquette phase at $`\lambda =1`$, at least for this coupling ratio. This also indicates that there is indeed an intermediate phase. Note that the critical index $`\nu `$ obtained here is about the same as that obtained from the dimer series, suggesting that this transition might belong to the same new universality class. Our analysis shows that the transition probably does not belong in the classical 3D Heisenberg universality class for $`J/J^{}1.1`$. Because the critical index $`\nu `$ for the minimum gap is much smaller than 1, we have obtained the dispersion results shown in Fig. 7 by performing the series extrapolation in a new variable

$$\delta =1(1\lambda /\lambda _c)^\nu .$$
(21)
in order to take the singular behavior into account.

Finally we consider the second plaquette expansion (PE2). The minimum triplet gap is no longer at $`𝐤=0`$, at least for small to moderate $`\lambda `$. This can be seen from the first few terms of the series

$$\mathrm{\Delta }(k_x,k_y)/J^{}=1\frac{4\lambda }{6}(\mathrm{cos}k_x+\mathrm{cos}k_y)+\frac{J\lambda }{3J^{}}\mathrm{cos}(k_x+k_y)$$
(22)
which gives the minimum gap at $`k_x=k_y=\mathrm{arccos}(J^{}/J)`$. It seems likely that the minimum triplet gap will remain at $`𝐤0`$, even for $`\lambda 1`$, although the series in this region are too erratic to confirm this. For this choice of unperturbed Hamiltonian we are also able to compute the excitation spectrum for the singlet excitation and we have obtained series to order $`\lambda ^7`$. The minimum singlet gap is at $`𝐤=0`$ and appears to be smaller than the triplet gap. An attempt to locate the critical point $`\lambda _c`$ by Dlog Padé approximants to this series was hampered by poor convergence. We have attempted other analysis procedures, but without great success. But since the ground state energy from the second plaquette expansion (PE2) is higher than that obtained from both the first plaquette expansion and the Néel ordered Ising expansion, the probability of having the second plaquette configuration as the intermediate phase is remote.

### D The energy gap in the columnar dimer expansion

Finally we discuss the triplet dispersion obtained from the columnar dimer expansion. Here for $`J/J^{}1.15`$, we find the minimum gap is located at momentum $`𝐤=(0,\pi )`$, as expected, since we expect to have a transition to the Néel ordered phase for small $`J/J^{}`$. Fig. 8 show the gap at $`𝐤=(0,0)`$, $`(\pi ,0)`$ and $`(0,\pi )`$ for $`J/J^{}=0.875`$ versus $`\lambda `$ obtained from the integrated differential approximants to the series. We can see that the minimum gap vanishes at a critical value $`\lambda _c<1`$. We also expect this transition (to Neél order) to lie in the same universality class as the classical 3D Heisenberg model. To determine the phase boundary of columnar dimer order, we can use Dlog Padé approximants to both the series for the minimum gap and for the antiferromagnetic susceptibility. To get more reliable estimates for the critical point, we assume the critical exponents to be $`\gamma 1.4`$, $`\nu 0.71`$. The results are shown in the phase diagram Fig. 9, where we can see that $`\lambda _c=0.514(6)`$ for $`J=0`$, and that $`\lambda _c`$ increases for increasing $`J`$. For $`J/J^{}=1.15`$, we estimate $`\lambda _c=0.97(5)`$, but unfortunately the biased Dlog Padé results here are not accurate enough to tell whether $`\lambda _c`$ is larger than 1.

For $`J/J^{}1.75`$, we find that the minimum gap is located at $`𝐤=(0,0)`$. Fig. 8 also shows the gap at $`𝐤=(0,0)`$, $`(\pi ,0)`$ and $`(0,\pi )`$ for $`J/J^{}=2`$ versus $`\lambda `$. Once again, one can locate the phase boundary by using the Dlog Padé approximants to the series for the minimum gap, and the results are also shown in Fig. 9. We see that as $`J/J^{}`$ decreases, $`\lambda _c`$ increases, but again, the analysis is not accurate enough to tell whether $`\lambda _c`$ is larger than 1 in the intermediate region.

In the analysis, we also note that the gap at $`𝐤=(\pi ,0)`$ shows some peculiar features, as shown in Fig. 8. For all values of $`J/J^{}`$, the majority of integrated differential approximants to the series show that the gap increases for small $`\lambda `$, then decreases dramatically to zero for $`\lambda 0.5`$. The majority of Dlog Padé approximants to the series also show that the series has a critical point at $`\lambda 0.5`$ with very small critical index ($`0.2`$). But we believe this peculiar behavious is an artefact of the short series.

For the most interesting region $`1.15J/J^{}1.75`$, the situation is quite complicated: the location of the minimum gap depends on both $`J/J^{}`$ and $`\lambda `$. We take the mid-point for the presumed intermediate phase, $`J/J^{}=4/3`$, as example. For $`\lambda `$ near 1, we find the minimum gap for $`J/J^{}=4/3`$ is at $`𝐤=(0.47,2.80)`$, slightly away from $`(0,\pi )`$. The standard Dlog Padé approximants to the series for the minimum gap gives a critical point $`\lambda _c=1.0(1)`$, with critical index $`\nu 0.4`$. This seems to be consistent with the indices obtained from the dimer and plaquette expansions, Since the standard Dlog Padé approximants cannot tell whether $`\lambda _c`$ is larger than 1, we examine the extrapolations for the gap itself obtained from integrated differential approximants. The results for the gap at $`𝐤=(0.47,2.80)`$, $`(0,0)`$ and $`(0,\pi )`$ are shown in Fig. 10. The gap at $`\lambda =1`$ for $`𝐤=(0.47,2.80)`$ is very small ($`0.1J^{}`$) when estimated from the direct $`\lambda `$-series. When analysed in term of variable $`\delta `$ (Eq. 21) the gap is compatible with zero. If the gap vanishes at or before $`\lambda =1`$, the columnar dimer phase would also be ruled out as intermediate phase for the Shastry-Sutherland model.

## IV Conclusions

We have attempted to further elucidate the nature of the phase diagram of the Shastry-Sutherland spin model at $`T=0`$, by series expansion methods. We have significantly extended previous series, computed by ourselves and others, and we have also derived a number of new series. The analysis of the various series allows us to draw some fairly firm conclusions, as well as others which are more tentative.

The first question is whether an intermediate phase between the Néel-ordered phase and the dimer phase exists at all. It seems rather clear that the singlet dimer phase, with a simple singlet-product ground state, persists from large $`J/J^{}`$ down to $`J/J^{}1.51.6`$. In our previous work we argued that there is a direct first-order transition from the dimer phase to Néel order at this point. Our present results show the staggered magnetization and perpendicular susceptibility in the Néel phase vanishing at $`J/J^{}1.2`$, with large uncertainty $`\pm 0.1`$. This appears more consistent with the Néel phase terminating at a second-order phase transition, so that there may indeed be an intermediate phase stretching over a very small range of coupling constants $`1.20J/J^{}1.5`$, in agreement with the suggestion of Koga and Kawakami.

Carpentier and Balents have recently discussed an effective mean-field theory approach to the generalized Shastry-Sutherland model. They argue that a direct continuous transition from the Néel phase to the dimer phase is not possible. If this argument is accepted, it means that the vanishing of the order parameter and susceptibility in the Néel phase cannot be associated with a second order transition directly to the dimer phase. It follows that an intermediate phase must exist, over at least a small range of couplings.

The next question is the nature of the transitions from the Néel and dimer phases, respectively, into the intermediate phase. The vanishing staggered magnetization and perpendicular susceptibility in the Néel phase indicate a 2nd order transition; but we are unable to estimate the critical exponents, since we do not have explicit series coefficients for an expansion in $`J/J^{}`$. The vanishing energy gap would indicate that the transition from the dimer phase is also 2nd order; but the result of Knetter et al. showing that the two-particle bound-state energy vanishes before the single triplet gap indicates that more work needs to be done to characterize this transition also. It is interesting that one-particle gaps appears to vanish with a similar exponent in the dimer expansion, the plaquette expansion, and also the columnar dimer expansion.

A further question concerns the nature of the intermediate phase. None of the suggested phases exhibit a ground-state energy which is distinctly lower than the Néel energy in the relevant coupling regime. Our results do not support the suggestion that the intermediate phase is plaquette ordered. The extrapolated ground-state energy from the Néel expansion consistently lies below that from the plaquette expansion; and the triplet gap from the plaquette expansion appears to vanish at or before the physical value $`\lambda =1`$ is reached, indicating instability in this phase. Our longer series thus appear to contradict the conclusions of Koga and Kawakami. Equally, our results do not support the suggestion that the intermediate phase is helically ordered. The extrapolated ground-state energy from the Néel expansion consistently lies below that of the helically-ordered state, for any value of the spin orientation $`q`$ other than $`\pi `$. We note that the mean-field approach of Albrecht and Mila is strictly only valid for $`d>2`$; while the $`1/N`$ expansion approach of Chung, Marston, and Sachdev is of doubtful validity for the spin-$`\frac{1}{/}2`$ case, and is certainly not quantitatively accurate.

We have also explored the possibility that the intermediate phase is a columnar dimer phase, as in the $`J_1J_2`$ Heisenberg model, or the spin-$`\frac{1}{2}`$ Heisenberg model on an anisotropic triangular lattice. The extrapolated ground state energy for this state is comparable with that of the Néel state, within error bars, and the triplet gap from this expansion shows a window $`1.15J/J^{}1.5`$ in which it may remain finite, though very small, at the physical value $`\lambda =1`$. This remains a possible candidate for the intermediate phase, although our data are not sufficiently accurate to allow a definitive conclusion.

The nature of an intermediate state in this model thus remains an open question. It could, for example, be a structureless spin liquid, in which case a signature of this phase is elusive. Carpentier and Balents have also discussed possible intermediate phases in the generalized model, including a WISDQ (weakly incommensurate spin-density wave) ordered state -i.e. a periodic modulation of the expectation value of the total spin on each dimer - and a ‘fractionalized’ state, with topological order and deconfined spin-$`\frac{1}{/}2`$ excitations (‘spinons’). We currently have no information to offer regarding these states. The fact that the energy gap drops to zero before or near the physical value $`\lambda =1`$ in every expansion we have tried appears to suggest that the intermediate phase is gapless, or near to it. Further support for this idea is given by the strong finite-size dependence of the exact diagonalization results of Miyahara and Ueda in this region.

The material SrCu2(BO3)2 appears to lie within the dimer phase of the Shastry-Sutherland model. Original estimates $`J/J^{}=1.47`$ put the material close to the phase transition point, and it was argued that some of the unusual properties of the material were due to this closeness. A more recent estimate, including interactions in the third direction, gives $`J/J^{}1.66`$, somewhat beyond even the point at which the two-particle gap vanishes. This puts the material clearly within the dimer phase and hence the question of intermediate phases is primarily of theoretical interest.

It is clear that much work remains to be done on this fascinating model. It would be interesting to confirm and extend the results of Knetter et al. concerning the transition from the dimer phase; and also to try and construct series coefficients for an expansion in $`J/J^{}`$ in the Néel phase. New numerical methods are badly needed to explore the intermediate phase also: the best nmethod to employ here remains a puzzle.

###### Acknowledgements.

This work forms part of a research project supported by a grant from the Australian Research Council. The computation has been performed on an SGI origin 2400 computer. We are grateful for the computing resources provided by New South Wales Australian Centre for Advanced Computing and Communication (ac3) center.