format: text

text:
# The structure of radiative shock waves

## 1 Introduction

Main indicators of shock waves in atmospheres of radially pulsating stars are enhanced hydrogen emission, double profiles of absorption lines and discontinuities in radial velocity curves. Most prominently these features are observed in W Vir (Lèbre & Gillet 1992), RV Tau (Gillet et al. 1989) and Mira type (Alvarez et al. 2000) pulsating variables. The ratio of the atmosphere scale height to the stellar radius in these stars is $`H/R10^2`$ and the linear theory of stellar pulsation admits the existence of progressive waves propagating through the stellar atmosphere (Unno 1965; Pijpers 1993). The running waves must inevitably transform into shock waves due to nonlinear effects and this conclusion is corroborated by hydrodynamic calculations of radial pulsations of low–mass late–type giants (Wood 1974; Fadeyev & Tutukov 1981; Bowen 1988; Fadeyev & Muthsam 1990).

Though the role of shock waves as the principal mechanism of mass loss from late–type stars is still disputable (Hoefner & Dorfi 1997) it is nevertheless undoubted that shock–driven gas flows are responsible for formation of circumstellar shells. It must be noted that the postshock velocity near the photosphere should not necessarily be greater than the local escape velocity, because periodic passage of shock waves distends the stellar atmosphere and, under some conditions, may lead to the gradual approach of the gas flow to the escape velocity (Willson & Hill 1979). Furthermore, the shock amplitude tends to increase while the shock wave propagates in a medium with decaying gas density.

The most important problem concerning the shock–driven mass loss is the rate of the shock decay during its passage through the stellar atmosphere. Various simplifying assumptions lead to estimates of mass loss rates that differ from one another by many orders of magnitude (Wood 1979). Therefore, an analysis must be based on the self–consistent solution of the radiation transfer, fluid dynamics and rate equations. Unfortunately, the description of the shock wave passage in terms of the Lagrangean approach cannot be applied, since it does not allow the postshock relaxation zones to be explicitly resolved with respect to the spatial coordinate (Klein et al. 1976, 1978). Another approach is based on the assumption of steady–state gas flow. Applicability of this assumption for shock waves in stellar atmospheres is justified by the small width (including both the radiative precursor and the postshock relaxation zone) of the shock wave in comparison with the atmosphere scale height.

Another necessary condition for applicability of the steady–state assumption implies that the characteristic time scale of the shock wave decay is $`t_{\mathrm{dec}}d/U`$, where $`d`$ and $`U`$ are the width and velocity of the shock wave. In atmospheres of radially pulsating stars this condition is obviously fulfilled because emission lines indicating the shock wave are observed during the major part of the pulsation period (Abt 1954; Willson 1976) and in some stars this interval is as long as a half (Fox et al. 1984) or even 0.8 (Gillet et al. 1985) of the pulsation period.

Even in the framework of the steady–state assumption the problem of the shock wave structure in the partially ionized gas remains tremendously difficult, due to the strong coupling between the radiative precursor and the shock wake. Indeed, the occupation number densities of bound and free levels in both atoms and molecules strongly depend on the radiation field produced in the shock wake and at the same time must be determined as a solution of rate equations that are stiff on the relevant time scales. Fortunately, in stars with effective temperatures of $`T_{\mathrm{eff}}3000\mathrm{K}`$, the gas is mostly in the atomic state and the principal dissipative processes in the postshock region are excitation and ionization of atoms.

In Paper I (Fadeyev & Gillet 1998) we proposed the method of global iterations that allows us to obtain the stable self–consistent solution of the equations of fluid dynamics, radiation transfer and rate equations for the steady–state plane–parallel shock wave propagating through the atomic hydrogen gas. In Paper II (Fadeyev & Gillet 2000) this method was employed for computations of the structure of shocks with upstream velocities of $`15\mathrm{km}\mathrm{s}^1U_170\mathrm{km}\mathrm{s}^1`$ and for the temperature and the density of the unperturbed hydrogen gas of $`T_1=6000\mathrm{K}`$ and $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$, respectively. Results of these calculations are obviously insufficient for astrophysical applications since the temperature and the density in atmospheres of pulsating stars vary in wide ranges.

Below we discuss the structure of shock waves propagating through partially ionized hydrogen gas with temperature and density ranges within $`3000\mathrm{K}T_18000\mathrm{K}`$ and $`10^{12}\mathrm{gm}\mathrm{cm}^3\rho _110^9\mathrm{gm}\mathrm{cm}^3`$, respectively. Basic equations and the method of their solution are described in Paper II. In comparison with our previous work, we improved the treatment of radiation transfer and extended the total spectral range to $`13\mathrm{log}\nu 16`$. Furthermore, we employed the much more correct treatment of the spectral line radiation transfer which is of great importance for atomic level populations obtained as solution of rate equations. All calculations were done for the hydrogen atom with $`L=4`$ bound levels and a continuum.

## 2 General description of the shock wave structure

In the comoving fluid frame the radiation–modified Rankine–Hugoniot relations are written (Marshak 1958; Mihalas & Mihalas 1984) as

$$\rho U=𝒞_0\dot{m},$$
(1)
$$\dot{m}U+P_\mathrm{g}+P_\mathrm{R}=𝒞_1,$$
(2)
$$\frac{1}{2}\dot{m}U^2+\dot{m}h+F_\mathrm{R}+U\left(E_\mathrm{R}+P_\mathrm{R}\right)=𝒞_2,$$
(3)
where $`\rho `$ is the gas density, $`U`$ is the gas flow velocity, $`h`$ is the specific enthalpy, $`P_\mathrm{g}`$ is the gas pressure, $`E_\mathrm{R}`$, $`F_\mathrm{R}`$ and $`P_\mathrm{R}`$ are the radiation energy density, radiation flux and radiation pressure, respectively.

In planar geometry the total energy flux $`𝒞_2`$ given by relation (3) is constant along the spatial coordinate $`X`$, so that while the parcel of gas passes the shock wave the decrease of the flux of kinetic energy $`F_\mathrm{K}=\frac{1}{2}\dot{m}U^2`$ is balanced by changes in the enthalpy flux $`F_\mathrm{h}=\dot{m}h`$ and the radiation flux $`F_\mathrm{R}`$. Throughout the shock wave the flux $`U(E_\mathrm{R}+P_\mathrm{R})`$ is significantly smaller than other terms of the left–hand side of relation (3) and can be ignored.

As in Paper II the shock wave model is represented by the comoving plane–parallel finite slab. The space coordinate $`X=0`$ is set at the viscous adiabatic jump which is treated as an infinitesimally thin discontinuous jump where hydrodynamic variables undergo an abrupt change. It is assumed that the space coordinate $`X`$ of the gas element increases while the gas flows through the shock wave, so that in the preshock region $`X<0`$ and in the postshock region $`X>0`$. The coordinates of the preshock and postshock outer boundaries are denoted as $`X_1`$ and $`X_N`$, respectively. The upstream radiation flux emerging from the preshock outer boundary is negative, whereas the downstream radiation flux emerging from the postshock outer boundary is positive, that is, $`F_{\mathrm{R1}}<0`$ and $`F_{\mathrm{R}N}>0`$.

On the upper panel of Fig. 1, we show the fluxes $`F_\mathrm{K}`$, $`F_\mathrm{h}`$ and $`F_\mathrm{R}`$ as a function of space coordinate $`X`$ for the shock wave model with $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$, $`T_1=6000\mathrm{K}`$ and $`U_1=60\mathrm{km}\mathrm{s}^1`$. For the sake of convenience we use the logarithmic scale along the spatial coordinate $`X`$, the preshock and postshock regions being represented by the left plot and the right plot, respectively. Below the upper panel of Fig. 1, we show as a function of $`X`$ the fractional number density of hydrogen atoms in 2–nd state $`n_2/n_\mathrm{H}`$, the hydrogen ionization degree $`x_\mathrm{H}=n_\mathrm{e}/n_\mathrm{H}`$, the electron temperature $`T_\mathrm{e}`$ and the temperature of heavy particles (i.e. of neutral hydrogen atoms and hydrogen ions) $`T_a`$, and on the lowest panel we give the plot of the compression ratio $`\rho /\rho _1`$.

We assume that at the preshock outer boundary $`X_1`$, the gas temperature and the gas density are the same as those of the unpertubed medium that are denoted as $`T_1`$ and $`\rho _1`$, respectively. At the same time, the upstream radiation flux emerging from the shock wake affects both the bound–level number densities and the ionization degree of hydrogen atoms. Therefore, the number densities of hydrogen atoms in $`i`$–th state $`n_i`$ and the number density of free electrons $`n_\mathrm{e}`$ at the preshock outer boundary are determined from the solution of the equations of statistical equilibrium. In our calculations we tried to set the preshock outer boundary $`X_1`$ as far as possible from the discontinuous jump in order to diminish the perturbing influence of the shock wave radiation field on the preshock gas.

The most prominent changes that the gas flow undergoes in the preshock region are due to the absorption of the Lyman continuum radiation. In Fig. 1 the radiative precursor is revealed by the growth of the enthalpy flux $`F_\mathrm{h}`$ and due to the corresponding changes of the radiation flux $`F_\mathrm{R}`$. The length of the radiative precursor $`\delta X_\mathrm{p}`$ is proportional to the mean free path of Lyman continuum photons and increases with decreasing density from $`\delta X_\mathrm{p}10^3\mathrm{cm}`$ at $`\rho _1=10^9\mathrm{gm}\mathrm{cm}^3`$ to $`\delta X_\mathrm{p}10^7\mathrm{cm}`$ at $`\rho _1=10^{12}\mathrm{gm}\mathrm{cm}^3`$.

Within the radiative precursor the divergence of the Lyman continuum flux is negative, that is,

