text:
# STABILITY OF STANDING ACCRETION SHOCKS, WITH AN EYE TOWARD CORE COLLAPSE SUPERNOVAE

## 1 INTRODUCTION

Accretion shocks occur in many settings in astrophysics, including core collapse supernovae, star formation, and accreting compact objects. In many cases these accretion shocks may exist in a quasi-steady state, with the postshock flow gradually settling onto the accreting object such that the radius of the accretion shock remains constant: a standing accretion shock, or SAS.

A particularly interesting example of a standing accretion shock arises in core-collapse supernovae. In classical models of such supernovae, an expanding shock front stalls at a radius $`100200`$ km, and remains relatively stationary for an eternity $`300`$ ms (see, for example, Mezzacappa et al., 2001, and references therein), during which time it is “revived” by an as yet unknown combination of factors including neutrino heating, convection, rotation, and magnetic fields.

Despite the ubiquity of standing accretion shocks in astrophysical lore, there is little discussion in the literature of their stability. This oversight may be attributed to the fact that strong, adiabatic shocks are known to be self-healing in that small perturbations to a planar shock front lead to postshock flow that acts against the perturbations (e.g., Whitham, 1974). The situation may be different in the case of a spherical accretion shock. The convergence of the postshock flow may amplify perturbations, and the confined postshock volume increases the opportunity for the postshock flow to feed back on the dynamics of the shock front. For example, an analysis of a converging blastwave suggests that such shocks are unstable for high Mach numbers (e.g., Whitham, 1974), in contrast with planar shocks. More recently, Foglizzo (2002) has described a vortical-acoustic instability of spherical accretion shocks in the context of accreting black holes. In this case the instability results from a feedback between acoustic waves perturbing the spherical shock and the perturbed shock exciting vorticity perturbations that advect downstream towards the accreting object. These vorticity waves generate acoustic waves that propagate back to the shock, completing the feedback cycle.

In this paper we investigate the stability of an SAS and, in particular, the role of oblique (aspherical) shocks in feeding turbulence in the confined postshock region. To this end we present one- and two-dimensional time-dependent hydrodynamics simulations of an idealized SAS. We begin by defining an idealized SAS and presenting an analytic solution for the postshock flow in Section 2. In Section 3 we describe the numerical method we use to evolve a time-dependent hydrodynamics simulation of an SAS, and in Section 4 we use these simulations to demonstrate the stability of an SAS in one dimension. In Section 5 we describe the unstable nature of an SAS in two dimensions. We conclude by discussing the relevance of an unstable accretion shock to the problem of core-collapse supernovae.

## 2 IDEALIZED MODEL OF A STANDING ACCRETION SHOCK

To model standing accretion shocks, we consider an idealized adiabatic gas in one dimension accreting onto a star of mass $`M`$. We begin by assuming that the infalling gas has had time to accelerate to free fall and that this velocity is highly supersonic. The fluid variables (velocity, $`u`$; density, $`\rho `$; pressure, $`p`$) just behind the standing accretion shock are then given by the Rankine-Hugoniot equations:

$$u=\frac{\gamma 1}{\gamma +1}(2GM)^{1/2}R_s^{1/2}$$
(1)
$$\rho =\frac{\gamma +1}{\gamma 1}\left(\frac{\dot{M}}{4\pi }\right)(2GM)^{1/2}R_s^{3/2}$$
(2)
$$p=\frac{2}{\gamma +1}\left(\frac{\dot{M}}{4\pi }\right)(2GM)^{1/2}R_s^{5/2}$$
(3)
where $`\gamma `$ is the adiabatic index and $`R_s`$ is the shock radius.

Below the shock front we assume radiative losses are negligible (we will discuss the appropriateness of this assumption, particularly in the case of core collapse supernovae, in more detail later) and the gas is isentropic, such that $`p/\rho ^\gamma `$ is constant. The assumption of a steady-state isentropic flow implies there is a zero gradient in the entropy of the postshock gas, which in turn implies this flow is marginally stable to convection. This allows us to separate effects due to convection from other aspects of the multidimensional fluid flow. In the more general case where a negative entropy gradient might drive thermal convection (Herant et al., 1994; Burrows, Hayes, & Fryxell, 1995; Janka & Mueller, 1996; Mezzacappa et al., 1998; Fryer & Heger, 2000), such convection could act as a seed for the instabilities we will examine in this paper.

Under these assumptions, the postshock flow is described by the Bernoulli equation:

$$\frac{1}{2}u^2+\frac{\gamma }{\gamma 1}\frac{p}{\rho }\frac{GM}{r}=0,$$
(4)
where the constant on the right hand side is set by the assumption of highly supersonic freefall above the accretion shock ($`P/\rho \frac{1}{2}u^2GM/r`$).

We normalize the problem such that $`GM=0.5`$, $`\dot{M}=4\pi `$, and the accretion shock is at $`R_s=1`$. Rewriting the Bernoulli equation in terms of velocity and radius, we have

$$ru^2+\alpha r^{32\gamma }u^{1\gamma }1=0,$$
(5)
where

$$\alpha =\frac{4\gamma }{(\gamma +1)(\gamma 1)}\left(\frac{\gamma 1}{\gamma +1}\right)^\gamma .$$
(6)
For $`\gamma =5/3`$, one can write down an analytic expression for the post-shock gas:

$$u(r)=\frac{\gamma 1}{\gamma +1}r^{1/2}.$$
(7)
In this special case, the Mach number of the postshock flow remains fixed at $`1/\sqrt{5}`$, as the gas accelerates inwards at a fixed fraction of the freefall velocity.

The behavior of the postshock gas is markedly different for smaller values of $`\gamma `$. Below a critical value, $`\gamma _c=1.522`$, the gas decelerates, rather than accelerates, as it drops below the standing accretion shock. This isentropic settling solution below the accretion shock is shown in Figure 1 for an adiabatic index $`\gamma =4/3`$.

The behavior of these settling solutions at small radii can be found by neglecting the first term in eqn. 5, which corresponds to the kinetic energy of the gas. This amounts to assuming that the local Mach number has become very small. The resulting solution represents an approximately hydrostatic atmosphere with a constant mass inflow rate. The velocity is then a power-law in radius:

$$u(r)r^{\frac{32\gamma }{\gamma 1}}.$$
(8)
For large values of $`\gamma `$ (but still below the critical value of 1.522) this is a relatively shallow power law, and the infalling gas slows down only when very close to $`r=0`$. For soft equations of state, however, the gas decelerates relatively quickly. To maintain the same pressure gradient required for (nearly) hydrostatic equilibrium, a gas with a softer equation of state must be compressed faster. This translates into a faster deceleration under the assumption of a constant mass inflow rate. For $`\gamma =4/3`$, the infall velocity drops linearly with radius; at $`r=0.1`$, the gas is settling inwards with a velocity of only 1% the local sound speed. These settling solutions are shown in Figure 2 for various values of $`\gamma `$. We note that these solutions are the same as those used by Chevalier (1989) to model the post-supernova fall-back of stellar gas onto the neutron star.

For comparison, in Figure 3 we plot our model SAS solution along with the results of a complete one-dimensional model of stellar core collapse, bounce, and the postbounce evolution, which includes Boltzmann neutrino transport and a realistic equation of state. These results show the structure of the stalled supernova shock at roughly 300 ms after core bounce, as described in Mezzacappa et al. (2001). Our SAS solution has been scaled to match the postshock density of the supernova model. For this particular scaling, the unit of time in our normalized SAS model corresponds to 6 ms in the supernova model.

## 3 NUMERICAL MODEL