$$𝐅_{\mathrm{Lyc}}=4\pi \underset{\nu _1}{\overset{\mathrm{}}{}}\left(\eta _\nu \kappa _\nu J_\nu \right)𝑑\nu <0,$$
(4)
where $`\eta _\nu `$ and $`\kappa _\nu `$ are the monochromatic emission and absorption coefficients, $`J_\nu `$ is the mean intensity and $`\nu _1`$ is the threshold frequency for ionization from the ground state of the hydrogen atom. Thus, the gas temperature increases due to absorption of the Lyman continuum radiation while the gas parcel approaches the discontinuous jump. As a result, the gas density $`\rho `$ just ahead of the discontinuity exceeds its unperturbed value $`\rho _1`$, whereas the gas flow velocity $`U`$ becomes smaller than $`U_1`$. The change of the gas density $`\rho `$ and the gas flow velocity $`U`$ in the radiative precursor depends not only on the upstream velocity $`U_1`$ but also on the temperature and density of the unperturbed preshock gas. At $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$ and $`U_1=75\mathrm{km}\mathrm{s}^1`$, the relative growth of the density within the radiative precursor ranges from $`\mathrm{\Delta }\rho /\rho _10.03`$ at $`T_1=3000\mathrm{K}`$ to $`\mathrm{\Delta }\rho /\rho _10.07`$ at $`T_1=7000\mathrm{K}`$.

It should be noted that when the hydrogen gas is not completely ionized in the preshock region, the increase in the enthalpy flux is mostly due to ionization of hydrogen atoms. The rates of collisional transitions, both bound–bound and bound–free, are negligible in comparison with those of radiative transitions, so that the increase in $`F_\mathrm{h}`$ ahead of the discontinuous jump results mostly from the photoionization of hydrogen atoms.

Across the discontinuity, the flux of the kinetic energy $`F_\mathrm{K}`$ undergoes an abrupt decrease which is exactly balanced by the same abrupt increase of the enthalpy flux $`F_\mathrm{h}`$. Because population number densities per unit mass $`n_i/\rho `$ and $`n_\mathrm{e}/\rho `$ do not change across the discontinuity, the increase in the specific enthalpy is only due to the abrupt increase in the temperature of heavy particles $`T_a`$ and the electron temperature $`T_\mathrm{e}`$. Here we assume that the temperature of heavy particles just behind the discontinuous jump is given by the solution of the Rankine–Hugoniot relations (1) – (3):

$`T_a^+`$ $`=`$ $`T_a^{}x_\mathrm{H}^{}\left(T_\mathrm{e}^+T_\mathrm{e}^{}\right)+{\displaystyle \frac{1}{5}}{\displaystyle \frac{\dot{m}U^{}}{n_\mathrm{H}^{}k}}{\displaystyle \frac{\eta ^21}{\eta ^2}}`$ (5)
$``$ $`{\displaystyle \frac{2}{5}}{\displaystyle \frac{F_\mathrm{R}^+F_\mathrm{R}^{}}{n_\mathrm{H}^{}kU^{}}}`$
$``$ $`{\displaystyle \frac{2}{5}}{\displaystyle \frac{(E_\mathrm{R}^++P_\mathrm{R}^+)\eta ^1(E_\mathrm{R}^{}+P_\mathrm{R}^{})}{n_\mathrm{H}^{}k}},`$
whereas the electron temperature increases due to adiabatic compression:

$$T_\mathrm{e}^+=\eta ^{\gamma 1}T_\mathrm{e}^{},$$
(6)
where $`\eta =\rho ^+/\rho ^{}`$ is the compression ratio across the discontinuous jump, $`\gamma =5/3`$ is the ratio of specific heats of the electron gas, and superscripts minus and plus refer to the variables defined ahead of and behind the discontinuous jump, respectively 111In relation (5) we corrected the typeset errors which appeared in relation (15) of Paper II..

The mean free path of photons in the vicinity of the discontinuous jump is several orders of magnitude larger than that of hydrogen atoms and the radiation flux will remain constant to a high order of accuracy, even if we take into account the finite thickness of the viscous adiabatic jump. In our calculations, the space interval between two adjacent cell centers located ahead of and behind the discontinuity was set equal to $`1\mathrm{cm}`$ for all models. Thus, the last two terms on the r.h.s. of (5) can be omitted without loss of accuracy because they refer to the difference between the upstream and downstream values of the radiation energy density, radiation flux and radiation pressure.

Excitation and ionization of hydrogen atoms behind the discontinuity are mostly due to radiative transitions. The timescale of photoexcitation of the lower bound levels is comparable with that of the electron temperature growth, whereas the photoionization of hydrogen atoms is much slower and begins with a substantial delay. The increase in the hydrogen ionization leads to a small increase in the enthalpy flux and, correspondingly, to the small decrease in the flux of kinetic energy because the radiation flux remains almost constant. However, the most promiment process accompanying hydrogen ionization is the increase in the compression ratio $`\rho /\rho _1`$. For the model presented in Fig. 1, the compression ratio increases from $`\rho /\rho _1=3.6`$ just behind the discontinuity, to $`\rho /\rho _1=6.8`$ in the layers where the hydrogen ionization reaches its maximum.

At larger distances from the discontinuity the kinetic energy of the gas flow initially stored as the ionization energy of the gas is converted into radiation. In Fig. 1 these layers are revealed as those of remarkable changes in $`F_\mathrm{h}`$ and $`F_\mathrm{R}`$, accompanied by the growth of the compression ratio $`\rho /\rho _1`$. For models considered in the present study, the compression ratio reached at the postshock outer boundary ranges from $`1<\rho _N/\rho _110^2`$ depending mostly on the upstream velocity $`U_1`$ and less on the preshock temperature $`T_1`$ and preshock density $`\rho _1`$.

The large increase in the postshock compression ratio $`\eta =\rho /\rho _1`$ is due to the ionization of hydrogen gas and its subsequent radiative cooling. In order to roughly estimate the contribution of these effects we can write the energy conservation relation (3) with the omitted radiation energy density and radiation pressure terms as

$$\begin{array}{c}\frac{U_1^2}{2}+\stackrel{~}{E}_{\mathrm{t1}}+\stackrel{~}{E}_{\mathrm{in1}}+\frac{P_{\mathrm{g1}}}{\rho _1}+\frac{F_{\mathrm{R1}}}{\dot{m}}=\hfill \\ =\frac{U^2}{2}+\stackrel{~}{E}_\mathrm{t}+\stackrel{~}{E}_{\mathrm{in}}+\frac{P_\mathrm{g}}{\rho }+\frac{F_\mathrm{R}}{\dot{m}},\hfill \end{array}$$
(7)
where $`\stackrel{~}{E}_\mathrm{t}=\frac{3}{2}(n_\mathrm{H}+n_\mathrm{e})kT/\rho `$ is the specific energy in the translational degrees of freedom and $`\stackrel{~}{E}_{\mathrm{in}}`$ is the specific energy of excitation and ionization of hydrogen atoms. The left–hand side of (7) is written for the preshock outer boundary, whereas its right–hand side represents an arbitrary layer of the postshock region.

Relation (7) and the momentum conservation relation (2) can be rewritten as

$$\eta =4+3\frac{\stackrel{~}{E}_{\mathrm{in}}\stackrel{~}{E}_{\mathrm{in1}}}{\stackrel{~}{E}_\mathrm{t}}+3\frac{F_\mathrm{R}F_{\mathrm{R1}}}{\dot{m}\stackrel{~}{E}_\mathrm{t}},$$
(8)
where we assumed that the specific energy in the translational degrees of freedom of the preshock gas is negligible in comparison with that of the postshock compressed gas, that is, $`\stackrel{~}{E}_{\mathrm{t1}}/\eta \stackrel{~}{E}_\mathrm{t}`$.

The first term in the right–hand side of (8) is the limiting compression ratio of the atomic hydrogen gas at the discontinuous jump. The second term in the right–hand side of (8) describes the compression of the postshock gas due to excitation and ionization of hydrogen atoms. For the model shown in Fig. 1, the maximum hydrogen ionization degree is $`x_\mathrm{H}=0.99`$ and the ratio of the specific energies is $`\stackrel{~}{E}_{\mathrm{in}}/\stackrel{~}{E}_\mathrm{t}2`$. Thus, the upper limit for the compression ratio beyond the layers of hydrogen ionization is $`\eta 10`$.

However the most important contribution to the growth of the compression ratio comes from the third term in the right–hand side of (7). Indeed, for the model shown in Fig. 1 the ratio of the sum of radiative fluxes (note that $`F_{\mathrm{R1}}`$ is negative) to the flux of the translational energy of the gas at the postshock outer boundary is $`(F_{\mathrm{R}N}F_{\mathrm{R1}})/(\dot{m}\stackrel{~}{E}_N)18`$. Therefore, the maximum compression ratio at the postshock outer boundary is $`\eta 64`$.

## 3 Effects of the preshock temperature

Just behind the discontinuity, the temperature of electrons, $`T_\mathrm{e}`$, is substantially lower than the temperature of heavy particles, $`T_a`$, and electrons acquire energy in elastic collisions with neutral hydrogen atoms and hydrogen ions. Because the elastic scattering crosssection of electrons and hydrogen ions is much larger than that of electrons and neutral hydrogen atoms, the efficiency of the energy exchange between heavy particles and free electrons is very sensitive to the hydrogen ionization degree, $`x_\mathrm{H}^{}`$, ahead of the discontinuous jump. In particular, at $`x_\mathrm{H}^{}10^2`$, the postshock equilibration between the electron temperature $`T_\mathrm{e}`$ and the temperature of heavy particles $`T_a`$ is due to elastic scattering of electrons by hydrogen ions and the equilibration rate very rapidly increases with increasing $`x_\mathrm{H}^{}`$. Also, in reverse, at $`x_\mathrm{H}<10^2`$, free electrons gain most of the energy from elastic collisions with neutral hydrogen atoms and the rate of temperature equilibation is much slower. This process, however, is perceptible only in shock waves with weak ionization in the radiative precursor (that is, at upstream velocities $`U_1<30\mathrm{km}\mathrm{s}^1`$), the preshock gas temperature being $`T_1<4000\mathrm{K}`$.

Thus, effects of the preshock gas temperature on the shock wave structure are mostly due to the preshock equilibrium hydrogen ionization. This is illustrated in Fig. 2 for shock wave models with $`T_1=4000\mathrm{K}`$ and $`T_1=8000\mathrm{K}`$. In both cases, the preshock gas density and the upstream velocity are $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$ and $`U_1=40\mathrm{km}\mathrm{s}^1`$, respectively, so that the postshock equilibration of $`T_\mathrm{e}`$ and $`T_a`$ is due to elastic scattering of electrons by hydrogen ions. As is seen from the upper panel of Fig. 2, the hydrogen ionization degree just behind the discontinuous jump is $`x_\mathrm{H}0.02`$ at $`T_1=4000\mathrm{K}`$ and $`x_\mathrm{H}0.46`$ at $`T_1=8000\mathrm{K}`$. As a result, at a preshock gas temperature $`T_1=8000\mathrm{K}`$, the rate of the energy gain by electrons in elastic collisions with hydrogen ions is larger by a factor of 60 than that at $`T_1=4000\mathrm{K}`$.

The increase in the preshock gas temperature $`T_1`$ is accompanied by the stronger radiation flux emerging from the shock wave. This is, obviously, due to the higher hydrogen ionization degree in the postshock region. The more gradual changes in the postshock electron temperature and of the ionization degree in shocks with higher preshock temperature lead to smaller compression ratios $`\rho /\rho _1`$ (see the lower panel of Fig. 2).

In Fig. 3 we give the plots of the final compression ratio $`\rho _N/\rho _1`$ which is reached at the postshock outer boundary $`X_N`$ in shock wave models with $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$. It should be noted that the final postshock gas density asymptotically tends to its limiting value $`\rho _{\mathrm{}}`$ with $`X\mathrm{}`$ and because in our study the shock wave model is represented by the finite slab, the compression ratio $`\rho _N/\rho _1`$ gives only the lower estimate for the limiting value $`\rho _{\mathrm{}}/\rho _1`$. In our calculations we tried to set the postshock outer boundary as far as possible from the discontinuous jump. Unfortunately, the convergence of global iterations is very sensitive to the spatial coordinate of the postshock outer boundary and for too large values of $`X_N`$, the iterations diverge. Thus, the coordinate $`X_N`$ was determined for each model from trial calculations as a compromise between the requirement of the postshock region and demands of convergence and accuracy. For this reason, the filled circles in Fig. 3, representing the models with the same value of $`T_1`$, do not lie on smooth curves.