We use the time-dependent hydrodynamics code VH-1 (http://wonka.physics.ncsu.edu/pub/VH-1/) to study the dynamics of a standing accretion shock. For this problem we have modified the energy conservation equation in VH-1 to evolve a total energy that includes the time-independent gravitational potential:

$$E_{tot}=\frac{1}{2}u^2+\frac{1}{\gamma 1}\frac{p}{\rho }+\mathrm{\Phi }$$
(9)
where $`\mathrm{\Phi }=GM/r`$ for the assumed gravitational field of a point mass. This approach ensures conservation of the total energy to within machine round-off error.

The numerical simulations are initialized by mapping the analytic solution onto a grid of 300 radial zones extending from $`0.1R_s`$ to $`2R_s`$. The radial width of the zones increases linearly with radius such that the local zone shape is roughly square in the two-dimensional simulations; i.e., $`\mathrm{\Delta }RR\mathrm{\Delta }\theta `$, to minimize truncation error. The two dimensional simulations used 300 angular zones covering a polar angle from $`0`$ to $`\pi `$, assuming axisymmetry about the polar axis. The fluid variables at the outer boundary are held fixed at values appropriate for highly supersonic free-fall at a constant mass accretion rate, consistent with the analytic standing accretion shock model.

For the purpose of our numerical simulations we must impose an inner boundary, which one would nominally associate with the surface of the accreting object. This boundary must allow a mass and energy flux off the numerical grid in order to maintain a steady state. We do this by imposing a leaky boundary condition: an inner boundary with a specified radial inflow velocity corresponding to the chosen radius of the inner boundary.

Our numerical code employs a Lagrange-remap algorithm, in which the Lagrangian fluid equations are evolved using solutions to the Riemann problem at zone boundaries. We fix the velocity at the inner boundary by replacing the Riemann solution for the time-averaged velocity at the inner boundary, with a fixed value. This value of the boundary flow velocity is chosen based on the analytic value at the boundary radius (see Figure 2), and, in effect, this value sets the height of the standing shock. The pressure gradient at the inner boundary is set to match the local gravitational acceleration, mimicking a hydrostatic atmosphere (the SAS is nearly hydrostatic far below the shock). The density at this leaky boundary is adjusted to force a zero gradient, allowing a variable mass flux off the inner edge of the numerical grid. This specific set of boundary conditions was influenced by the desire to construct time-dependent solutions that are stable in one dimension (see Section 4).

To test the dependence of our results on this inner radial boundary condition we also ran simulations in which we replaced the leaky boundary with a hard reflecting surface and added a cooling layer above it. This cooling layer is meant to mimic the energy loss due to neutrino cooling above the surface of the proto-neutron star in a more realistic model. Stationary accretion shock solutions with a cooling layer are described in Houck & Chevalier (1992). We have chosen a cooling function described by $`\alpha =\beta =2.5`$ in their notation. These values correspond to a radiation-dominated gas ($`pT^4`$) loosing energy through neutrino cooling with a negilible density dependence but a strong temperature dependence ($`\dot{E}T^{10}`$). We adjusted the magnitude of the cooling function to place the standing shock at a radius of unity.

For the two-dimensional simulations, reflecting boundary conditions were applied at the polar boundaries and, the radial boundaries were treated the same as in the one-dimensional simulations. The tangential velocity was initialized to zero everywhere in the computational domain. The tangential velocity at the inner radial boundary was set to zero throughout the simulations.

The PPM algorithm (Colella & Woodward, 1984) used in VH-1 is a good method to study the stability of accretion shocks because shock transitions are typically confined to only two numerical zones. However, a drawback of this algorithm is the generation of postshock noise when a strong shock aligned with the numerical grid moves very slowly across the grid. This is exactly the situation encountered in this problem. In one dimensional simulations this produces small entropy waves that advect downstream from the shock with little effect on the flow. In multiple dimensions, the zone-to-zone fluctuations in pressure and density transverse to the shock front generate considerable postshock noise (Leveque, 1998).

To minimize the effects of this numerical noise we have added a small amount of dissipation by “wiggling” the computational grid in both the radial and azimuthal directions (Colella & Woodward, 1984). The radial grid “wiggling” is confined to zones near the location of the steady-state shock so that our inflow and outflow boundaries are not affected. This approach reduces the entropy variations in one-dimensional flows by an order of magnitude, and in two dimensions the flow remains stable (nearly spherical) for many sound crossing times (defined here as the time for a sound wave to cross the radius of the accretion shock). In the two-dimensional case, the smooth flow is eventually disturbed by noise generated at the corners of the numerical grid, where the polar axes meet the inner radial boundary. This numerical noise ultimately seeds the feedback mechanism discussed in Section 5, and the flow becomes unstable.

## 4 STABLE FLOW IN ONE DIMENSION

Before jumping to the effects of aspherical accretion shocks, we first examine the stability of our idealized model in one dimension. Our goals here are two-fold: (1) to test the ability of our numerical model (in particular, our inner boundary conditions) to simulate a standing accretion shock and (2) to determine the stability properties of our solutions under radial perturbations.

If the initial conditions are set to match the inner boundary conditions (i.e., the initialized flow velocity at the inner boundary equals the imposed value enforced by the boundary conditions), our simulations are relatively uneventful. After initial transients die away, the solutions settle into a steady-state that closely matches the analytic solution.

To test the stability of these models, we tried various ways of perturbing the initial solution. In one set of experiments we created a mismatch between the initial conditions and the imposed flow velocity at the inner boundary. If the fixed boundary velocity is set smaller than the flow velocity in the initial analytic solution, the accretion flow is choked off at the inner boundary, and the buildup of pressure initially drives the accretion shock to larger radii. In another set of experiments, we added a large density perturbation in the supersonic free-fall portion of the initial conditions. When this over-dense shell reaches the accretion shock, the increased ram pressure initially drives the shock down to smaller radii.

The results of one of these attempts to perturb the SAS are shown in Figure 4. In all cases the shock front oscillated a few times, but these oscillations were always strongly damped. For the $`\gamma =4/3`$ solution shown here, the time for material to drift from the accretion shock to the inner boundary is $`\tau _d=17`$, compared to a sound crossing time of $`\tau _s=1.6`$ and a freefall time of $`\tau _f=0.65`$. The shock front rebounds from a perturbation on the relatively short sound-crossing timescale, $`\tau _s`$. For large perturbations this rebound typically overshoots the equilibrium shock radius and subsequent oscillations decay away after several $`\tau _s`$. The accretion shock radius always settled back to within $`0.5\%`$ of the original steady-state radius after several $`\tau _s`$. The mass accretion rate takes longer to settle down to the steady-state value, as the perturbed postshock flow gradually advects down through the leaky boundary. It is only after a few $`\tau _d`$ that the full effects of the perturbation have advected out of the region and the mass accretion rate at the inner boundary has relaxed to the steady-state value.

These simulations demonstrate that SAS models are stable to radial perturbations. As long as the flow remains spherically symmetric, the accretion shock remains at a fixed radius.

## 5 INSTABILITY IN TWO DIMENSIONS

To study the stability of standing accretion shocks under the assumption of axisymmetry we ran additional simulations in two dimensions with a variety of initial perturbations, including small density perturbations in the infalling material, weak pressure waves emanating from near the inner boundary, and small random velocity fluctuations in the postshock flow. To quantify the behavior of these accretion simulations our numerical code kept track of several global metrics, including the total kinetic energy (and the component in the angular direction), the thermal and gravitational potential energy of the shocked gas, the total vorticity, the average radius of the accretion shock, and the angular variance in the accretion shock radius.

All of our two-dimensional SAS simulations proved to be unstable. To illustrate the typical behavior of an axisymmetric accretion shock we will focus on a model with $`\gamma =4/3`$ and initialized with two rings (a density enhancement of 20%) placed asymmetrically in the preshock flow (see Figure 6). For sufficiently small perturbations, the early evolution is characterized by sound waves bouncing throughout the interior of the accretion shock. These waves are nearly spherically symmetric, but gradually take on an asymmetry characterized by a spherical harmonic with $`l=1`$. The strength of these asymmetric waves grows with time, and eventually they begin to significantly impact the global properties of the accretion shock. This evolution is shown in Figure 5, where we plot the total energy in angular motion (which we assume to be proportional to the local turbulent energy) and the angle-averaged shock radius. We also show the results for an unperturbed simulation to show that our numerical techniques can hold a SAS stable for at least several dynamical times. From this, we conclude that the introduction of asymmetric perturbations leads to a growing postshock turbulence in standing accretion shocks.

Also shown in Figure 5 is a similar simulation but with different inner boundary conditions: a reflecting wall plus a thin cooling layer. The initial radial solution is nearly identical to the adiabatic solution except close to the inner boundary, where the cooling becomes strong and the gas quickly settles into a thin layer on top of the hard surface. As shown in Figure 5, the evolution is qualitatively the same as in the adiabatic model with a leaky boundary. This comparison verifies the independence of these results from the details of the boundary conditions.

To further illustrate this early evolution we show a time series of the two-dimensional flow in Figure 6. We use the gas entropy to visualize the flow because a steady, spherical shock would produce a postshock region of constant entropy. Thus, any deviation from a spherical shock shows up as a change in the postshock entropy. Aside from the overall trend of growing perturbations in the postshock entropy, this figure illustrates the dominance of low-order modes. This is a generic result of our two-dimensional simulations. While the wave modes in the early evolution depend on the initial perturbations, the $`l`$ = 1 mode always grows to dominate the dynamics.

To better understand the origin of the SAS instability we show a spacetime diagram of angle-averaged quantities in Figure 7. Here we see two dominant trends: (1) pressure waves rising up from small radii at the local speed of sound to perturb the shock front and (2) regions of high transverse velocity originating at the shock and advecting downwards at the local flow speed (slower than the sound speed). It is the feedback between these two waves (at least for low values of $`l`$) that leads to amplification of the initial perturbations. Aspherical pressure waves rising up to the accretion shock distort the shock front from its equilibrium spherical shape. The slight obliquity of the shock front relative to the radial pre-shock flow then creates non-radial post-shock flow that advects inward, feeding turbulent flow at small radii. This growing turbulence in turn drives stronger aspherical pressure waves that propagate upwards and further distort the accretion shock. This feedback was discovered independently in the context of accretion onto a black hole by Foglizzo (2002), who dubbed it a “vortical-acoustic” cycle.

To contrast this with planar flow, any vorticity generated by perturbations to the planar shock front are advected downstream, unable to impact the shock dynamics. In the case of an SAS, however, the vorticity is trapped in a finite volume. With nowhere to go, the vorticity grows with time and eventually becomes strong enough to affect the shock front.

Once the turbulent energy density in the postshock flow becomes significant, the volume of the shocked gas begins to grow with no apparent bound. The origin of this increase in shock radius can be traced to the distribution (kinetic versus thermal versus gravitational) of total energy below the shock. Note that this is strictly adiabatic flow—there is no source of external heating to assist in driving the shock outwards. The total energy on the computational grid remains relatively constant (the total energy of the gas flowing off the inner edge of the grid is small and varies slightly due to the growing turbulence, while the total energy of the gas falling onto the grid from the outer boundary is approximately zero), but the distribution of energy changes. As the nominally spherical accretion shock becomes more distorted, the postshock gas retains more kinetic energy at the expense of thermal energy. This gradual change can be seen in a plot of the various components of the total energy (the three terms in eqn. 9) integrated over the volume of shocked gas. This is shown in Figure 8. Moreover, this redistribution of energy among kinetic, thermal, and gravitational components results in a redistribution of the total energy per unit mass in the postshock gas. In Figure 9 we see that, as the instability described here develops, an increasing fraction of the postshock gas attains positive energies at the expense of gas that becomes more bound. There is no new source of energy in this system, only a redistribution of energy.

To follow the long-term evolution of an SAS we continued the simulation of our standard model, with $`\gamma =4/3`$, but allowed the numerical grid to expand radially to accommodate the increasing shock radius. We followed the expansion out to a radius of 15 in this fashion. The evolution can be broken into three distinct epochs. At early times ($`t<40`$) the shock remains roughly spherical, and the volume of shocked gas is approximately constant while the interior turbulent energy grows steadily. The intermediate epoch is characterized by large oscillations in the postshock flow that lead to an aspect ratio for the accretion shock that varies between one and two. The volume of shocked gas exhibits a growing oscillation, reminiscent of a breathing mode. This phase of the evolution is dominated by the $`l`$ = 1 mode; i.e., is distinctly oscillatory. Beyond a time of 170, the flow takes on a strikingly self-similar form dominated by the $`l`$ = 2 mode; i.e., is distinctly bipolar. The aspect ratio becomes remarkably constant at a value $`2`$ (2.3 to be precise), but the scale of this aspherical shock grows roughly linearly with time. From $`t=165`$ to $`t=225`$, the angle-averaged shock radius increases from 5 to 15. This corresponds to an expansion velocity of order 50% the escape velocity at the shock. This self-similar expansion is illustrated in Figure 10, where we plot images of the gas entropy at three different times chosen such that the shock has roughly doubled in size in each interval. When scaled to the same size, the images look qualitatively similar.

The dynamics of this self-similar state are illustrated in Figure 11. The global features include a dominant accretion flow in the equatorial region and mildly supersonic, quasi-periodic outflows along the two poles. These outflows are fed by accreting gas that “misses the target.” The oblique shocks near the equator lead to low-entropy gas advecting in (as material hits the shock at an increasingly oblique angle, it does not suffer as large a change in entropy when traversing the shock).

The only two parameters in this idealized SAS model are the adiabatic index, $`\gamma `$, and the radius of the inner boundary, $`R_i`$. We ran several simulations varying these parameters and found qualitatively similar results in virtually all the simulations, but with two clear trends, illustrated in Figure 12. Simulations with a larger value of $`R_i`$, and hence less volume for postshock turbulence to develop, exhibited a slower growth in the total turbulent energy. Herant, Benz, & Colgate (1992) pointed out that the scale of convection in a confined spherical shell scales with the ratio of the inner to the outer boundaries. By analogy we might then expect the dominant mode of instability in the accretion shock to scale with $`R_i`$: the larger the value of $`R_i`$, the larger the value of $`l`$ of the dominant mode. While the growth of the instability was slower, all of our models did eventually evolve into a strong bipolar mode that drove the shock out to larger radii. Even the simulation with $`R_i=0.4`$ rapidly grew to a turbulent energy of 0.6 by a time of 80. We also found slower growth of the SAS instability with smaller values of $`\gamma `$. Only in the case of a very soft equation of state ($`\gamma =1.25`$) did the shock remain marginally stable for several flow times. Eventually, this model became as unstable as the others, with a growing bipolar accretion shock. While all of the SAS models are unstable, a decrease in the effective adiabatic index and spatial extent of the postshock flow does appear to slow the growth of the instability. Note that both of these trends might work against the growth of the instability in the evolution of core-collapse supernovae, where the effective adiabatic index is lower and the postshock volume is affected by neutrino cooling.

## 6 CONCLUSION AND DISCUSSION

Our primary conclusion is that relatively small-amplitude perturbations, whether they originate from inhomogeneities in the infalling gas, aspherical pressure waves from the interior region, or perturbations in the postshock velocity field can excite perturbations in a standing accretion shock that lead to vigorous turbulence and large-amplitude variations in the shape and position of the shock front.

It is important to note that this behavior is not driven by convection: the initial conditions in our SAS model are marginally stable to convection. Furthermore, during the course of our simulations the angle-averaged entropy remains relatively constant with radius. In Figure 6 we can see the development of convection on small scales from local entropy gradients in the flow below the accretion shock. However, there is no global gradient extending over the entire postshock region, driving convection on that scale.

Our models clearly show an evolution to low-order modes, as first seen in the early two-dimensional supernova models of Herant, Benz, & Colgate (1992). Independent of the initial seed, our models were eventually dominated by either an $`l=1`$ or $`l=2`$ mode after several flow times; these were the only modes that ever became highly nonlinear. Moreover, once these low-order modes began to dominate, the average shock radius began to grow, altering the global behavior from a standing accretion shock to an expanding aspherical blastwave.

Our work also suggests that large-scale ($`l`$ = 1 or 2) modes may dominate the late-time evolution of core-collapse supernovae, although our models cannot directly address this possibility because the assumption of supersonic freefall above the shock is no longer valid once the shock encounters the subsonic regions in the stellar core. Nonetheless, it is interesting to note that when our outflows become self similar, they achieve an aspect ratio $`2`$. This is within the range suggested by spectropolarimetry data (e.g., Höflich, Khokhlov, & Wang, 2001). Could it be that the polarization uniformly observed in core collapse supernovae is a signature of the unstable growth of $`l=1,2`$ modes in the initial explosion? Moreover, if the SAS instability is responsible for generating the large aspect ratios in observed explosions, it provides at least one example where a “conventional” neutrino-driven explosion model can be constructed that is not in conflict with these observations, a concern expressed by Wang et al. (2001).

While our results are relatively independent of the initial conditions, we note that the amplitude and scale of our seed perturbations are consistent with those expected to be produced by the late stages of shell burning in presupernova stars (Bazan & Arnett, 1998) or in the core collapse process itself. Bazan and Arnett followed the evolution of shell oxygen burning in a presupernova star and found that unstable nuclear burning can lead to significant density perturbations. When this layer of the progenitor star falls through the accretion shock, it can excite the instabilities found here. Alternatively, entropy-driven postshock convection (Herant et al., 1994; Burrows, Hayes, & Fryxell, 1995; Janka & Mueller, 1996; Mezzacappa et al., 1998; Fryer & Heger, 2000) could provide the necessary seeds for this shock instability. Yet another possibility is the growth of lepton-driven convection (Keil, Janka, and Mueller, 1996) or neutron fingers (Wilson & Mayle, 1993) deep in the proto-neutron star. If these grow to sufficient amplitude, they would drive aspherical pressure waves out into the accreting postshock gas, exciting this shock instability.

This instability exists! We have used time-dependent hydrodynamic simulations to show that a simple adiabatic model of a standing accretion shock that is stable in one dimension is unstable in two dimensions. The physical mechanism responsible for this instability appears to be the feedback between aspherical-shock–induced vorticity and the pressure waves generated by this vorticity. This feedback was also discovered in the context of accreting black holes and dubbed a “vortical-acoustic” feedback by Foglizzo (2002). We note that Houck & Chevalier (1992) performed a linear stability analysis of steady, spherical accretion flows with cooling. We have used one of these models to test the effect of different inner radial boundary conditions on our simulation outcomes. For $`\gamma =4/3`$, our one-dimensional simulations are in agreement with their results. In particular, we find that the shock is stable for small values of the shock stand-off distance (when the shock is less than 20 times the radius of the inner boundary surface) and unstable for larger values. However, for $`\gamma =4/3`$ Houck & Chevalier (1992) did not extend their analysis to the nonradial case; therefore, we cannot directly compare their results with our findings in this case.

Many questions must be answered before we can determine whether this vortical-acoustic instability plays an important role in core-collapse supernovae. Will neutrino cooling near the proto-neutron star surface dampen the feedback mechanism by damping the pressure waves? Will neutrino cooling below the shock drive the flow away from conditions under which an instability can develop by reducing the postshock volume? In the models considered here, radiative losses were neglected, although mimicked in one case by using a soft ($`\gamma =1.25`$) equation of state. In this case the instability was weakened but still developed into an aspherical flow.

Neutrino cooling and heating will have several important consequences for the development of the instability described here: (1) As mentioned above, neutrino cooling may affect the generation of pressure waves at the base of the density cliff associated with the surface of the proto-neutron star. (2) Neutrino cooling and heating will, in part, determine the postshock volume between the density cliff and the supernova shock, which in turn will, in part, determine which $`l`$ modes may grow in the postshock flow. (3) Neutrino heating may selectively drive modes other than the fastest growing, $`l=1,2`$, modes. Moreover, the effects due to neutrino cooling and heating will be time dependent. Once an explosion is underway, the cooling region will be confined to the surface of the proto-neutron star, although pressure waves may still be damped, and the heating region will increase in volume, allowing lower-$`l`$ modes to develop—perhaps the $`l=1,2`$ modes.

One can already begin to assess the relevance of the SAS instability to the core-collapse supernova problem by comparing our models with previous multidimensional core-collapse supernova models. In Mezzacappa et al. (1998), once the shock front is distorted by the initial neutrino-driven convection, the SAS instability does appear to affect the postshock flow a few hundred ms after bounce. The shock front is significantly distorted from the initial spherical shape, and this oblique shock feeds strong interior turbulence. However, unlike in our SAS models, the $`l=1`$ mode in fact did not come to dominate the postshock flow. Rather, $`l`$few modes dominated at the shock’s peak radius and later were “squeezed out” as the shock receded. In this case, owing to neutrino cooling and the lack of sufficient neutrino heating, the separation between the shock front and the surface of the proto-neutron star was never sufficient for low-$`l`$ modes of the SAS instability to develop. To test this hypothesis we ran a two-dimensional SAS model with an inner boundary at 0.6, comparable to the effective surface of the proto-neutron star relative to the shock front at late times in the simulations of Mezzacappa et al. (1998). While low-$`l`$ modes did exist, they were unable to grow significantly in such a confined geometry. This dependence on the shock radius with respect to the surface of the proto-neutron star further emphasizes the need for realistic neutrino heating and cooling—i.e., neutrino transport—in multidimensional supernova models.

The simulations of Burrows, Hayes, & Fryxell (1995) clearly exhibit strong turbulence fed by an aspherical shock at late times. However, they included only one quadrant in their explosion simulation, and, therefore, neither an $`l`$ = 1 or 2 mode was allowed. In the simulations described in Janka & Mueller (1996), one can again identify a close link between the deformed accretion shock and strong large-scale vorticies in the postshock flow. Furthermore, by the end of their simulations, at $`700`$ ms, it appears that the $`l=1`$ mode is beginning to dominate. An interpretation of the SPH models of Herant et al. (1994) in our context is complicated by the fact that their simulation is also confined to a $`90^o`$ wedge, cutting off the lowest order modes. Furthermore, their explosion time scale is relatively short, providing insufficient time for the long wavelength SAS instability to develop. This short explosion timescale is most likely the result of their neutrino transport implementation.

It is interesting to note that the instability we identify in this paper may provide a somewhat revised paradigm for powering core collapse supernovae. In the models we presented here, explosions were obtained in the absence of neutrinos, powered by gravitational binding energy. The SAS instability serves to redistribute the energy of the flow into kinetic energy at the expense of internal and gravitational energy. Thus, both the instability and neutrinos may serve as vehicles to convert gravitational binding energy into the kinetic energy of explosion. If the SAS instability does in fact play a role in the explosion mechanism in more realistic models, it would be interesting to compute, if possible, the contributions to the kinetic energy in the explosion that result from these two vehicles. Moreover, it has been pointed out that the energies of explosions powered by neutrino heating may be limited to $`13\times 10^{51}`$ erg (Janka, 2001). If the SAS instability sets in and also serves to power the supernova, we may end up with a more energetic explosion, one only partially powered by neutrinos. On another front, it is interesting to note that the SAS instability may help alleviate the precarious balance that is struck during neutrino shock reheating—between accretion luminosity, mass accumulation in the gain region, and ram pressure. As pointed out by Janka (2001), it is necessary to maintain sufficient neutrino luminosity and mass in the gain region to power an explosion, and this is achieved by maintaining a sufficient level of mass accretion through the shock and onto the proto-neutron star. On the other hand, this accretion is the source of the ram pressure the shock combats in its attempt to propagate outward in the stellar core. If neutrino heating is able to sustain a sufficiently large shock radius for a period of time, the SAS instability may act to initiate the explosion, removing some of the burden from the neutrinos.

There are other factors that may influence the role of shock instability in core-collapse supernovae. How will rotation and/or magnetic fields couple to this mechanism? For example, will rotation provide sufficient distortion of the accretion shock to help drive low-$`l`$ modes? And will this lead to an explosion where otherwise (in the absence of rotation) one would not occur? It may be that a combination of near-explosive conditions and rotation will be required before the instability we describe can develop. If it does, it will certainly contribute to the explosion and should be considered an ingredient in the explosion mechanism.

Our results are certainly influenced by the imposed axisymmetry of our two-dimensional simulations. For example, our bipolar outflows may be an artifact of our boundary conditions, although the aspect ratio we observe when our flows become self similar is certainly an independent feature; we could have obtained bipolar outflows with very different aspect ratios. Generally speaking, will the SAS instability exist in three dimensions? Will the low-order modes still dominate the flow? Will rotation serve the role in three dimensions that the imposed axisymmetry serves in two? Three-dimensional simulations are underway, and we hope to report on them soon.

J.M.B. and A.M. are supported in part by a SciDAC grant from the U.S. DoE High Energy and Nuclear Physics Program. A.M. is supported at the Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the U.S. Department of Energy under contract DE-AC05-00OR22725. C.D. was supported by an NSF REU program in the Department of Physics at North Carolina State University. We thank the North Carolina Supercomputing Center and the Center for Computational Sciences at ORNL for their generous support of computing resources. We acknowledge many useful interactions with Raph Hix, Alan Calder, Jeff Knerr, and Dana Paquin.


format: text

text:
# Spectroscopic signatures of a bandwidth-controlled Mott transition at the surface of 1T-TaSe2

## Abstract

High-resolution angle-resolved photoemission (ARPES) data show that a metal-insulator Mott transition occurs at the surface of the quasi-two dimensional compound $`1T`$-TaSe2. The transition is driven by the narrowing of the Ta 5d band induced by a temperature-dependent modulation of the atomic positions. A dynamical mean-field theory calculation of the spectral function of the half-filled Hubbard model captures the main qualitative feature of the data, namely the rapid transfer of spectral weight from the observed quasiparticle peak at the Fermi surface to the Hubbard bands, as the correlation gap opens up.

Electronic correlations can modify the electronic structure of solids not only quantitatively, but also qualitatively, inducing new broken-symmetry phases which exhibit charge, spin or orbital-order, and more exotic states in low dimensions. One of the most notable consequences of electronic correlations is the much studied metal-insulator (M-I) Mott transition . Recently, new theoretical approaches have considerably extended our understanding of this fundamental problem .

Many physical properties indirectly reflect the dramatic rearrangement of the electronic structure at the transition. Photoelectron spectroscopy, which probes the single-particle spectral function, can provide a direct view of such changes . However, comparing samples with different compositions faces materials problems like stoichiometry, defects, and disorder. A quantitative analysis is further complicated by the known surface sensitivity of the technique . An ideal experiment would record the energy and momentum-dependent spectrum, while tuning the crucial ($`W/U`$) parameter ($`U`$ is the on-site Coulomb correlation energy; $`W`$ is the bandwidth) in the same single crystal sample. Remarkably, it is possible to approach this ideal situation exploiting the occurrence of modulated structures (charge-density-waves; CDWs) in appropriate low-dimensional systems. There, the lattice distortion modulates the transfer integrals and therefore modifies the bandwidth. In materials that are close enough to a Mott transition, the reduced bandwidth may lead to an instability. There are strong indications for this scenario in the layered chalcogenide $`1T`$-TaS2, which presents a sharp order-of-magnitude increase of the resistivity at T=180 K , with a strong rearrangement of the electronic states . However, the complex phase diagram of the CDW in $`1T`$-TaS2 affects the electronic transition, which cannot be considered as a typical Mott transition.

Isostructural and isoelectronic $`1T`$-TaSe2 exhibits a similar CDW, but only one phase below TC=475 K. Its electrical resistivity remains metallic - albeit rather large - to very low temperatures , suggesting that the Se compound lies further from the instability than the S analog. Nevertheless, a transition could still occur at the crystal surface, where the $`U/W`$ ratio is expected to be larger as a result of smaller screening and coordination. The surface sensitivity of ARPES is ideal to investigate such an inhomogeneous state. In this Letter we present evidence for a Mott transition at the surface of $`1T`$-TaSe2. High-resolution temperature-dependent ARPES shows for the first time the disappearance of the coherent quasiparticle signatures at the Fermi surface and the opening of a correlation gap. These results are qualitatively well described by a dynamical mean-field (DMFT) calculation, and provide new insight into the spectral properties of the MI transition.

$`1T`$-TaSe2 has a layered structure, with the $`d^1`$ Ta atoms in a distorted octahedral environment. Adjacent layers interact weakly through van der Waals gaps, and all physical properties exhibit a strong anisotropy. A threefold CDW develops below $`T_C=475K`$, with a commensurate $`\sqrt{1}3\times \sqrt{1}3`$ superstructure, analogous to the much studied low-temperature CDW phase of $`1T`$-TaS2. In real space the CDW corresponds to a modulation of the atomic positions, within a 13 Ta atoms unit, in the so-called ‘star-of-David’ configuration. Extended Huckel calculations  suggest that the CDW splits the Ta d conduction band into subbands which contain a total of 13 electrons per unit cell. Two subbands, carrying 6 electrons each, are filled and lie below the Fermi level. The Fermi surface is formed by a half-filled subband carrying the 13th electron. The opening of a correlation gap in this subband is responsible for the resistivity jump in (bulk) $`1T`$-TaS2.

We performed high-resolution ARPES measurements in Lausanne and at the PGM beamline of the SRC, University of Wisconsin. The energy and momentum resolution were $`\mathrm{\Delta }`$E=10 meV and $`\mathrm{\Delta }k=\pm 0.02`$Å-1, and the Fermi level location was determined with an accuracy of $`\pm `$ 1 meV by measuring the metallic edge of a polycrystalline gold reference. Single crystal samples grown by the usual iodine transport technique were characterized by Laue diffraction and resistivity measurements, which confirmed the assignment to the 1T polytype. They were mounted on the tip of a closed-cycle refrigerator and cleaved at a base pressure of $`1\times 10^{10}`$ torr. We did not observe any sign of surface degradation or contamination over a typical 8 hours run.

Figure 1 shows ARPES intensity maps (h$`\nu `$ = 21 eV) measured at 300 K and 70 K along the high-symmetry $`\mathrm{\Gamma }`$M direction of the hexagonal Brillouin zone (BZ). A higher photon energy (h$`\nu `$ = 50 eV) yields similar results. The maps correspond to the same CDW phase, as confirmed by low-energy electron diffraction (LEED), but exhibit remarkable differences. At 300 K the narrow topmost Ta d subband crosses the Fermi level at $`k_{F_{1,2}}`$ $`\pm `$ 1/4 $`\mathrm{\Gamma }M`$, in good agreement with band structure calculations . The filled CDW subbands are visible at $``$ 0.3 eV and, with lower intensity at larger binding energies ($``$ 0.8 eV) and momenta. The overall dispersion of the Ta d band is influenced by the CDW superlattice, as will be discussed elsewhere. The parabolic band with a maximum at $`\mathrm{\Gamma }`$ and $``$0.5 eV is a Se p band. At 70 K the Ta d spectral weight is narrower and clearly removed from $`E_F`$, and a gap has appeared.

Bulk sensitive properties and surface sensitive LEED data rule out a structural phase transition between 300 K and 70 K. The large spectral changes are therefore the consequence of an electronic surface transition. Signatures of a surface gap in $`1T`$-TaSe2 were previously observed at 70 K by scanning tunneling spectroscopy (STS) . We characterized this transition by temperature-dependent measurements of the shallow Ta $`4f`$ core levels, which exhibit a CDW-induced fine structure . We find (not shown) that between 300 K and 70 K the CDW-split components sharpen, and their energy separation increases by 40 meV, in agreement with lower-resolution data . We conclude that the CDW amplitude and the corresponding lattice distortion are larger at the lower temperature.

The increased distortion reduces the overlap between the ‘cluster orbitals’ which form the basis of the band structure in the CDW phase. Calculations which explicitely account for this effect are not available, but the ($`W/U`$) ratio is certainly reduced, possibly below the critical value for the M-I transition. This is confirmed by an inspection of the ARPES signal at $`k=k_{F_1}`$ (Fig. 2a), which reveals a sudden loss of intensity near $`E_F`$ below $``$260 K, and the appearence of a strong signal centered at 0.26 eV. All spectra were normalized to the same integrated area. The intensity at $`E_F`$ (Fig. 2b) exhibits a sharp step around 260 K, and a further linear decrease at lower temperature. The intensity and the spectral lineshape (not shown) are recovered upon heating, but only at higher temperatures. The large hysteresis ($`80K`$) suggests a first-order transition .

In Fig. 3, following common practice in ARPES work on the cuprates, we symmetrized the spectra of Fig. 2a around $`E_F`$. This procedure removes the perturbing effect of the Fermi distribution on the intrinsic temperature dependence of the spectral function $`A(k_F,\omega ,T)`$. The 300 K spectrum exhibits a broad ($``$1 eV) incoherent background and a weak quasiparticle (QP) feature at $`E_F`$. The QP signal disappears at lower temperature, and spectral weight is transferred to the lower and - as inferred by symmetry - upper sidebands, representing the lower (LHB) and upper (UHB) Hubbard subbands. The integrated intensity is conserved. We emphasize that the width ($``$100 meV) of the central peak is much larger than expected for a coherent quasi-particle at the Fermi surface of a ‘good’ metal. Clearly, in the 220-260 K temperature range, the corresponding excitations are heavily scattered, and their lifetime is short. This ‘bad metal’ character is consistent with the broad maximum and large value of the electrical resistivity of bulk $`1T`$-TaSe2 at 250-300 K .

The spectra of Fig. 3 are qualitatively consistent with the changes expected at the Mott transition. In order to substantiate this, we have calculated $`A(k_F,\omega ;T)`$ for a one-band Hubbard model at half-filling, within the dynamical mean-field theory (DMFT) framework . The “iterated perturbation theory” approximation  was used, and checks were made using Quantum Monte Carlo and the maximum entropy method. The Coulomb term $`U`$ was set at 0.52 eV, equal to the energy separation between the LHB and UHB features in Fig. 3. A semi-circular density of states was used, with a bandwith $`W`$ assumed to depend linearly on temperature between 300 K ($`W`$=0.50 eV) and 70 K ($`W`$=0.36 eV). These values are only indicative, and are not the result of a specific attempt to find an optimum fit to the data, but we note that the overall magnitude of $`W`$ is consistent with the dispersions observed in Fig.1.

The high-temperature (300 K, 260 K) calculated spectra (Fig. 4) correspond to a correlated metal in the incoherent regime, with a broad low-energy peak and two intense Hubbard sidebands. The central peak is strongly reduced at 220 K, and at 70 K it has disappeared completely, leaving two sharp features centered at $`\pm U`$/2 and separated by a real gap. The overall shape of the spectra is in good qualitative agreement with the data, as well as the dramatic transfer of spectral weight that takes place between the central peak and the Hubbard bands as the temperature is lowered. We observe that in both theory and experiment the spectral weight accumulated in the insulator near the maximum of the HB sidebands comes mainly from the QP central peak, but that a small fraction also comes from energies above 0.5 eV. This results in two energies at which all spectra approximately cross. The data provide the first direct momentum-resolved observation of two of the key predictions of DMFT regarding the one-particle spectrum through the Mott transition, namely the three-peak structure in the metallic state  and the large transfers of spectral weight from the metal to the insulator

The comparison of theory and experiment also reveals some differences. Unlike the calculated spectra, the intensity at $`E_F`$ is never completely suppressed in the experimental spectra. This signal could originate from the underlying metallic bulk, or from surface inhomogeneities (metallic ‘patches’) observed by STS . It could also indicate that the tails of the LHB and UHB overlap slightly, a passibility which was considered in a different context . Another quantitative difference is the sharper separation between the QP peak and Hubbard bands in the calculation. This is a known feature of DMFT which is likely to be weakened as dimensionality is lowered. Recent work in particular  suggests that long-wavelength charge modes partially fill the preformed gap.

From a theoretical standpoint, the main discrepancy concerns the value of the critical temperature at which the first-order metal insulator transition is observed. Indeed, with the observed value of $`U`$, the simple one-band Hubbard model treated within DMFT would have a first-order transition at T$`U/8090`$ K, a factor of three smaller than observed experimentally. The trajectory in the (T,($`W/U`$)) space used in the theoretical calculation does not intersect the first-order transition line, so that theory would interpret the spectral changes as due to a rapid crossover between a bad metal and a Mott insulator. On the other hand the experimental observation of a hysteresis does suggest a real transition. In a purely electronic model, it is known that orbital degeneracy does lead to increased $`T_c`$ values . However, DMFT model calculations including the whole $`d`$ manifold, show a significant increase in $`T_c`$ only when the first subband is much closer to the Fermi level than in the experiment. Thus, it is unlikely that orbital degeneracy could explain the increased $`T_c`$. The inadequacy of a purely electronic model is also suggested by the large width of the low-energy QP. In a purely electronic model, this can only be accounted for if the temperature is significantly larger than the critical temperature for the metal-insulator transition.

The observation of a broad peak and a true transition strongly suggests that the coupling to lattice degrees of freedom plays an important role, as in fact expected in a CDW compound. Indeed, it has been shown in a simple model that the coupling to the lattice can lead to an increase in $`T_C`$ . Also, it is possible that the observed hysteresis results from differences in the pinning of the CDW upon heating and cooling. Obviously, further investigation - both theoretical and experimental - is required to clarify these issues.

In summary, we presented momentum-resolved high-resolution ARPES data which illustrate the spectral consequences of a bandwidth-controlled surface Mott transition in $`1T`$-TaSe2. The transition from a bad-metal, characterized by a largely incoherent spectrum, to a correlated insulator is qualitatively captured by a DMFT calculation for the half-filled Hubbard model. Quantitative differences between theory and experiment suggest that the model should be extended to include the coupling to lattice degrees of freedom, in order to provide a more accurate description of electronic transitions in such CDW materials.

We acknowledge R. Gaal for the resistivity measurements, correspondence with P. Aebi, E. Canadell, F. del Dongo, and the support of G. Margaritondo. This work has been supported by the Swiss NSF through the MaNEP NCCR. A.G is grateful to F.Mila for discussions and hospitality at UNIL-Lausanne. A.G, S.F and S.B are most grateful to KITP-UCSB for the warm hospitality during the final stage of this work, where it was supported in part by NSF under Grant PHY99-07949, and also acknowledge computing time at IDRIS-CNRS under project 011393. The Synchrotron Radiation Center, University of Wisconsin-Madison, is supported by NSF under Award No. DMR-0084402.


text:
# Mid-rapidity 𝜋^±, 𝐾^±, and 𝑝̄ spectra and particle ratios from STAR

## Abstract

Results are presented on $`\pi ^\pm `$, $`K^\pm `$, and $`\overline{p}`$ transverse mass spectra and particle multiplicity ratios at mid-rapidity in Au+Au collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$=130 and 200 GeV and in p+p collisions at $`\sqrt{s}`$=200 GeV. Comparisons are made to results from lower energies. The bulk properties of the collision inferred from these results are discussed.

New states in which the nucleon mass vanishes or quarks and gluons are deconfined have been predicted to exist at high densities over an extended volume . Such states may be created in relativistic heavy-ion collisions . Pions, kaons, and antiprotons are abundantly produced particles in relativistic heavy-ion collisions, reflecting the bulk properties of the collision: the yield of pions reflects on total entropy; kaons carry a significant fraction of the total strangeness produced; and the antiproton yield is a measure of baryon production. Measurements of these particles thus provide diagnostics to possible formation of the predicted new states .

We present mid-rapidity transverse mass ($`m_{}`$) spectra of charged pions ($`\pi ^\pm `$), charged kaons ($`K^\pm `$), and antiproton ($`\overline{p}`$) in Au+Au collisions at RHIC. We study spectral shapes in the picture of collective transverse radial flow and production of kaons and antiprotons relative to that of pions. We investigate the systematics of these results with respect to the collision centrality and in the context of results from lower energies to search for systematic changes which could result from possible changes in the collision dynamics.

The data presented here are from minimum bias Au+Au collisions at nucleon-nucleon center-of-mass energy of $`\sqrt{s_{_{\mathrm{NN}}}}`$=130 and 200 GeV and p+p collisions at $`\sqrt{s}`$=200 GeV by the STAR experiment . Tracks were reconstructed in the STAR Time Projection Chamber (TPC). The magnetic field was 0.25 and 0.5 Tesla for the 130 and 200 GeV data, respectively. The primary vertex of the interaction was found by fitting the tracks to a common point of origin. Tracks used in the analysis were required to come from within 3 cm of the primary vertex. Corrections were made for the energy loss of particles in the detector material. Particles were identified by measuring the specific energy loss ($`dE/dx`$) of charged particles in the TPC gas. The $`dE/dx`$ resolution was estimated to be 11% which imposed upper momentum limits on the identified particle spectra. Corrections were applied to account for losses due to limited acceptance, decays, tracking inefficiencies, and hadronic interactions. The centrality bins used in the analysis were obtained by subdividing the measured mid-rapidity charged particle multiplicity in the TPC .

Figure 1 shows the mid-rapidity $`\pi ^{}`$, $`K^{}`$, and $`\overline{p}`$ spectra for all centrality bins in 200 GeV Au+Au collisions. Systematic errors on the spectra are estimated to be 10%. The $`\pi ^+`$ and $`K^+`$ spectra are similar to the $`\pi ^{}`$ and $`K^{}`$ spectra, respectively. All the spectra are similar to the corresponding ones at 130 GeV . We fit the pion spectra with a Bose-Einstein function, the kaon spectra with a $`m_{}`$ exponential function, and the antiproton spectra with a $`p_{}`$ gaussian function, to extract the mean transverse momenta, $`p_{}`$, and the $`dN/d\mathrm{y}`$ yields. The fit results are superimposed in Fig. 1. The measured yields are about 65% of the extrapolated $`dN/d\mathrm{y}`$ for pions, 50-65% for kaons, and 50-75% for antiproton, respectively. Similar to , we use other functional forms to estimate the systematic errors on the extrapolated $`p_{}`$ and $`dN/d\mathrm{y}`$. The total systematic errors are 10% on $`p_{}`$, and 15% on pion and kaon $`dN/d\mathrm{y}`$ and 15-25% on antiproton $`dN/d\mathrm{y}`$.

Figure 1 also shows the mid-rapidity $`\pi ^{}`$, $`K^{}`$, and $`\overline{p}`$ spectra for minimum bias p+p collisions at 200 GeV. Primary vertex finding inefficiency, becoming more significant for lower multiplicity events, is not corrected for and contribute to the majority of the estimated systematic error of $`{}_{35}{}^{}{}_{}{}^{+10}`$% on the spectra. Same functional forms (as for the Au+Au data) are used to fit the p+p spectra to extract the $`p_{}`$ and $`dN/d\mathrm{y}`$. Possible systematic effects of the vertex finding inefficiency on $`p_{}`$ are investigated by studying multiplicity dependence of the spectra. The total systematic errors on $`p_{}`$ and $`dN/d\mathrm{y}`$ are estimated to be 10% and $`{}_{40}{}^{}{}_{}{}^{+15}`$%, respectively.

Figure 2 (left) shows the extracted $`p_{}`$ of $`\pi ^{}`$, $`K^{}`$, and $`\overline{p}`$ as a function of the mid-rapidity $`\pi ^{}`$ multiplicity density, $`dN_\pi ^{}/d\mathrm{y}`$, used as a measure of the collision centrality. A systematic increase with centrality is observed in the kaon and antiproton $`p_{}`$, consistent with collective transverse radial flow being built up in non-peripheral collisions. No significant change in $`p_{}`$ is observed from 130 to 200 GeV. Figure 2 (right) shows $`p_{}`$ versus particle mass for five selected systems. The p+p and most peripheral A+A results indicate no transverse radial flow; the increase in $`p_{}`$ with mass in these systems is due to a trivial mass effect. The central collision results deviate from p+p and are consistent with transverse radial flow, which appears to be stronger at RHIC than SPS. Transverse radial flow is a property inferred from particle spectral shapes at kinetic freeze-out. It remains a question at what stage of the collision this collective flow is built up.

Figure 3 (left) compares the mid-rapidity $`K^{}/\pi ^{}`$ ratio as a function of centrality at RHIC to those at the AGS  and SPS . The p+p results are also shown. At RHIC energies, little centrality dependence is found for $`K^{}/\pi ^{}`$ once a certain centrality is reached. The increase from p+p to A+A is achieved in peripheral collisions. This behavior was not seen in lower energy collisions. Both the AGS and SPS data dispatch a linear increase of $`K^{}/\pi ^{}`$ with $`dN_\pi ^{}/dy`$, suggesting that the physics might be different at RHIC than at lower energies. It has been suggested that the longitudinal thickness of the colliding nuclei may play a significant role at low energies and become less so at RHIC for medium-central to central collisions . Low energy RHIC data, expected in the coming years, should shed more light onto the physics.

Figure 3 (right) shows the $`\overline{p}/\pi ^{}`$ ratio as a function of centrality at AGS , SPS , and RHIC. Similar to $`K^{}/\pi ^{}`$, the $`\overline{p}/\pi ^{}`$ ratio at RHIC is constant over centrality from medium-central to central collisions. Unlike $`K^{}/\pi ^{}`$, however, the $`\overline{p}/\pi ^{}`$ ratios at the AGS and SPS are also flat, or perhaps decrease, with centrality . This might be related to antiproton absorption at later stages of the collision at the low energies. On the other hand, due to the relatively large production threshold, there is a dramatic increase in antiproton production from AGS to SPS to RHIC.

In summary, we have presented results from STAR on mid-rapidity $`\pi ^\pm `$, $`K^\pm `$, and $`\overline{p}`$ transverse mass spectra and $`K^{}/\pi ^{}`$ and $`\overline{p}/\pi ^{}`$ multiplicity ratios in Au+Au collisions at 130 and 200 GeV and p+p collisions at 200 GeV. No significant changes were observed from 130 to 200 GeV, and there is a smooth evolution from p+p to A+A collisions. The following bulk properties of the collision may be inferred from these results: (1) At kinetic freeze-out, significant transverse radial flow has developed in central Au+Au collisions. The magnitude of radial flow increases with the collision centrality and is larger at RHIC than SPS. (2) At chemical freeze-out, the $`K^{}/\pi ^{}`$ and $`\overline{p}/\pi ^{}`$ multiplicity ratios are higher in A+A than in p+p, with the increase achieved quickly but smoothly from p+p to peripheral A+A collisions, and remains relatively constant in medium-central to central collisions. This behavior is distinctly different from those observed at the AGS and SPS, suggesting a possible change in the production mechanisms of kaons and antiprotons with respect to pions from low energies to central collisions at RHIC energies. The present data may put stringent constraints on and provide inputs to the modeling of relativistic heavy-ion collisions.


text:
# The Formation of Stellar Clusters in Turbulent Molecular Clouds: Effects of the Equation of State

## 1 INTRODUCTION

Fragmentation of molecular clouds is key to the formation of star clusters, yet the conditions for fragmentation are poorly understood. There have been a few analytical approaches (Jeans 1902; Penston 1966; Low & Lynden-Bell 1976), and numerical investigations of the effects of various physical processes on the collapse, such as geometry and rotation of the clouds (see Bonnell & Bastien 1993 for a review), and magnetic fields (Galli et al. 2001). Recently, the effects of turbulence have been studied extensively in a series of 3D simulations (Klessen, Burkert, & Bate 1998, Paper I; Klessen, Heitsch, & Mac Low 2000, Paper II; Heitsch, Mac Low, & Klessen 2001, Paper III; and Klessen 2001, Paper IV). However, these results are based on isothermal models with $`\gamma =1`$. The effect of the equation of state (EOS), which is essential in understanding the physical structure and stability of the turbulent clouds, remains uncertain.

The balance of heating and cooling in a molecular cloud can be approximately described by the polytropic EOS, $`P=K\rho ^\gamma `$, where $`K`$ is a constant, and $`P,\rho `$ and $`\gamma `$ are thermal pressure, gas density and polytropic index, respectively. The density structures generated by supersonic turbulence depends on the EOS (Scalo et al. 1998; Passot & Vazquez-Semadeni 1998; Passot & Vazquez-Semadeni 1999). A detailed analysis by Spaans & Silk (2000) suggests that $`0.2<\gamma <1.4`$ in the interstellar medium.

In this paper, we perform a 3D numerical survey over this range of $`\gamma `$, to investigate the effects of the EOS on the fragmentation of turbulent clouds and the subsequent formation of a protostellar cluster.

## 2 SIMULATIONS

We carry out smoothed particle hydrodynamics (SPH) simulations to determine the effects of EOS on fragmentation. The code is implemented with periodic boundary conditions (Klessen 1997) and “sink” particles (Bate, Bonnell, & Price 1995) to replace high-density cores. The driven turbulence is set up following the $`k=`$1–2 model in Mac Low (1999). For detailed descriptions of the code, see Paper II. The treatment of the polytropic cloud in the code is the same as in Bonnell & Bastien (1991). We vary $`\gamma `$ in steps of 0.1 in otherwise identical simulations. Self-gravity is included in the calculation at $`t=2.0`$ after turbulence is fully established (see Paper II). Our models use 200,000 particles. The 13 simulations were performed simultaneously using a serial version of the code on processors of the Parallel Computing Facility of the AMNH.

### 2.1 Fragmentation

Figure 1 shows that the fragmentation properties and the accretion of the sinks strongly depend on the value of $`\gamma `$. The lower the $`\gamma `$, the earlier fragmentation occurs, and the more fragments form. It takes longer for gas with high $`\gamma `$ to fragment. The slope of the accretion curves are roughly the same, but the number of sinks is very different for different $`\gamma `$ at the same time, which suggests that individual protostellar cores is different for different $`\gamma `$. In low-$`\gamma `$ gas, there are more cores, but the mass of each core is much smaller than in high-$`\gamma `$ gas.

The Jeans mass for a polytropic gas sphere $`M_J=(K\pi /G)^{\frac{3}{2}}\gamma ^{\frac{3}{2}}\rho ^{\frac{3}{2}(\gamma \frac{4}{3})}`$ so

$$\frac{M_J}{\rho }\frac{3\gamma 4}{2}\rho ^{(3\gamma 6)/2},$$
(1)
which implies that $`\gamma =4/3`$ is a critical value during collapse. When $`\gamma <4/3`$, the Jeans mass decreases as density increases, promoting fragmentation, but when $`\gamma >4/3`$, the Jeans mass increases with density, stabilizing collapsing gas. Indeed, our model with $`\gamma =1.4`$ shows no collapse.

### 2.2 The Mass Spectra

Observations show that the mass distribution of gas clumps in various clouds follows a simple power law, $`dN/dM=M^\nu `$, with $`\nu 1.5`$ (Blitz 1993), while the stellar IMF follows $`\nu =2.33`$ at the high-mass end (Salpeter 1955).

Figure 2 shows that the spectra of both the clumps and the cores change with $`\gamma `$, but that the ones of cores change more dramatically. In low-$`\gamma `$ models, the mass distribution of the cores at the high-mass end is roughly log-normal. As $`\gamma `$ increases, fewer but more massive cores emerge. When $`\gamma >1.0`$, the distribution is dorminated by high mass cores only, and the spectrum tends to flatten out. It is no longer fitted by either a log-normal or a power-law. The clump mass spectra, on the other hand, do show power-law behavior at the high mass side, even for $`\gamma >1.0`$.

Our results suggest that stars tend to form in clusters in a low-$`\gamma `$ environment. Protostellar cores are low mass. The apparent lack of power-law behavior for the cores in the protostellar cluster might imply that simple accretion is unable to generate as many high-mass stars as predicted by the observations, which further suggests that other mechanisms such as collisions (Bonnell, Bate & Zinnecker 1998) may be at work to produce the massive stars in a cluster. Higher resolution models will be necessary to confirm this, however.

On the other hand, our results also imply that massive stars can form in small groups or alone in gas with $`\gamma >1.0`$. Spaans & Silk (2000) suggest that a stiffer EOS ($`\gamma >1.0`$) leads to a peaked IMF, biased toward massive stars, while an EOS with $`\gamma <1.0`$ results in a power-law IMF, in general agreement with our simulations.

## 3 ISOLATED STARS

The formation of isolated massive stars is of great interest, since in most cases, massive stars are found in clusters. But recently, Lamers et al. (2002) reported observations of isolated massive stars or very small groups of a few isolated massive stars in the bulge of M51. Massey (2002) also reported many field massive stars in both the Large and Small Magellanic Clouds. From our simulations, we see that when $`\gamma >1`$, only very few or possibly only one fragment occurs. This cloud core is massive, and would result in the formation of a high-mass star.

High resolution simulations by Abel, Bryan & Norman (2002) of the formation of the first star suggest that initially only one massive metal-free star forms per pregalactic halo. In the early Universe, inefficient cooling due to the lack of metals may result in high $`\gamma `$. Our models would then suggest weak fragmentation, resulting in the formation of only one massive star per cloud.

We thank M. Fall, H. Lamers, and H. Zinnecker for valuable discussions, and D. Janies and W. Wheeler for their work on the Parallel Computing Facility. YL thanks the AIP for its warm hospitality, and the Kade Foundation for support of her visit there. This research was partially funded by NASA grant NAG5-10103, and by NSF CAREER grant AST99-85392.


text:
# Computer simulations of domain growth in off-critical quenches of two-dimensional binary mixtures

## I Introduction

Fluid systems can be described at different scales. Molecular dynamics simulations have the advantage of taking into account elementary forces and microscopic details providing a realistic description of the fluid. On the other hand, coarse grained models, where for example a coarse-grained variable describes the concentration of the fluid over a mesoscopic region, are more convenient from a computational point of view, allowing simulations on time scales sometimes unreachable with molecular dynamics.

In phase separation of binary mixtures, for example, after the system is quenched by an initially disordered state into the coexistence region, the typical size of the domains of the different fluid components grows with power laws which are generally better measured in continuum coarse grained models. Due to this reason, in this paper we will study the phase separation of binary mixtures subjected to a shear flow describing the system by a continuous field that represents the difference in concentration between the two phases. Assuming local equilibrium, a chemical potential derived by a free-energy functional of the concentration field will be the source of the dynamical evolution. (This model is sometimes referred in literature as model B). A convection term due the presence of the shear flow will complete the description of the system.

The problem of phase separation of binary mixtures in a shear flow is relevant for many practical reasons especially concerning polymer mixtures and dilute polymer solutions. New recent results have shown that the shear flow greatly affects the morphology of the growing domains. Not only the growth is anisotropic with interfaces aligning with the direction of the flow, but for each direction two typical lenghtscales have been measured. This last feature makes a deep difference with what occurs in phase separation without imposed flows where the whole domain evolution is characterized by a single lengthscale $`R(t)`$.

In this paper we consider the case not studied before of mixtures with off-symmetric composition. Our model will be properly defined in Section II. The results will be described in Section III. Domains with two lenghtscales for each direction will be found also for the off-symmetrical mixtures. The role of the flow in the first stage of the phase separation will be examined together with the behaviour of the excess viscosity. Some conclusions will end the paper.

## II The model

In the following we model a binary mixture by coupling a diffusive field $`\phi `$, representing the normalized concentration difference between the two components of the mixture, to an applied velocity field. This approach neglects hydrodynamic effects. For weakly sheared polymer blends with large polymerization index and similar mechanical properties of the two species the present model is expected to be satisfactory in a preasymptotic time domain when velocity fluctuations are small. When hydrodynamic effects become important, instead, the full dynamical model, where the time evolution equation for the field $`\phi `$ is coupled to the Navier-Stokes equation, must be considered.

The mixture is described by the Langevin equation

$$\frac{\phi }{t}+\stackrel{}{}(\phi \stackrel{}{v})=\mathrm{\Gamma }^2\frac{\delta }{\delta \phi }+\eta $$
(1)
where $`\mathrm{\Gamma }`$ is a mobility coefficient, $`\stackrel{}{v}`$ is the external velocity field describing plane shear flow with profile $`\stackrel{}{v}=\gamma y\stackrel{}{e}_x`$, $`\gamma `$ and $`\stackrel{}{e}_x`$ being, respectively, the shear rate and the unit vector in the $`x`$ direction (flow direction). $`\eta `$ is a gaussian white noise, representing the effects of thermal fluctuations with zero mean and correlations that, according to the fluctuation-dissipation theorem, are given by

$$\eta (\stackrel{}{r},t)\eta (\stackrel{}{r}^{},t^{})=2T\mathrm{\Gamma }^2\delta (\stackrel{}{r}\stackrel{}{r}^{})\delta (tt^{})$$
(2)
where $`T`$ is the temperature of the mixture, and $`\mathrm{}`$ denotes the ensemble average. The term $`\stackrel{}{}(\phi \stackrel{}{v})`$ describes the advection of the field $`\phi `$ by the velocity and $`\mathrm{\Gamma }^2(\delta /\delta \phi )`$ takes into account the diffusive transport of $`\phi `$. The free-energy functional generally used in phase separation studies of binary mixtures is

$$=𝑑\stackrel{}{r}\left[\frac{a}{2}\phi ^2+\frac{b}{4}\phi ^4+\frac{\kappa }{2}(\phi )^2\right]$$
(3)
The first two terms are related to the bulk properties of the fluid. While the parameter $`b`$ is always positive, the sign of $`a`$ distinguishes between a disordered and a segregated mixture with phase-separated components. The case with $`a>0`$ gives a polynomial with one minimum in the origin corresponding to the a disordered state with $`\phi =0`$ everywhere while two symmetric minima are present when $`a<0`$ corresponding to the two pure phases with $`\phi =\pm \sqrt{a/b}`$. These are the equilibrium values of the order parameter $`\phi `$ in the small temperature limit. In the following we will consider a deep quench at $`T=0`$ into the two-phase coexistence region. The qualitative behaviour in the non-zero temperature case is similar: Thermal fluctuations are responsible for the increased roughness of the interfaces with respect to $`T=0`$ and for the thermal excitations inside the ordered domains . The gradient term is related to the interfacial properties and takes into account the energy cost for interface formation. The equilibrium profile between the two coexisting bulk phases is (in a one-dimensional system) $`\phi (x)=\mathrm{tanh}\sqrt{{\displaystyle \frac{1}{2\kappa }}}x`$, giving a surface tension equal to $`{\displaystyle \frac{2}{3}}\sqrt{2\kappa }`$ and an interfacial width proportional to $`\sqrt{2\kappa }`$.

The thermodynamic properties of the fluid follow from the free energy (3). The chemical potential difference between the two fluids is given by

$$\frac{\delta }{\delta \phi }=a\phi +b\phi ^3\kappa ^2\phi .$$
(4)
Equation (1) can be cast in a dimensionless form after a redefinition of time, space and field scales. Then we have chosen $`\mathrm{\Gamma }=|a|=b=\kappa =1`$.

The main observable for the study of the growth kinetics is the structure factor $`C(\stackrel{}{k},t)=\psi (\stackrel{}{k},t)\psi (\stackrel{}{k},t)`$, namely the Fourier transform of the real-space equal time correlation function. $`\psi (\stackrel{}{k},t)`$ is the Fourier transform of the concentration fluctuations with respect to its average value: $`\phi \phi `$. From the knowledge of the structure factor one computes the average size of domains. In the unsheared case ($`\gamma =0`$ in Eq. (1)) this quantity can be defined as

$$R(t)=\pi \frac{𝑑kC(k,t)}{𝑑kkC(k,t)}$$
(5)
where

$$C(k,t)=\frac{1}{N_k}\underset{|\stackrel{}{k}|=k}{}C(\stackrel{}{k},t)$$
(6)
is the spherical average of the structure factor and $`N_k`$ is the number of lattice vectors of length $`k`$ in the reciprocal space of momenta. In the case with shear, due to anisotropy, we measure

$$R_x(t)=\pi \frac{𝑑\stackrel{}{k}C(\stackrel{}{k},t)}{𝑑\stackrel{}{k}|k_x|C(\stackrel{}{k},t)}$$
(7)
and analogously for the other directions. An excess of viscosity $`\mathrm{\Delta }\eta `$ can be measured in sheared mixtures due to the work spent in stretching the domains. $`\mathrm{\Delta }\eta `$, which is of experimental interest, is defined as

$$\mathrm{\Delta }\eta =\frac{1}{\gamma }\frac{d\stackrel{}{k}}{(2\pi )^d}k_xk_yC(\stackrel{}{k},t).$$
(8)
where $`d`$ denotes the number of spatial dimensions.

We have simulated Eq. (1) in $`d=2`$ by a first-order Euler discretization scheme. Periodic boundary conditions have been implemented in the $`x`$ direction; Lees-Edwards boundary conditions were used in the $`y`$ direction (shear direction). These boundary conditions, originally developed for molecular dynamic simulation of fluids in shear, consist of moving the top and bottom periodic images of the lattice with respect to it so that they require the identification of the point at $`(x,0)`$ with the one located at $`(x+\gamma L\mathrm{\Delta }t,L)`$, where $`L`$ is the size of the lattice (the same in all the directions) and $`\mathrm{\Delta }t`$ is the time discretization interval.

## III Results

Simulations were run using lattices of size $`L=1024`$ with grid size $`\mathrm{\Delta }x=0.5,1`$. We do not observe significative differences between these two choices of $`\mathrm{\Delta }x`$. The results here shown were obtained with $`\mathrm{\Delta }x=1`$, $`\mathrm{\Delta }t=0.01`$, $`T=0`$ and three different values of $`\gamma =0,0.00488,0.0488`$. The system was initialized in a high temperature disordered state with $`\phi (\stackrel{}{r},0)=\phi _0+f(\stackrel{}{r})`$, where $`f(\stackrel{}{r})`$ is a random number in the range $`(0.01,0.01)`$ such that $`\phi (\stackrel{}{r},0)=\phi _0`$. The evolution was studied with $`\phi _0=0.4`$, corresponding to a composition $`70:30`$ of the mixture.

### A Unsheared quench

In this section we report the results for the unsheared quench. In Fig. 1 we show a sequence of configurations at consecutive times. The minority phase is made of isolated droplets whose number decreases, while their average radius increases in time. The main physical mechanism for coarsening in the present case, when hydrodynamics is neglected and thermal fluctuations are absent, is the evaporation-condensation mechanism proposed by Lifshitz and Slyozov. In this case larger domains grow at the expenses of smaller domains due to diffusion of atoms. This gives a growth exponent $`\alpha =1/3`$. To check this prediction we computed the average size of domains $`R(t)`$ using Eq. (5). The plot showing the time dependence of $`R(t)`$ is presented in Fig. 2. After an initial regime when well defined domains are forming, corresponding to the fact that $`\phi `$ is approaching its equilibrium values either $`1`$ or $`1`$, $`R(t)`$ shows from $`t500`$ onwards an asymptotic growth consistent with an exponent $`\alpha =1/3`$. Due to the large system used we could eliminate finite-size effects and observe the $`1/3`$ regime for two decades in time having a convincing evidence of the mentioned growth mechanism.

### B Sheared quench

The effect of the advective term, due to shear, on the segregation process is investigated in this section. We used two values of the shear rate $`\gamma =0.00488,0.0488`$. In Figs. 3-4 a time sequence of configurations for each value of $`\gamma `$ is shown. Snapshots are at the same values of the shear strain $`\gamma t`$ in order to have an immediate comparison between the two runs. For strains less than or of the order of one the effects of shear on the shape of bubbles are negligible: the mean sizes of the domains for a given shear rate remain directionally independent. However for the smaller shear rate the system generates bigger droplets because there is more time, keeping the strain fixed, for growth before the shear can influence the morphology.

This behaviour can be observed quantitatively computing the typical sizes of domains according to Eq. (7). The evolution of $`R_x,R_y`$ vs. $`\gamma t`$ is shown in Figs. 5-6. The ratio $`R_x/R_y`$ stays constant ($`1`$) until $`\gamma t1`$ for both the shear rates, then a crossover is observed. Moreover, we see that at $`\gamma t1`$, $`R_y(\gamma =0.00488)3.7`$ and $`R_y(\gamma =0.0488)2.6`$.

When the strain is larger than one, there is a competition between coarsening and shear-induced deformations of domains. For the smaller shear rate case, the diffusion-evaporation mechanism is still important and only domains which could grow in size are affected by the flow, while smaller ones evaporate. When the shear rate is larger, on the other hand, patterns grow smaller compared to those in the smaller shear rate case and diffusion from small to large droplets is inhibited as it can be seen in Fig. 4 at $`\gamma t=3`$, where the minority phase is more dispersed, consisting of a large number of domains with smaller surfaces. In both the cases the growth is faster in the flow direction and domains assume a striplike shape aligned with the flow. As the elongation increases, non-uniform patterns appear. This phenomenon was observed also in the case of sheared critical quenches and was interpreted as due to the presence in the system of two typical length scales in each direction, which compete in time. At $`\gamma t=8`$ there is an abundance of stretched domains which is also reflected in the local minima of $`R_y`$ (see Figs. 5-6). A convincing evidence of such a behaviour is given by the structure factor $`C(\stackrel{}{k},t)`$ which is contour-plotted at $`\gamma t=8`$ in Fig. 7 for the case $`\gamma =0.00488`$ (a similar behaviour is observed also for the other value of $`\gamma `$). It is separated in two distinct foils, each of them having two peaks. The presence of a four-peaked structure factor in sheared mixtures is a quite general feature which can be also shown mathematically . Due to the symmetry $`C(\stackrel{}{k},t)=C(\stackrel{}{k},t)`$, the foils are symmetric and one can consider only the two peaks of one foil. Let us consider the peaks which are in the half-plane $`k_x<0`$. The peaks are located at $`(k_{x_1},k_{y_1})(0.035,0.146)`$ and $`(k_{x_2},k_{y_2})(0.080,0.310)`$. The interpretation of a peak in the structure factor as the signature of the existence in the system of a characteristic length scale, indicates that there are two relevant length scales for each direction in the present case proportional to $`\pi /|k_{x_{1,2}}|`$ and $`\pi /|k_{y_{1,2}}|`$ in the flow and in the shear direction, respectively. This gives a value of around 2.3 and 2.1 for the ratio between the characteristic sizes of domains in the flow and in the shear direction, respectively. In the case of sheared quench of a symmetric mixture these ratios were 2.5 and 2.1 . Considering the behaviour in the shear direction, we see that there are domains with two different thicknesses in the system (the behaviour in the flow direction is analogous). At $`\gamma t=8`$ the peak which prevails is the one located at $`(k_{x_2},k_{y_2})`$, clearly indicating a larger abundance of thin domains, as already observed. Then, when the strain increases, domains are stretched and a cascade of ruptures occurs in those regions where the stress is higher and thicker domains, not yet broken, prevail as at $`\gamma t=15`$. This mechanism produces an oscillation in time in the behaviour of the peaks of the structure factor and now it is the other peak, located at $`(k_{x_1},k_{y_1})`$, to prevail. This is shown in Fig. 8 where a three-dimensional plot of one foil of the structure factor is pictured at two consecutive times. This mechanism is similar to the one observed in the case of symmetric mixtures under shear flow . Such an oscillation affects also $`R_y`$ whose amplitude is reduced at higher shear rates. We investigated also the existence of a scaling regime, analogously to the one observed in the unsheared case. Due to the inevitable finite-size effects, simulations cannot give a convincing evidence of such a scaling regime. We showed in a previous paper, using a renormalization group scheme, that if there is scaling, then $`R_x\gamma t^{4/3}`$ and $`R_yt^{1/3}`$. This means that the growth in the shear direction is not affected by the external velocity while in the flow direction there is an extra contribution of $`1`$ to the growth exponent coming from the advective term in Eq. (1). Our runs do not allow to check this prediction.

We also considered the rheological properties of the mixture. The excess viscosity was measured according to Eq. (8). The plots of $`\mathrm{\Delta }\eta `$ as a function of $`\gamma t`$ are shown in Figs. 9-10. Starting from zero, $`\mathrm{\Delta }\eta `$ shows a net increase up to a global maximum at $`\gamma t5`$, independent of the shear rate. This results from the stretching of domains caused by shear, which requires work against the surface tension. In the case with $`\gamma =0.00488`$ (Fig. 9) this increase is not monotonic and $`\mathrm{\Delta }\eta `$ shows a first peak at $`\gamma t<1`$. This first peak disappears by increasing $`\gamma `$. The height of the main peak depends on the shear rate being higher at smaller shear rates. This is due to the morphology of domains which have more interfaces for smaller $`\gamma `$ (cfr. Figs. 3-4 at $`\gamma t=5`$), giving a larger contribution to $`\mathrm{\Delta }\eta `$. For larger times $`\mathrm{\Delta }\eta `$ decreases due to the breakup of domains which dissipates energy previously stored. An oscillation is completed at $`\gamma t15`$, which is more evident for $`\gamma =0.00488`$, and it is related to the corresponding oscillation in $`R_y`$. If scaling were obeyed, one would expect $`\mathrm{\Delta }\eta `$ to scale as $`\mathrm{\Delta }\eta R_x^1R_y^1`$. Given the above discussed behaviour of $`R_x`$ and $`R_y`$, it should be $`\mathrm{\Delta }\eta t^{5/3}`$. This power law is roughly consistent with our data, especially for $`\gamma =0.00488`$.

## IV Conclusions

We studied the effects of a shear flow on the phase separation of binary mixtures with asymmetric composition. The binary mixture is described by a continuum free-energy functional and we have solved numerically the related model B. Simulations were carried out on two-dimensional systems; lattices are large enough to minimize finite-size effects.

Without shear we observe an asymptotic scaling regime with a growth exponent $`\alpha =1/3`$ consistent with previous theoretical and numerical studies. Thanks to the size of the lattice, we could observe such an exponent over two time decades.

More interesting and rich is the phenomenology when an external flow is applied during the quench. The patterns are stretched and aligned along the flow direction. Besides this behaviour, domain growth competes with shear-induced deformations and at high shear rates the minority phase appears more dispersed and made of smaller domains. In the limit of very high shear rates one could observe a complete demixing. Moreover, domains with two different thickness are clearly visible in the simulations: we observe a larger abundance first of thin and, then, of thicker domains, analogously to what observed in sheared critical quenches. The existence of two typical length scales in each spatial direction is related to a four-peaked structure factor. The alternate predominance of these length scales produces an oscillation between the peaks of the structure factor and in the physical observables of the system.

The rheological properties were investigated through the excess viscosity. This quantity increases up to a maximum and then decreases. This behaviour is associated to the stretching of domains, which requires work against the surface tension, and to the following ruptures of patterns which dissipates the stored energy. The position of the main peak as a function of the strain does not depend on the shear rate but its height decreases increasing $`\gamma `$ due to the more dispersed morphology of the minority phase.

Finally, we want to mention that we could not give a definitive answer to the question, which is still matter of an open debate, whether growth continues indefinitely with shear or a steady state is reached. Related to this question there is the problem of a reliable determination of the growth exponents.

Acknowledgments

G.G. acknowledges partial support by INFM PRA-HOP 1999. A.L. acknowledges partial support by INFM PA-G01-1.