At a fixed upstream velocity, $`U_1`$, the Mach number $`M_1=U_1/a_1`$ decreases with increasing $`T_1`$ due to the temperature dependence of the adiabatic sound speed $`a_1`$ at the preshock outer boundary. In a hydrogen gas with negligible ionization ($`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$, $`T_17000\mathrm{K}`$) the adiabatic sound speed is $`a_1\sqrt{T_1}`$ and the compression ratio

$$\rho _{\mathrm{lim}}/\rho _1=\gamma M_1^2$$
(9)
corresponding to the isothermal shock wave is inversly proportional to $`T_1`$. On the other hand, as is seen from Fig. 3, the final compression ratio $`\rho _N/\rho _1`$ at the postshock outer boundary also decreases with increasing preshock temperature though this dependence is not as prominent as that of $`\rho _{\mathrm{lim}}/\rho _1`$.

In Fig. 4 we show the ratio $`\rho _N/\rho _{\mathrm{lim}}`$ for shock wave models with $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$, $`3000\mathrm{K}T_17000\mathrm{K}`$, $`20\mathrm{km}\mathrm{s}^1U_190\mathrm{km}\mathrm{s}^1`$. Each model on this plot is represented by a filled circle and for models with $`U_1=20\mathrm{km}\mathrm{s}^1`$ and $`U_1=90\mathrm{km}\mathrm{s}^1`$ we give the values of the upstream Mach number $`M_1`$.

At a fixed preshock gas temperature $`T_1`$ the ratio $`\rho _N/\rho _{\mathrm{lim}}`$ decreases with increasing upstream velocity since the final compression ratio $`\rho _N/\rho _1`$ grows with $`U_1`$ more slowly than $`\rho _{\mathrm{lim}}/\rho _1`$ (compare dependencies of $`\rho _N/\rho _1`$ and $`\rho _{\mathrm{lim}}/\rho _1`$ shown for $`T_1=3000\mathrm{K}`$ in Fig. 3). However the most interesting feature is that the increase in the preshock gas temperature $`T_1`$ is accompanied by the gradual approach of the final compression ratio $`\rho _N/\rho _1`$ to its upper limit $`\rho _{\mathrm{lim}}/\rho _1`$ corresponding to the isothermal shock wave. This is due to the decrease of the Mach number $`M_1`$ with increasing $`T_1`$.

## 4 Effects of the preshock density

The two upper panels of Fig. 5 show the postshock hydrogen ionization degree $`x_\mathrm{H}`$ and the postshock temperatures $`T_\mathrm{e}`$ and $`T_a`$ as a function of distance from the discontinuous jump for shock wave models with $`10^{12}\mathrm{gm}\mathrm{cm}^3\rho _110^9\mathrm{gm}\mathrm{cm}^3`$, $`T_1=6000\mathrm{K}`$ and $`U_1=50\mathrm{km}\mathrm{s}^1`$. The two lower panels of Fig. 5 for the same models show the radiation flux $`F_\mathrm{R}`$ and the compression ratio $`\rho /\rho _1`$. Note that we use the logarithmic scale for the plot of radiative flux and the deep minima of $`\mathrm{log}|F_\mathrm{R}|`$ correspond to layers with $`F_\mathrm{R}0`$.

The main results of these calculations are as follows. First, the width of the postshock relaxation zone increases with decreasing preshock gas density $`\rho _1`$. This is due to the fact that the mean collision time of particles is proportional to the gas density. Second, the radiative flux produced by the shock wave decreases with decreasing $`\rho _1`$ because the kinetic energy of the gas flow at the preshock outer boundary is proportional to $`\rho _1`$. Third, more gradual ionization and recombination of hydrogen atoms at the lower preshock gas density $`\rho _1`$ leads to smaller postshock compression ratios.

## 5 Radiative flux of the shock wave

As seen from plots on the upper panel of Fig. 1, the upstream radiative flux $`F_{\mathrm{R1}}`$ emerging from the preshock outer boundary and the downstream radiation flux $`F_{\mathrm{R}N}`$ emerging in the opposite direction from the postshock outer boundary are not equal in absolute value. This difference increases with increasing upstream velocity. For example, at $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$ and $`T_1=6000\mathrm{K}`$ the ratio of these fluxes ranges from $`|F_{\mathrm{R1}}|/F_{\mathrm{R}N}1.01`$ at $`U_1=20\mathrm{km}\mathrm{s}^1`$ to $`|F_{\mathrm{R1}}|/F_{\mathrm{R}N}1.16`$ at $`U_1=60\mathrm{km}\mathrm{s}^1`$.

However the much stronger asymmetry of the radiation field is revealed at frequencies $`\nu >\nu _1`$ because the Lyman continuum flux is directed mostly upstream. The region of the effective transport of the Lyman continuum radiation encompasses the near vicinity of the discontinuous jump. Ahead of the discontinuity, this is the radiative precursor and behind the discontinuity the Lyman continuum flux becomes negligible where the total radiation flux changes its sign (i. e.: in the layers with $`F_\mathrm{R}0`$). The contribution of the Lyman continuum to the total radiation flux within this region rapidly increases with increasing upstream velocity. For example, for shock waves propagating in a gas with $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$ and $`T_1=3000\mathrm{K}`$, the ratio of the Lyman continuum flux to the total radiation flux increases from $`F_{\mathrm{Lyc}}/F_\mathrm{R}10^2`$ at $`U_1=40\mathrm{km}\mathrm{s}^1`$ to $`F_{\mathrm{Lyc}}/F_\mathrm{R}0.2`$ at $`U_1=60\mathrm{km}\mathrm{s}^1`$. The growth of the ratio $`F_{\mathrm{Lyc}}/F_\mathrm{R}`$ ceases at an upstream velocity $`U_175\mathrm{km}\mathrm{s}^1`$, corresponding to almost full hydrogen ionization in the radiative precursor. At larger upstream velocities the fraction of the Lyman continuum radiation does not exceed $`40\%`$ of the total radiation flux, but the region of the effective transport of the Lyman continuum radiation becomes wider, both ahead of and behind the discontinuity.

The plots of the ratio $`F_{\mathrm{Lyc}}/F_\mathrm{R}`$ for shock wave models with upstream velocities $`50\mathrm{km}\mathrm{s}^1U_180\mathrm{km}\mathrm{s}^1`$ are shown in Fig. 6. It should be noted that the positive value of the ratio $`F_{\mathrm{Lyc}}/F_\mathrm{R}`$ implies that the total radiation flux and the Lyman continuum flux are both negative and are directed upstream. Negative values of this ratio imply that the Lyman continuum flux is directed downstream, whereas the total radiation flux is still upstream. Oscillations of the ratio $`F_{\mathrm{Lyc}}/F_\mathrm{R}`$ in the postshock region are due to vanishing values of $`F_\mathrm{R}`$.

Thus, behind the discontinuous jump, the Lyman continuum flux $`F_{\mathrm{Lyc}}`$ changes its sign at smaller distance than the radiative flux

$$F(\nu <\nu _1)=\underset{0}{\overset{\nu _1}{}}F_\nu 𝑑\nu ,$$
(10)
transported at frequencies below the Lyman continuum edge $`\nu _1`$. In Fig. 7 we depict the zones of the postshock region with upstream and downstream total radiative flux $`F_\mathrm{R}`$ in shock wave models with $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$ and $`T_1=6000\mathrm{K}`$. The zone with downstream Lyman continuum flux ($`F_{\mathrm{Lyc}}>0`$) and upstream flux $`F(\nu <\nu _1)`$ is shown on this diagram by the shaded area. Space coordinates $`X[F(\nu <\nu _1)=0]`$ and $`X(F_\mathrm{R}=0)`$ coincide because most of the radiation is transported in these layers at frequencies below the Lyman continuum edge.

In Fig. 7 we also show the space coordinates of the layers where half of the ionized hydrogen atoms are recombined. The close values of $`X(F_\mathrm{R}=0)`$ and $`X(\frac{1}{2}\mathrm{max}x_\mathrm{H})`$ clearly illustrate that the main mechanism of the shock wave energy dissipation is the ionization of hydrogen atoms.

At upstream velocity $`U_1<60\mathrm{km}\mathrm{s}^1`$, the postshock gas remains partially ionized and the maximum ionization degree increases with increasing $`U_1`$, while the space coordinate of the ionization maximum decreases. Correspondingly, the hydrogen recombination also occurs at a smaller distance from the discontinuous jump.

At upstream velocities of $`U_1>60\mathrm{km}\mathrm{s}^1`$, the hydrogen is fully ionized and the width of the zone with $`x_\mathrm{H}1`$ increases with increasing $`U_1`$. Correspondingly, the hydrogen recombines at larger distances from the discontinuous jump and the coordinate $`X(F_\mathrm{R}=0)`$ increases with increasing $`U_1`$. At full hydrogen ionization, the zone of the effective energy transport in the Lyman continuum spreads both upstream and downstream with increasing upstream velocity and at $`U_1>80\mathrm{km}\mathrm{s}^1`$ both the Lyman continuum flux and the total radiation flux change sign nearly in the same layers.

Another interesting conclusion which follows from the diagram in Fig. 7 is that the width of shock waves with partial postshock hydrogen ionization decreases with increasing upstream velocity $`U_1`$, whereas the width of shock waves with full postshock hydrogen ionization increases with increasing upstream velocity.

## 6 Frequency dependent radiation field

The radiation field produced by the shock wave is remarkably non–equilibrium and any attempts to use an assumption of thermal equilibrium inevitably lead to large errors. In order to evaluate the degree of departure from thermal equilibrium, it is instructive to compare the monochromatic mean intensity $`J_\nu `$ with the local Planck function $`B_\nu (T_\mathrm{e})`$. In Fig. 8 the plots of $`J_\nu `$ and $`B_\nu (T_\mathrm{e})`$ are shown as a function of frequency $`\nu `$ for four distinct layers of the shock wave model with $`\rho =10^{10}\mathrm{gm}\mathrm{cm}^3`$, $`T_1=6000\mathrm{K}`$ and $`U_1=60\mathrm{km}\mathrm{s}^1`$. For the sake of convenience each pair of plots is shifted with respect to others.

The upper pair of plots represents $`J_\nu `$ and $`B_\nu (T_\mathrm{e})`$ at $`X10^4\mathrm{cm}`$. As is seen from Fig. 1, the photoionization of hydrogen atoms in these layers due to absorption of the Lyman continuum radiation becomes perceptible and this layer can roughly be considered as the boundary of the radiative precursor. Excess radiation in comparison with $`B_\nu (T_\mathrm{e})`$ at frequencies $`\nu >\nu _1`$ is due to the fact that the optical depth of the radiative precursor decreases with increasing frequency $`\nu `$.

The second and third pairs of plots represent the cell just behind the discontinuous jump ($`X=0.5\mathrm{cm}`$) and the layer where the hydrogen ionization degree reaches its maximum ($`X40\mathrm{cm}`$). Between these layers, various relaxation processes redistribute the energy of heavy particles among other degrees of freedom and the radiation field is furthest from equilibrium.

The lower pair of plots represent the layers ($`X1.210^5\mathrm{cm}`$) where most of the hydrogen atoms are recombined and the gas temperature gradually approaches its unperturbed value $`T_1`$. The optical depth at frequencies of the Lyman continuum is so large that $`J_\nu `$ and $`B_\nu (T_\mathrm{e})`$ coincide at $`\nu >\nu _1`$.

All shock wave models demonstrate two common features. First, we see as excess of Lyman continuum radiation in the preshock region, due to the existence of the radiative precursor. Second, there is a lower mean intensity of $`J_\nu `$ in comparison with $`B_\nu (T_\mathrm{e})`$, due to the small optical depth of the shock wave at frequencies $`\nu <\nu _1`$.

The plots of the mean intensity $`J_\nu `$ displayed in Fig. 8 reveal the presence of six spectral line features. These are Ly$`\alpha `$, Ly$`\beta `$, Ly$`\gamma `$, H$`\alpha `$, H$`\beta `$ and Pa$`\alpha `$. In comparison with our previous work described in Paper II, we improved the treatment of the spectral line radiation transfer and considered each line of frequency $`\nu _0`$ within the frequency interval as wide as $`[\nu _0(1\delta ),\nu _0(1+\delta )]`$, where $`\delta =510^4`$. The use of wide frequency intervals is necessary because of the strong Doppler broadening of line profiles behind the discontinuous jump. In particular, the radiation flux transported in spectral lines is not negligible and contribution of the spectral line radiation to the total radiation flux increases with increasing upstream velocity. The most important are the Ly$`\alpha `$ and H$`\alpha `$ lines. For example, in the shock wave model with $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$, $`T_1=6000\mathrm{K}`$ and $`U_1=60\mathrm{km}\mathrm{s}^1`$ the upstream radiation flux in Ly$`\alpha `$ and H$`\alpha `$ lines emerging from the preshock outer boundary is $`2\%`$ of the total radiation flux.

In Fig. 9 we show the plots of the radiation flux integrated within frequency intervals of the Ly$`\alpha `$, Ly$`\beta `$, Ly$`\gamma `$ (upper panel), H$`\alpha `$, H$`\beta `$ and Pa$`\alpha `$ (lower panel) lines. Comparing with plots of $`n_2/n_\mathrm{H}`$ and $`x_\mathrm{H}`$ displayed in Fig. 1 we can conclude that the spectral line radiation is produced mostly in the layers of hydrogen recombination. Another interesting feature is that the upstream and downstream Ly$`\alpha `$ radiation fluxes are nearly equal in absolute value, whereas the radiation flux in Balmer lines emerges mostly upstream.

## 7 Radiative losses of the shock wave

In order to estimate the irreversible loss occurring in the shock wave energy it is instructive to compare the values of each term of the left–hand side of relation (3) at both boundaries of the slab. In Fig. 10 these quantities are shown for shock wave models with $`\rho _1=10^{10}\mathrm{gm}\mathrm{cm}^3`$, $`T_1=6000\mathrm{K}`$ and $`20\mathrm{km}\mathrm{s}^1U_180\mathrm{km}\mathrm{s}^1`$. As is seen by comparing the upper and lower panels representing the preshock and postshock outer boundaries, respectively, we see that most of the kinetic energy of the gas flow is irreversibly lost, since the ratio of the fluxes of kinetic energy at both boundaries of the slab ranges within $`10^2F_{\mathrm{K1}}/F_{\mathrm{K}N}710^3`$. At the same time, the increase of the enthalpy flux after the passage of the parcel of gas through the slab does not exceed 20%. Thus, most of the energy of the shock wave is lost due to radiation.

The role of radiation in irreversible energy losses of the shock wave can be evaluated from comparison of the radiation flux emerging from the boundary of the slab with the total energy flux $`𝒞_2`$. In Fig. 11 the plots of the ratio $`F_{\mathrm{R}N}/𝒞_2`$ at the postshock outer boundary are shown for two sequences of shock wave models with $`T_1=3000\mathrm{K}`$ and $`T_1=6000\mathrm{K}`$. As is seen, the radiative losses increase very rapidly with increasing upstream velocity $`U_1`$, though there is also the dependence on the preshock gas temperature $`T_1`$.

## 8 Conclusion

In this work we improved the treatment of radiation transfer and considered the structure of shock waves propagating through partially ionized hydrogen gas with densities from $`10^{12}\mathrm{gm}\mathrm{cm}^3`$ to $`10^9\mathrm{gm}\mathrm{cm}^3`$. Our calculations showed that while the unpertubed hydrogen gas is partially ionized, the preshock gas temperature affects the shock wave structure, mostly due to the preshock hydrogen ionization because the postshock equilibration of translational degrees of freedom is very sensitive to the number density of hydrogen ions. In general, a higher temperature in the preshock gas will lead to larger postshock ionization and to a stronger radiation field produced by the shock wave. At the same time, the increase in the preshock gas temperature is accompanied by the slight decrease of the postshock gas compression.

The preshock gas density $`\rho _1`$ affects the shock wave structure mostly due to the fact that the preshock flux of kinetic energy as well as the radiation flux emerging from the shock wave are proportional to $`\rho _1`$. Therefore, in the gas with lower density, the shock wave produces weaker radiation flux and less compression of the postshock gas. Moreover, at lower gas density, the width of the shock wave increases. The width of the radiative precursor increases due to the smaller absorption in the Lyman continuum, whereas the width of the postshock relaxation zone increases because of the proportionality of the rates of relaxation processes to the gas density.

The increase in the shock wave velocity leads to higher rates of relaxation processes in the compressed gas behind the discontinuous jump and to higher postshock ionization. While the maximum ionization degree of the postshock gas is less than unity, the width of the shock wave decreases with increasing shock wave velocity. However, when the postshock gas is fully ionized, the increase in the shock wave velocity leads to the extension of the postshock zone of hydrogen ionization and, therefore, to an increase in the shock wave width.

Strong radiative cooling of the postshock gas leads to compression ratios as high as $`\rho /\rho _110^2`$. Such large increases in the postshock gas density obviously can favour the condensation of dust grains in outer atmospheres of radially pulsating late–type giants.

The isothermal approximation undoubtedly overestimates the radiative losses of the shock wave energy by at least a factor of two for shocks with an upstream velocity of $`U_130\mathrm{km}\mathrm{s}^1`$ and can provide good accuracy at upstream velocities of $`U_1>80\mathrm{km}\mathrm{s}^1`$, that is at Mach numbers $`M_110`$.

Our calculations demonstrated that the monochromatic radiation flux at frequencies of hydrogen lines significantly exceeds the flux of the background continuum, so that our models reproduce the emission spectrum observed in radially pulsating late–type supergiants. Moreover, the fraction of the energy transported in hydrogen lines is not negligible and can be as large as one or two percent of the total radiation flux because of the strong broadening of the line profiles behind the discontinuous jump.

In this paper, we did not discuss the details of the emission line spectra because current results are insufficient for comparison with observations. The change in the velocity of the gas during its passage through the shock wave must inevitably influence the line profiles. Calculations of shock wave structure based on the solution of the transfer equation which takes into account effects of the gas flow velocity gradient were beyond the scope of the present work and will be presented in the forthcoming paper.

###### Acknowledgements.

The work of YAF was done in part under the auspices of the Aix-Marseille I University and CNRS. YAF acknowledges also the support from the Russian Foundation for Basic Research (grant 98–02–16734).


text:
# Spin and Current Correlation Functions in the 𝐝-density Wave State of the Cuprates

## I Introduction

Recently, a new order termed $`d`$-density wave (DDW) was proposed for the underdoped regime of the high-$`T_c`$ cuprate superconductors. It was argued that the observed anomalies in the pseudogap gap phase are due to the DDW and its competition with the $`d`$-wave superconductor (DSC). Recent elastic neutron scattering , $`\mu `$SR , and NMR/NQR experiments appear to have found some direct evidence for the existence of DDW order in YBCO by observing the magnetic fields associated with this order parameter – in other words, by measuring the order parameter itself. In this paper, we wish to elucidate the nature of the competition between the DDW and DSC order paramters by explicitly computing the spin-spin and the current-current correlation functions in two space-dimensions in the Hartree-Fock approximation. Our aim is to provide complementary, indirect evidence in favor of the DDW scenario by showing in some detail that it provides a consistent explanation of the existence and behavior of the resonant peak seen in inelastic neutron scattering and also of the evolution of the superfluid density with doping and temperature, at low temperatures.

Since we deal with broken symmetry states with order parameters, it is possible to argue that their descriptions are adequately described by the Hartree-Fock approximation, at least deep within the ordered phases at zero temperature. The Hartree-Fock approximation by itself cannot be adequate if collective modes are important, but these modes could be described in terms of small oscillations about the ordered state by suitable random phase approximations (RPA). However, critical fluctuations close to a quantum critical point cannot be treated by either Hartree-Fock or RPA. Nonetheless, it is useful to explore the consequences of the simplest possible Hartree-Fock and RPA theories to see if any robust features arise due to the existence of the new order parameter, DDW, and its competition with DSC.

For underdoped bilayer superconductors , the dynamic spin-spin correlation function is peaked above $`T_c`$ at the in-plane wavevector $`𝐐=(\pi ,\pi )`$ (in this paper we shall take the lattice constant to be unity and will also set $`\mathrm{}=c=k_B=1`$) in both even and odd channels with respect to the layers in the bilayer complex. Below $`T_c`$, there is a resonant peak at the same wavevector and approximately the same energy in the odd channel. A resonant peak is also found below $`T_c`$ in optimal and overdoped bilayer high-$`T_c`$ superconductors. There are many approximate theoretical calculations of the dynamic structure factor. Here we concentrate only on the aspects germane to our discussion of DDW. We wish to stress that this analysis is simply a consistency check of the DDW; we do not claim that our explanation of the resonant peak in inelastic neutron scattering is unique or better than the others. The real test of the DDW hypothesis is the direct observation of the order parameter referred to earlier. We find that the experimental observations of the neutron dynamic structure factor are consistent with the hypothesis that the neutron scattering peak – and, hence, the pseudogap – is due to DDW order, while the resonance peak in the underdoped superconducting state is due to the coexistence of DDW and DSC orders. The suppression of the resonant peak by a magnetic field perpendicular to the plane suggests that superconducting pairing is important for its formation.

The behavior of the superfluid density in the underdoped regime of the high-$`T_c`$ materials is another interesting physical quantity. Many experiments have indicated a rapid collapse of the zero-temperature superfluid density as the doping is decreased below optimal. The zero-temperature penetration depth, for example, grows rapidly in the pseudogap regime, correlating with the suppression of $`T_c`$ with underdoping, yet saturates at overdoping in a way reminiscent of a traditional BCS superconductor. This, as we will show, can be explained by the competition between DDW and DSC orders. At finite temperatures the superfluid density is suppressed from its zero temperature value at that doping fraction. This suppression is linear in temperature ($`T`$) at asymptotically low temperatures, with a slope that is independent of doping. This behavior is captured in our calculations by the thermal excitation of nodal quasiparticles into the upper quasiparticle band. For heavily underdoped samples, there is an intermediate low temperature regime in which we predict this suppression to be proportional to $`\sqrt{T}`$. There is some experimental evidence for this behaviorin YBa2Cu4O8, although the explanation of these authors involves a proximity model of alternating stacked superconducting and normal layers. Our analysis is simpler and follows from the nodal exciations in the mixed DDW and DSC state. This prediction can be tested in future experiments in the regime in which pseudogap dominates the DSC gap.

## II Order Parameter

Our Hartree-Fock analysis of the physical quantities merely requires us to specify a mean field Hamiltonian with the proposed order parameters. The actual microscopic Hartree-Fock analyses to obtain these order parameters have been discussed in the past on numerous occasions and need not be repeated here; for a recent set of calculations, see Refs. , , and references therein. Hamiltonians with short range repulsion and superexchange have the DDW ordered state (both the singlet and the triplet variety) as one of their many possible saddle points. If one includes correlated hopping terms, then within a reasonable range of parameters, the DDW saddle point can be stabilized over the other possibilities. In the present paper, we shall parametrize the order parameters phenomenologically, and our main conclusions are independent of any microscopic self-consistent Hartree-Fock equations.

The singlet DDW state is defined by the order parameter

$$c_{𝐤\alpha }^{}c_{𝐤+𝐪\beta }=i\frac{\mathrm{\Phi }_𝐪}{2}(\mathrm{cos}k_x\mathrm{cos}k_y)\delta _{\alpha \beta },$$
(1)
where $`c_{𝐤\alpha }`$ is the fermion destruction operator for wavevector $`𝐤`$ and spin-index $`\alpha `$, and $`𝐪`$ is the ordering wavevector. This order parameter is a generalization of the familiar charge density wave (CDW) order parameter to the case of angular momentum two, in exact analogy to the generalization of the BCS $`s`$-wave superconducting order parameter to its $`d`$-wave (DSC) counterpart.

For $`𝐪=𝐐`$, the case of interest to us, the underlying bipartite square lattice band structure is equivalent under the transformation $`𝐐𝐐`$ and this forces the DDW order parameter to be imaginary . Thus, time-reversal symmetry is broken (i.e. the system exhibits magnetism), and the ground state has an array of bond currents, which alternate in direction (clockwise, anticlockwise) in the neighboring plaquettes of the 2D bipartite square lattice. The corresponding $`s`$-wave CDW order parameter is, of course, real.

## III Spin Dynamics

### A Pure DDW state

The Hartree-Fock DDW Hamiltonian, using (1), is

$$H_{\mathrm{DDW}}=\underset{𝐤\sigma }{}(ϵ_𝐤\mu )c_{𝐤\sigma }^{}c_{𝐤\sigma }+\underset{𝐤\sigma }{}iW_𝐤c_{𝐤\sigma }^{}c_{𝐤+𝐐\sigma }+\text{h.c.},$$
(2)
where the DDW gap is given by

$$W_𝐤=\frac{W_0}{2}(\mathrm{cos}k_x\mathrm{cos}k_y),$$
(3)
$`ϵ_𝐤=2t(\mathrm{cos}k_x+\mathrm{cos}k_y)`$ gives the band structure, and $`\mu `$ is the chemical potential. A more realistic band structure would include the effect of next-neighbor hopping $`t^{}`$. This does not affect our results here, so we drop it for simplicity. At half filling, the chemical potential $`\mu =0`$, while $`\mu `$ takes non-zero negative values as we introduce holes into the system.

At $`\mu =0`$ the zero-temperature spin-spin correlation function $`𝒮`$, at momentum transfer $`𝐪=𝐐`$, is given by

$$𝒮(𝐐,\omega ,\mu =0)=3\pi \underset{𝐤\mathrm{rbz}}{}\delta \left(\omega 2\sqrt{W_𝐤^2+ϵ_𝐤^2}\right),$$
(4)
where the integration is over the magnetic or reduced Brillouin zone (rbz). This function is peaked at an energy equal to twice the maximum value of the gap, $`2W_0`$. As we dope holes into the system, at nonzero values of $`\mu `$, Eq. 4 is changed to

$`𝒮(𝐐,\omega ,\mu )=3\pi {\displaystyle \underset{𝐤\mathrm{rbz}}{}}`$ $`\delta `$ $`\left(\omega 2\sqrt{W_𝐤^2+ϵ_𝐤^2}\right)`$ (5)
$`\times `$ $`\theta \left(\mu +\sqrt{W_𝐤^2+ϵ_𝐤^2}\right)`$ (6)
The only effect of $`\mu `$ is the depletion of the integration region. Hence, as long as $`|\mu |<W_0`$, the peak at $`𝐪=𝐐`$ exists and stays at the same energy $`2W_0`$.

Note that $`𝒮(\omega )`$ also shows a peak for momentum transfer $`𝐪=(\pi ,0)`$  at the Hartree-Fock level. From a relation analogous to that of Eq. 6, we can see that the peak-energy for $`𝐪=(\pi ,0)`$ scales with $`t`$. This is much higher in energy than the $`(\pi ,\pi )`$ peak; it is likely to have an enormous width (since there is an enormous amount of phase space into which it can decay) and be unobservable. This is important because no peak has been observed at $`𝐪=(\pi ,0)`$ in the spin-fluctuation spectrum of the pseudogap regime of the high-$`T_c`$ materials. Early mean-field decouplings of the Heisenberg model  led to an effective $`tJ`$, so that the peaks at $`𝐪=(\pi ,\pi )`$ and $`𝐪=(\pi ,0)`$ were comparable in energy, in contradiction to the experiments.

### B Pure DSC state

Below the superconducting transition temperature $`T_c`$, a DSC state is defined by the gap

$$\mathrm{\Delta }_𝐤=\frac{\mathrm{\Delta }_0}{2}(\mathrm{cos}k_x\mathrm{cos}k_y)$$
(7)
corresponding to the DSC order parameter. This order parameter, as stressed in reference, can compete and coexist with the singlet DDW order parameter. In fact, at half filling, the system can be rotated continuously from a pure DSC order to a pure singlet DDW order and vice versa, without ever having to close the quasi-particle gap . (The symmetry between the $`s`$-wave counterparts of these two types of orders is exact in the negative $`U`$ Hubbard model .)

Using (6), the mean-field one-body DSC Hamiltonian is

$$H_{\mathrm{DSC}}=\underset{𝐤\sigma }{}(ϵ_𝐤\mu )c_{𝐤\sigma }^{}c_{𝐤\sigma }+\underset{𝐤}{}\mathrm{\Delta }_𝐤c_𝐤^{}c_𝐤^{}+\text{h.c.}$$
(8)
The spin-spin correlation function, at $`\mu =0`$ and momentum transfer $`𝐪`$, is

$`𝒮(𝐪,\omega ,\mu =0)`$ $`=`$ $`{\displaystyle \frac{3\pi }{2}}{\displaystyle \underset{𝐤\mathrm{rbz}}{}}\left(1{\displaystyle \frac{ϵ_𝐤ϵ_{𝐤+𝐪}+\mathrm{\Delta }_𝐤\mathrm{\Delta }_{𝐤+𝐪}}{E_𝐤E_{𝐤+𝐪}}}\right)`$ (9)
$`\times `$ $`\delta (\omega E_𝐤E_{𝐤+𝐪}),`$ (10)
where $`E_𝐤=\sqrt{\mathrm{\Delta }_𝐤^2+ϵ_𝐤^2}`$. For $`𝐪=𝐐`$ the coherence factor equals two and Eq. 10 becomes identical to Eq. 4. Hence, just like the DDW, at $`\mu =0`$ the DSC spin-spin correlator also shows a peak at twice the gap maximum, $`2\mathrm{\Delta }_0`$.

But for a DSC, at non-zero values of $`\mu `$, $`𝒮(𝐪,\omega ,\mu )`$ is given by

$`𝒮(𝐪,\omega ,\mu )`$ $`=`$ $`{\displaystyle \frac{3\pi }{2}}{\displaystyle \underset{𝐤\mathrm{rbz}}{}}\left(1{\displaystyle \frac{\stackrel{~}{ϵ}_𝐤\stackrel{~}{ϵ}_{𝐤+𝐪}+\mathrm{\Delta }_𝐤\mathrm{\Delta }_{𝐤+𝐪}}{E_𝐤E_{𝐤+𝐪}}}\right)`$ (11)
$`\times `$ $`\delta (\omega E_𝐤E_{𝐤+𝐪}),`$ (12)
where $`E_𝐤`$ now incorporates $`\mu `$ in the usual way by $`ϵ_𝐤\stackrel{~}{ϵ}_𝐤=ϵ_𝐤\mu `$. Equation 12 shows a peak at $`𝐪=𝐐`$ that is shifted in energy from $`2\mathrm{\Delta }_0`$. The peak-energy shift is small for small $`|\mu |`$ , but can be easily checked to scale with $`2|\mu |`$ as $`|\mu |`$ becomes large. Hence, one can clearly see that the peak-energy behaves differently with $`\mu `$ for DDW and DSC.

These Hartree-Fock calculations for the resonance peak in the normal state becom e interesting only if the chemical potential remains pinned close to $`\mu =0`$ in the underdoped regime. (This is not true for the calculation of the superfluid density reported in the following section, which is insensitive to such assumptions for $`\mu `$.) There is some evidence that this is the case in the photoemission measurements of Ino et al. in $`\text{La}_{2x}\text{Sr}_x\text{CuO}_4`$. They find that while the chemical potential shift is large in the overdoped samples, it is largely suppressed in the underdoped regime. Numerical studies of the 2D Hubbard model also suggest that the chemical potential does not shift much for small doping-fractions. In the Monte Carlo study performed in reference , the calculated shift of $`\mu `$ follows $`\mathrm{\Delta }\mu =\mu x^2`$ for small values of doping, $`x`$. The data of reference suggest that the chemical potential varies with $`x`$ in the overdoped regime as one would expect for a Fermi liquid, but not in the underdoped regime. Ideas involving charge ordered stripe states are suggestive of the segregation of doped holes in the boundaries of antiferromagnetic domains, thus pinning the chemical potential. The presence of charge inhomogeneity, such as stripes, or impurity-disorder may pin the chemical potential in other multilayer high-$`T_c`$ cuprate systems as well.

### C Coexisting DDW and DSC

From the order parameter competition picture proposed in reference (See, Fig. 1), one notes that to the right of $`x0.05`$ and below $`T_c`$, the two orders, singlet DDW and DSC, coexist and compete for the same regions on the Fermi surface. The mean-field Hamiltonian for the system in the folded Nambu basis is

$$H\mu N=\underset{𝐤\mathrm{rbz}}{}\mathrm{\Psi }_𝐤^{}A_𝐤\mathrm{\Psi }_𝐤,$$
(13)
where $`\mathrm{\Psi }_𝐤^{}=(c_𝐤^{},c_{𝐤+𝐐}^{},c_𝐤,c_{𝐤𝐐})`$ and the matrix $`A_𝐤`$ is

$$A_𝐤=\left(\begin{array}{cccc}(ϵ_𝐤\mu )& iW_𝐤& \mathrm{\Delta }_𝐤& 0\\ iW_𝐤& (ϵ_𝐤+\mu )& 0& \mathrm{\Delta }_𝐤\\ \mathrm{\Delta }_𝐤& 0& (ϵ_𝐤\mu )& iW_𝐤\\ 0& \mathrm{\Delta }_𝐤& iW_𝐤& (ϵ_𝐤+\mu )\end{array}\right)$$
(14)
The matrix $`A_𝐤`$ has four eigenvalues, $`\pm E_1(𝐤)`$ and $`\pm E_2(𝐤)`$, where $`E_1(𝐤)=[(\sqrt{ϵ_𝐤^2+W_𝐤^2}\mu )^2+\mathrm{\Delta }_𝐤^2]^{1/2}`$ and $`E_2(𝐤)=[(\sqrt{ϵ_𝐤^2+W_𝐤^2}+\mu )^2+\mathrm{\Delta }_𝐤^2]^{1/2}`$.

The spin-spin correlation function is obtained by first evaluating the imaginary-time ordered correlator ($`\omega _m`$ is the bosonic Matsubara frequency)

$`T_\tau `$ $`𝐒(𝐪,i\omega _m)𝐒(𝐪,i\omega _m)=`$ (16)
$`{\displaystyle \frac{3}{4\beta }}{\displaystyle \underset{n,𝐤\mathrm{rbz}}{}}\mathrm{Tr}\left[G(𝐤,i\omega _n)G(𝐤+𝐪,i\omega _n+i\omega _m)\right]`$
where $`𝐒`$ denotes the spin operator and the $`G`$’s are the 4$`\times `$4 matrix Green functions computed from the Hamiltonian in Eq. 13. After analytically continuing the frequency to the real axis, we extract the imaginary part and use the fluctuation-dissipation theorem to get the spin-spin correlation function.

The zero-temperature spin-spin correlation function thus obtained, at non-zero $`\mu `$, for momentum transfer $`𝐪=𝐐`$, is given by

$`𝒮(𝐐,\omega ,\mu )={\displaystyle \frac{3\pi }{2}}{\displaystyle \underset{𝐤\mathrm{rbz}}{}}`$ $`\left(1+{\displaystyle \frac{ϵ_𝐤^2+\mathrm{\Delta }_𝐤^2+W_𝐤^2\mu ^2}{E_1(𝐤)E_2(𝐤)}}\right)`$ (17)
$`\times `$ $`\delta \left(\omega E_1(𝐤)E_2(𝐤)\right)`$ (18)
When $`W_𝐤=0`$, the coherence factor in Eq. 18 clearly matches with the coherence factor in Eq. 12, where we have to put $`𝐪=𝐐`$, and the two expressions become identical. On the other hand, for $`\mathrm{\Delta }_𝐤=0`$, one notes that $`E_1(𝐤)|\sqrt{ϵ_𝐤^2+W_𝐤^2}\mu |`$, and $`E_2(𝐤)|\sqrt{ϵ_𝐤^2+W_𝐤^2}+\mu |`$. Keeping in mind that $`\mu `$ is negative in our hole-doped system, we note that for $`\sqrt{ϵ_𝐤^2+W_𝐤^2}+\mu >0`$ the coherence factor in Eq. 18 is two, while the other choice for $`\sqrt{ϵ_𝐤^2+W_𝐤^2}+\mu <0`$ reduces the coherence factor to zero. Thus the expression in Eq. 18 becomes identical to the expression in Eq. 6 for $`\mathrm{\Delta }_𝐤=0`$.

One can easily check that $`𝒮(𝐐,\omega ,\mu )`$, as a function of $`\omega `$, has a peak in the spin-spin correlation function . The peak is located at $`2\sqrt{W_0^2+\mathrm{\Delta }_0^2}`$ for $`\mu =0`$. For small $`\mu `$ the peak-energy is shifted to higher values, initially by a small amount, but finally scaling as $`2|\mu |`$ for large values of $`|\mu |`$. We should mention here that one can in fact control the behavior of the peak-energy with $`\mu `$ by adjusting the relative strengths of the order parameters. When the singlet DDW order overshadows the DSC order, the peak-energy will tend to be pinned to the total order parameter magnitude $`2\sqrt{W_0^2+\mathrm{\Delta }_0^2}`$, a behavior characteristic of the singlet DDW.

### D RPA Corrections

If we add a reduced Coulomb repulsion $`\overline{U}`$, which is assumed to be renormalized due to particle-particle correlations, we can drive the system toward antiferromagnetism. To model this effect, we have to go beyond the Hartree-Fock, and we use a crude RPA form for the susceptibility $`\chi (𝐪,\omega )`$ given by,

$$\chi (𝐪,\omega )=\frac{\chi _0(𝐪,\omega )}{1\overline{U}\chi _0(𝐪,\omega )}$$
(19)
to describe the spin dynamics in the system. Here $`\chi _0(𝐪,\omega )`$ is the Hartree-Fock susceptibility. Extracting $`𝒮(𝐪,\omega ,\mu )`$ from $`\chi (𝐪,\omega ,\mu )`$, in Fig. 2 we plot the results as a function of frequency for five different values of $`\overline{U}`$.

The first thing one notices is that after the introduction of $`\overline{U}`$, the correlation function $`𝒮`$, which was logarithmically singular when derived from $`\chi _0(𝐪,\omega )`$ alone, now becomes broad. Also, as $`\overline{U}`$ is progressively increased, the spectral-weight appears at smaller and smaller energies. For $`\overline{U}`$ large enough, the correlator starts peaking up again and finally evolves into a divergence at zero energy as $`\overline{U}_c`$ is approached. All of this behavior is qualitatively consistent with what happens to the spin-fluctuation spectrum as the doping is reduced in the underdoped regime of the cuprates. As the hole density is reduced, the repulsive Coulomb energy becomes more and more important and gradually approaches the critical value $`\overline{U}_c`$ for which the peak in the spin-spin correlation function evolves into a divergence at zero energy, signalling antiferromagnetic order in the system.

This analysis applies both below and above $`T_c`$ in the high-$`T_c`$ cuprate materials. Above the superconducting transition temperature $`T_c`$, $`\mathrm{\Delta }_𝐤`$ will go to zero, but $`W_𝐤`$ will survive up to the pseudogap temperature $`T^{}`$. Hence, in the spin-spin correlator the DDW will continue to maintain its signature intact in the form of a peak at twice the maximum of its gap.

## IV Superfluid Density

The anomalous behavior of the superfluid density in the underdoped regime of the high-$`T_c`$ cuprates has been discussed extensively. Here, we address it in the context of the competition of the order parameters that we have been discussing. We will make a phenomenological assumption about the doping dependence of the DDW and single particle gaps, respectively, $`W_𝐤`$ and $`\sqrt{\mathrm{\Delta }_𝐤^2+W_𝐤^2}`$. Following , we hypothesize that the DDW order parameter vanishes above a critical doping $`x_c`$ (about 0.2) and follows the formula $`W_0/k_B\alpha J^{}(1x/x_c)`$ where $`J^{}=980K`$, $`x`$ is the doping-fraction, and $`\alpha 0.6`$ is a constant arbitrarily chosen to set a scale. In other words,

$$W_0.049(1x/x_c)$$
(20)
in eV; in Fig. 3, we test another model in which $`W_0(1x/x_c)^2`$.

We also posit that the single-particle gap remains constant at its value at $`x=.05`$, where superconductivity starts developing; then

$$W_0^2+\mathrm{\Delta }_0^2\mathrm{\Phi }_0^2=(.037)^2$$
(21)
We calculate the superfluid stiffness for such a model using the formula

$$\frac{n_s}{m^{}}=k_x\mathrm{\Lambda }_{xx}(q_x=0,q_y0,i\omega _m=0)$$
(22)
which gives the response to a transverse vector potential $`A_x(t)`$. Here, $`n_s`$ is the superfluid density and $`m^{}`$ is the effective mass, $`k_x`$ is the kinetic energy per site per lattice dimension, and $`\mathrm{\Lambda }_{xx}(𝐪,i\omega _m)`$ is the paramagnetic current-current correlation function.

The paramagnetic current response is obtained from the relation

$`\mathrm{\Lambda }_{xx}(𝐪,i\omega _m)`$ $`=`$ $`T_\tau j_x(𝐪,i\omega _m)j_x(𝐪,i\omega _m)`$ (23)
$`=`$ $`{\displaystyle \frac{4t^2}{\beta }}{\displaystyle \underset{n,𝐤\mathrm{rbz}}{}}\mathrm{sin}k_x\mathrm{sin}(𝐤+𝐪)_x`$ (27)
$`\mathrm{Tr}[G(𝐤,i\omega _n)MG(𝐤+𝐪,i\omega _n+i\omega _m)M]`$
$`{\displaystyle \frac{4W_0^2}{\beta }}{\displaystyle \underset{n,𝐤\mathrm{rbz}}{}}\mathrm{sin}k_x\mathrm{sin}(𝐤+𝐪)_x`$
$`\mathrm{Tr}[G(𝐤,i\omega _n)NG(𝐤+𝐪,i\omega _n+i\omega _m)N]`$
In Eq. 27, the $`G`$’s are the same matrix Green functions used in the calculation of the spin-spin correlation function and $`M`$ and $`N`$ are the matrices

$`M`$ $`=`$ $`\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right),`$ (32)
$`N`$ $`=`$ $`\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right)`$ (37)
The superfluid stiffness (and therefore the inverse square penetration depth) at $`T=0`$, evaluated in the entire doping range, is given by

$`{\displaystyle \frac{n_s}{m^{}}}`$ $`=`$ $`32t^2{\displaystyle \underset{𝐤\mathrm{rbz}}{}}\mathrm{sin}^2k_xW_𝐤^2\left(1+{\displaystyle \frac{ϵ_𝐤^2+\mathrm{\Delta }_𝐤^2+W_𝐤^2\mu ^2}{E_1(𝐤)E_2(𝐤)}}\right)`$ (38)
$`\times `$ $`{\displaystyle \frac{1}{[E_1(𝐤)+E_2(𝐤)]^3}}(1+𝒪(W_0^2/t^2))`$ (39)
$`+`$ $`8t^2\mathrm{lim}_\beta \mathrm{}{\displaystyle \underset{𝐤\mathrm{rbz}}{}}\mathrm{sin}^2k_x{\displaystyle \frac{ϵ_𝐤^2}{W_𝐤^2+ϵ_𝐤^2}}`$ (40)
$`\times `$ $`\left({\displaystyle \frac{df(E_1(𝐤))}{dE_1(𝐤)}}+{\displaystyle \frac{df(E_2(𝐤))}{dE_2(𝐤)}}\right)(1+𝒪(W_0^2/t^2))`$ (41)
$``$ $`4t{\displaystyle \underset{𝐤\mathrm{rbz}}{}}\mathrm{cos}k_x\left(1+{\displaystyle \frac{ϵ_𝐤^2+\mathrm{\Delta }_𝐤^2+W_𝐤^2\mu ^2}{E_1(𝐤)E_2(𝐤)}}\right)`$ (42)
$`\times `$ $`{\displaystyle \frac{ϵ_𝐤}{E_1(𝐤)+E_2(𝐤)}}(1+𝒪(W_0^2/t^2))`$ (43)
where $`E_1(𝐤)`$ and $`E_2(𝐤)`$ are the two energy-values given before. The first two terms in Eq. 43 come from the paramagnetic current response, and the last term is the kinetic energy (or the diamagnetic term). At finite, but low temperatures, one can show that

$`{\displaystyle \frac{n_s}{m^{}}}(T)`$ $``$ $`{\displaystyle \frac{n_s}{m^{}}}(0)8t^2{\displaystyle \underset{𝐤\mathrm{rbz}}{}}\mathrm{sin}^2k_x{\displaystyle \frac{ϵ_𝐤^2}{W_𝐤^2+ϵ_𝐤^2}}`$ (44)
$`\times `$ $`\left({\displaystyle \frac{df(E_1(𝐤))}{dE_1(𝐤)}}+{\displaystyle \frac{df(E_2(𝐤))}{dE_2(𝐤)}}\right)(1+𝒪(W_0^2/t^2)).`$ (45)
The $`𝒪(W_0^2/t^2)`$ terms derive from the modification of both the kinetic energy and current operators in the low-energy effective (Hartree-Fock) Hamiltonian associated with DDW order. The modified current operator gives rise to the matrix $`N`$ in (27). Since these terms are small, $`𝒪(W_0^2/t^2)`$, we neglect them in the following discussion.

Equation 43 has several interesting properties. First, for $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_00`$ and $`W_𝐤=W_0`$, that is, for competing CDW and $`s`$-wave superconductivity, the terms involving the Fermi functions do not contribute at zero-temperature. At $`\mu =0`$, $`E_1(𝐤)=E_2(𝐤)`$, and the same coherence factor in the other two terms becomes two. Then, by a partial integration one can show that (denoting the first term by $`K^{\mathrm{para}}`$ and the third term by $`K^{\mathrm{dia}}`$), $`K^{\mathrm{para}}=\frac{W_0^2}{\mathrm{\Phi }_0^2}K^{\mathrm{dia}}`$, and the full kernel is given by $`K^{\mathrm{total}}=K^{\mathrm{para}}+K^{\mathrm{dia}}=K^{\mathrm{dia}}(1\frac{W_0^2}{\mathrm{\Phi }_0^2})`$. When the two orders are of $`s`$-wave type , $`K^{\mathrm{total}}`$ yields a superfluid density which is maximum when $`W_0`$=0 (i.e, $`\mathrm{\Delta }_0=\mathrm{\Phi }_0`$) and zero when $`W_0`$=$`\mathrm{\Phi }_0`$ (i.e, $`\mathrm{\Delta }_0`$=0).

Returning to $`d`$-wave order parameters, at $`\mathrm{\Delta }_0=0`$, Eq. 43 reduces to the superfluid density of the DDW state, which is, of course, zero. The second term, which contains the derivatives of the Fermi functions is crucial for the cancellation in this case. For $`\mathrm{\Delta }_0`$ and $`\mu `$ finite, the second term does not contribute at zero temperature (except for a contribution from a single point in k-space; but the finite and almost constant contribution from this is ignored here for simplicity). The system now acquires a finite superfluid density.

The superfluid density derived from Eq. 43 can be cast as a function of a single variable, $`x`$, by expressing the three parameters $`W_0`$, $`\mathrm{\Delta }_0`$ and $`\mu `$ in terms of the doping fraction. For the dependence of $`\mu `$ on $`x`$, we use $`\mu x^2`$ (in units in which $`2t=1`$) as in reference, but other reasonable dependences should give similar results (cf. the discussion in the previous section). Hence, the results of this section are valid and relevant for the cuprates even if $`\mu `$ is not pinned close to zero in the underdoped regime. For $`W_0`$ and $`\mathrm{\Delta }_0`$ we choose the relations given in Eq. 20 and Eq. 21 only because they have received some recent experimental support . But it’s worth-emphasizing that the qualitative features that we extract from Eq. 43 for the behavior of the superfluid density with $`x`$ are fairly independent of the precise functional forms of $`W_0`$ and $`\mathrm{\Delta }_0`$. The only input that is needed is the existence of DDW order, with diminishing strength with $`x`$, and complementary development of the DSC order. The DDW order eats away part of the superfluid density from an otherwise pure DSC system even within the superconducting dome in the $`Tx`$ phase diagram. With increasing $`x`$, the DDW order weakens, hence the superfluid density increases in a model-independent manner.

In Fig. 3 we plot $`n_s/m^{}`$ against the doping-fraction $`x`$ using Eq. 20 and Eq. 21 for $`W_0`$ and $`\mathrm{\Delta }_0`$. To illustrate the qualitative robustness of our result, we also show, by the dashed curve, the same plot for $`W_0`$ dispersing quadratically with $`x`$ between the same two end points as in Eq. 20. As can be seen from the figure the zero-temperature superfluid density, in a model-independent manner, shows a rapid drop in the underdoped regime, similar to that observed in experiments.

The leading temperature dependence of the superfluid stiffness can be evaluated from Eq. 45. At temperatures much smaller than the relevant energy scales $`W_0`$ and $`\mathrm{\Delta }_0`$, only the nodal regions close to the points $`(\pi /2,\pi /2)`$ and symmetry-related points will contribute in the suppression of the superfluid weight. By expanding around those points, one can see that the leading temperature dependence is indeed linear for the optimally and moderately doped samples in a fairly wide range of temperatures. For these doping concentrations, where $`\mathrm{\Delta }_0`$ is larger than or comparable to $`W_0`$, $`W_0`$ plays a subleading role to $`\mathrm{\Delta }_0`$ in determining the temperature dependence of the suppression. This is due to the peculiar band structure of the problem.

On the other hand, for the heavily underdoped samples the situation is quite different. Critical fluctuations are clearly important to determine any low-temperature property of the system, for it is close to a quantum critical point. Even if we ignore these fluctuations, and rely strictly on our mean-field results, the conclusions are very different from moderate or optimally doped samples. Though in the asymptotically low temperature regime, the depletion of the superfluid density is linear in temperature, there is an intermediate temperature range over which the suppression of the superfluid density actually behaves as $`\sqrt{T}`$. As the DDW gap is much larger than the superconducting gap in these heavily underdoped samples, $`W_0`$ crosses over to produce the leading contributions in the expansion of Eq. 45 around the nodes, and is eventually responsible for the $`\sqrt{T}`$ behavior of the suppression. All of this is summarized in Fig. 4, where we have plotted the temperature-dependent part of the superfluid stiffness with temperature scaled by the superconducting gap, for six values of the doping fraction $`x`$. The plot shows $`\sqrt{T}`$ behavior for the heavily underdoped systems for the experimentally relevant temperatures. For asymptotically low temperatures for these doping values, and for a fairly wide range of temperatures for moderate or optimal doping, the data show exact collapse on a single straight line signifying a unique low-temperaure slope.

## V Discussion

The cuprate phase diagram proposed in reference is reproduced here in Fig. 1 for easy reference; we only show the relevant portion of the phase diagram, ignoring the complex set of competing charge ordered states in the underdoped regime as well as the spin glass phase.

From our calculation, the behavior of the $`T=0`$ superfluid density with $`x`$ is as expected from experiments in the underdoped regime of the high-$`T_c`$ cuprates. The striking rapid suppression of the superfluid density (and sharp increase of the London penetration depth) below $`x0.2`$ is naturally explained by the emergence of the singlet DDW at that doping-fraction. The behavior of the superfluid density near $`T_c`$ or near $`x0.05`$ (the lowest doping at which superconductivity occurs) will be strongly affected by thermal or quantum phase fluctuations, so our Hartree-Fock results will be suspect in those regimes. However, away from these critical points, we might expect our calculation to be on solid footing, and it is encouraging to see that it agrees with experimental measurements of $`n_s`$. The temperature dependence also captures the striking linearity shown in Fig. 4, and the doping independence of the slope, along with a $`\sqrt{T}`$ dependence in an intermediate temperature regime in the heavily underdoped regime akin to the experiments.

The neutron scattering peaks observed in the pseudogap and superconducting regimes are broadly consistent with our calculations. There are peaks at $`(\pi ,\pi )`$ at energies which are controlled by the single-particle gap and the doping. There are no observable peaks at $`(\pi ,0)`$ and symmetry-related points because, from our analysis, they would be at high energies controlled by the band dispersion where they are likely to be strongly damped. This resolves an earlier puzzle of the analysis of the staggered flux phase. Though we have not presented here results for finite tempeartures, from the phase diagram and our zero temperature theory, the qualitative aspects at finite temperatures can still be explained. However, the details of the peak position in energy and its doping dependence is beyond a simple Hartree-Fock calculation.

If one approaches along the path labeled (1) in Fig. 1, at first one will find no structure in the spin-spin correlator, typical of a Fermi liquid. Below $`T^{}`$, the correlator peaks at an energy $`2W_0`$, twice the maximum of the DDW gap at the wavevector $`𝐐`$. This is consistent with experiments. Below $`T_c`$, the intensity in the spin-spin correlator is amplified as the DSC order develops, as discussed in Ref. . Due to the coexistence of $`DDW`$ and $`DSC`$, the peak will be at a higher energy — $`2\sqrt{W_0^2+\mathrm{\Delta }_0^2}`$, shifted by an amount that depends on the doping. To see this, recall that, for coexisting DDW and DSC, the peak-energy shifts to higher values with $`|\mu |`$, which increases with doping.

Consider now the path (2) in Fig. 1. In the DDW state, and for small chemical potential, the neutron scattering intensity should exhibit a peak at wavevector $`𝐐`$ and energy $`2W_0`$; recall that, for DDW order, the peak-energy shows no movement with $`|\mu |`$, but is destroyed by a large $`|\mu |>W_0`$. Since the DDW gap varies with $`x`$ as in Eq. 20, the peak would shift to higher energies as $`x`$ is decreased, if $`\chi _0`$ alone were responsible for the spin-fluctuation spectrum. However, the reduced Coulomb interaction $`\overline{U}`$ is important since it shifts the peaks to lower energies, broadening them simultaneously. For smaller $`x`$, due to the lower density of the holes, the effective $`\overline{U}`$ increases, and, as a result, as shown in Fig. 2, the peak moves to lower energies. Consequently, the peak-energies in this part of the phase diagram would be influenced by these two competing effects.

Consider the path (3) in Fig. 1. In the coexisting region, DDW plus DSC, the peak energy is $`2\sqrt{W_0^2+\mathrm{\Delta }_0^2}`$, shifted by $`|\mu |`$, but now, in contrast to $`T>T_c`$, the magnitude of the total order parameter remains constant as per Eq. 21. Of course, $`\overline{U}`$ must be important as well, modifying this conclusion.

###### Acknowledgements.

We would like to thank Qiang-Hua Wang for pointing out a minor error in an earlier version of this paper. S. C. and S. T. are supported by NSF-DMR-9971138, and C. N. by NSF-DMR-9983544 and the A.P. Sloan Foundation. The work of H. -Y. K. was conducted under the auspices of the DOE, supported by funds provided by the University of California for the conduct of discretionary research by Los Alamos National Laboratory.