text:
# Modelling X-ray variability in the structured atmospheres of hot stars

## 1 Introduction

The X-ray emission from hot stars has proved to be an important “window” for investigating mass loss in early type stars, and the advent of the latest suite of orbital X-ray telescopes are providing new insights. Here we develop a model to explore the X-ray variability of early type stars, first to understand the near absence of significant variability observed so far, and second to determine the requisite S/N and time resolutions necessary to measure the expected X-ray variability with current and future instrumentation. The objects of interest are hot, luminous stars with supersonic winds that carry significant mass, such as OB and Wolf-Rayet (WR) stars. The main parameters characterising these stellar winds are the terminal velocity ($`v_{\mathrm{}}`$) and mass loss rate ($`\dot{M}`$). Typical O star mass loss rates are $`\dot{M}10^7`$$`M_{}`$yr-1 with terminal speeds $`v_{\mathrm{}}10003000`$$`\mathrm{km}\mathrm{s}^1`$. For Wolf-Rayet (WR) stars $`\dot{M}`$ is larger at around $`10^{4.5}`$$`M_{}`$yr-1 but with terminal speeds like those of O stars. Such strong outflow results in a predominance of anomalously strong and broad emission lines in the spectra of WR stars. The high mass loss rate during the WR phase has a significant effect on its evolution and the chemical enrichment and mechanical energy input to the interstellar medium.

It is widely accepted now that stellar winds of OB stars are maintained by radiation pressure in numerous lines of heavy elements (Castor et al. 1975; Pauldrach et al. 1986). But the acceleration of high mass loss WR winds is though to be due to multi-line scattering of photons (Lucy & Abbott 1993; Springmann 1994; Gayley et al. 1995). Non-LTE, moving atmosphere modelling codes have been developed (e.g., Hamann & Koesterke 2000; Hillier & Miller 1998), which adopt radiative and statistical equilibrium and assume a monotonic velocity law a priori. In the framework of the so-called “standard model” the emission lines are formed in a spherically symmetric, time-independent, dense, smooth stellar wind that is photoionized by the hot core. However, there has long been growing evidence that at least some assumptions of the standard model are oversimplified. A pointed example is the X-ray emission from hot star winds, first detected as discrete X-ray sources with Einstein (0.2–4.0 keV) (Harnden et al. 1979; Seward et al. 1979; Seward & Chlebowski 1982).

Early UV observations of hot stars revealed the appearance of superionisation, especially the presence of Ovi which would not be expected in the winds of hot stars. Cassinelli & Olson (1979) suggested that Auger ionisation by X-rays in the context of a coronal model could explain the existence of these highly ionised species. However, the hypothesis of a hot corona in the outer parts of a stellar atmosphere was ruled out due to the fact that the soft X-ray flux is not sufficiently absorbed at low energies (e.g., Cassinelli & Swank 1983). According to a phenomenological model proposed by Lucy & White (1980) and developed by Lucy (1982), dynamical instabilities should arise in line driven stellar winds leading to shock generation throughout the flow. Cassinelli & Swank (1983) found that the shock model of Lucy (1982) produced too low an X-ray flux and also that the X-rays were too soft to explain the Einstein observations of the O stars. They also pointed out that the spherically symmetric shock model would lead to strong variability of the X-rays, which was not seen. To explain the lack of clear evidence for variability of the X-rays from hot stars they proposed that the X-rays are not from spherically symmetric shocks but instead from shock fragments in the wind. Rosat (0.2–2.4 keV) observations have brought new insight in the properties of X-ray emission from hot stars. Kudritzki et al. (1996) have found interesting relations between the X-ray emission and the wind properties from an analysis of 42 O star spectra. For example, they confirm the relation between X-ray luminosity and stellar bolometric luminosity as $`L_\mathrm{X}10^7L_{\mathrm{Bol}}`$ that was previously known, and also find that $`L_\mathrm{X}`$ scales with the wind kinetic luminosity $`0.5\dot{M}v_{\mathrm{}}^2`$. For the WR stars, the data is of lower quality and consists almost entirely of passband fluxes. However, based on Einstein observations, Pollock (1987) was able to find that on average, N-rich WR stars (WN) tend to be more luminous than the C-rich WR stars (WC), possibly attributable to the difference in abundances of the two types. The Rosat data have also revealed that, unlike their progenitors the O stars, the X-ray luminosities of single WR stars are not correlated with $`L_{\mathrm{Bol}}`$, wind momentum $`\dot{M}`$$`v_{\mathrm{}}`$, wind kinetic luminosity $`0.5\dot{M}v_{\mathrm{}}^2`$, or WR subtype (Wessolowski 1996; Ignace & Oskinova 1999). Recent data from Chandra and XMM-Newton of OB stars (Schulz et al. 2000; Kahn et al. 2001; Waldron & Cassinelli 2001) and WR stars (Maeda & Tsuboi 1999) are just becoming available. These newer data, mostly of O stars, are showing rich X-ray line spectra, the analysis of which are yielding some unexpected results. We refer the reader to the above mentioned papers for a description of these data and here concentrate on the topic of X-ray variability in single stars, a topic that has not much been addressed from a theoretical perspective.

Dynamical instabilities in line-driven winds have been studied extensively for OB stars (Owocki & Rybicki 1984; Owocki et al. 1988; Rybicki et al. 1990). Model computations predict shock velocity-jumps ranging from 500 km s-1 to 700 km s-1, implying post-shock temperatures and emission measures that are marginally adequate to explain the observed thermal X-ray emission (e.g. Hillier et al. 1993, Feldmeier et al. 1997a). Gayley & Owocki (1995) have shown that similar to OB stars, the instability mechanism should also operate in the denser WR winds.

In terms of observed X-ray variability, Crowther & Willis (1996) have presented a review of available data based on Rosat observations. For $`\zeta `$ Ori, Berghöfer & Schmitt (1994) claimed to detect X-ray variability over a 2-day period, reaching about 30% in hard X-rays (0.64-2.38 keV), while the soft band data (0.16-0.5includegraphics1 keV) remained about constant. For $`\zeta `$ Pup, Berghöfer et al. (1996) found evidence for a modulation in the harder (0.9–2.0 keV) X-ray emission and in the optical data with a 16.7 hour period, and suggested that this reflects periodic changes in the base wind density of the O4 f star. On the other hand, the study of $`\sigma `$ Ori showed no significant X-ray variability. Except for WR 6 (HD 50896) and WR 1 (HD 4004), no X-ray variability studies for single WR stars have been undertaken. For WR 6, Willis & Stevens (1996) reported variability at the $``$ 30% level on timescale of $`1`$ day, together with larger epoch changes. No significant spectral shape changes were found. In the case of WR 1, no evidence of significant X-ray variability was reported by Wessolowski & Niedzielski (1996). Thus far, no short-time scale X-ray variability on the order of the wind flow time ($`R_{}/v_{\mathrm{}}`$ hours) has been observed in O and WR stars.

In this paper we consider the application of stochastic inhomogeneous wind models, as widely accepted necessary components to explain spectral variations at optical and UV wavelengths, to the X-ray emission from single O and WR type stars. We restrict our analysis to shock models of X-ray formation in hot stars and do not consider alternative models (see discussion in Waldron & Cassinelli 2001). In Sect. 2, we describe a stellar wind model consisting of numerous zones filled by hot X-ray emitting gas. Two following sections use two different approximations of the optical depths of the cool material, which absorbs X-ray emission. In Sect. 3 we apply the so-called exospheric approximation treating the X-ray emitting zones as spherical shells. In Sect. 4 we generalise our model by using formal integration of the radiative transfer equation for randomly distributed wind shocks. We investigate the crucial parameters of the variability and their relation to the input parameters of the model. A summary discussion of the results is presented in Sect. 5.

## 2 Stochastic model for the distribution of X-ray emitting material.

We consider a spherically symmetric and time -independent stellar wind that is a mix of “cool” and “hot” gas in dynamical equilibrium. The minor hot gas component gives rise to X-ray emission. We assume that hot gas is present in the form of spatially separated compact volumes characterised by temperature $`T_\mathrm{X}`$ and density which are different from the temperature and density of the ambient cool wind gas.

It is common in the literature to attribute wind “clumps” to dense formations with temperatures much lower then the temperatures of the X-ray emitting material, sometimes referred to as DWEEs (discrete wind emission elements, Lépine & Moffat 1999). These clumps are presumably responsible for observed line profile variability in the optical spectra of WR stars, or the formation of DACs (discrete absorption components) observed in the optical and UV spectra of some O stars (Brown et al. 1995). Such cold clumps and hot X-ray emitting zones may well have related origins. Feldmeier et al. (1997a) considered the formation of X-ray emitting zones as a consequence of shell collisions due to the propagation of reverse shocks in the stellar wind. In their hydrodynamical simulations they demonstrated the existence of dense shells moving according to the stationary wind velocity law and small-mass, high-velocity blobs. Copious X-rays are produced when a fast blob rams into an outer shell. These X-rays thus originate over a fairly narrow range in radius. In this respect, recent observations (Schulz et al. 2000) seem to confirm such a scenario, except that the expected strong variability is not observed.

The general description of a source with a continuous spatial distribution of temperature (e.g., a radiative cooling zone behind the front of a shock wave) is given by the differential emission measure which includes integration over all surfaces with the same temperature, such as over multiple independent shocks in a stellar wind (Craig & Brown 1976). But it is necessary to note that an overall consistent picture of dynamical, structured, thick stellar winds is still far from being complete. Although detailed models do not predict the X-ray sources to be isothermal with unique temperature $`T_\mathrm{X}`$ amongst them, the isothermal approximation for the overall wind structure should be fairly close to the actual case because of the short radiative cooling times in these dense stellar winds (Feldmeier et al. 1997b). The two-temperature approximation with isothermal X-ray emitting zones embedded in a cooler outflowing wind is able to reasonably well reproduce the observed X-ray spectra from hot stars (Hillier et al. 1993) and might be used to place constraints on physical parameters such as the volume filling factor (e.g., Oskinova et al. 2001).

Therefore, we will consider a spherically symmetric and, on average, time independent outflow that is a homogeneous mix of “cool” and “hot” gas in dynamical equilibrium. For isothermal optically thin hot X-ray zones, the particular choice of $`T_\mathrm{X}=10^7K`$ tends to maximise the X-ray emission in the ROSAT band (see Ignace et al. 2000). We use the Raymond-Smith cooling function ($`\mathrm{\Lambda }_{\mathrm{RS}}`$) for optically thin plasmas (Raymond & Smith 1977). For WR stars, a rough correction of $`\mathrm{\Lambda }_{\mathrm{RS}}`$ for non-solar abundances is given in Ignace et al. (2000). The ambient cool medium ($`10^5`$ K) is as described by the standard model and is optically thick to X-rays.

In the two-temperature approximation, the ratio of emission measures of hot and cool components of the wind, $`\mathrm{EM}_\mathrm{X}/\mathrm{EM}_\mathrm{w}`$, is proportional to the ratio of volumes filled by hot ($`V_\mathrm{X}`$) and cold ($`V_\mathrm{w}`$) material and to the square of their density ratio $`D_\mathrm{s}`$. Thus, the “filling factor” is

$$f_\mathrm{X}=\frac{\mathrm{EM}_\mathrm{X}}{\mathrm{EM}_\mathrm{w}}D_\mathrm{s}^2\left(\frac{\mu _\mathrm{e}^\mathrm{w}}{\mu _\mathrm{e}^\mathrm{X}}\right)^2\frac{V_\mathrm{X}}{V_\mathrm{w}},$$
(1)
where $`\mu _\mathrm{e}^\mathrm{w}`$ and $`\mu _\mathrm{e}^\mathrm{X}`$ are the mean mass per free electron of the cool or standard wind component. The X-ray emission arises primarily from two-body processes–collisional excitation and radiative recombination. Therefore, the volume emissivity $`j_\nu `$ (energy per unit time per unit mass) is proportional to the square of the density,

$$4\pi j_\nu =f_\mathrm{X}\rho _\mathrm{w}^2\frac{\mathrm{\Lambda }_\nu (T_\mathrm{X})}{\mu _\mathrm{e}^\mathrm{X}\mu _\mathrm{i}^\mathrm{X}m_\mathrm{H}^2}$$
(2)
where $`\mu _\mathrm{i}^\mathrm{X}`$ is the mean molecular weight per ion of the hot gas and $`m_\mathrm{H}`$ is the proton mass. $`\mathrm{\Lambda }_\nu `$ (erg cm3s-1) is the cooling function at the energies $`E`$ concerned and $`\rho _\mathrm{w}`$ is the density of the standard wind is

$$\rho _\mathrm{w}(r)=\frac{\dot{M}}{4\pi r^2v(r)}.$$
(3)
By the definition of filling factor in Eq. (1), the total specific luminosity (erg s-1) emerging from the wind

$$L_\mathrm{X}(E)=f_\mathrm{X}\mathrm{\Lambda }_\nu (T_\mathrm{X})\mathrm{EM}_\mathrm{w}.$$
(4)
This then is the monochromatic X-ray luminosity at any instant in time, since our hypothesis is essentially that the emission measure of X-ray emitting material (or correspondingly the filling factor) is time variable.

It was pointed out by Cassinelli et al. (1996) that the X-rays from a wind shock will be seen only after that shock has moved into a region where the cool wind attenuation is optically thin to the X-rays. As we shall show, the wind is quite opaque at most X-ray energies up until the flow has reached a substantial fraction ($`>50\%`$) of its terminal speed. Thus, we can assume that the X-rays emerge primarily from distances where $`v(r)v_{\mathrm{}}`$. For the wind attenuation, K-shell absorption by metals in the cool wind is the dominant opacity source. Assuming photo-electric absorption by K-shell electrons

$$\kappa _\mathrm{w}(E)=\sigma _\nu (E)/\mu _\mathrm{i}^\mathrm{w}m_\mathrm{H}=\frac{1}{\mu _\mathrm{i}m_\mathrm{H}}\underset{\mathrm{j}}{}\frac{n_\mathrm{j}}{n_\mathrm{i}}\sigma _\mathrm{j}(E),$$
(5)
where $`\mu _\mathrm{i}^\mathrm{w}`$ is the mean molecular weight per ion of the cool wind and $`\sigma _\mathrm{j}(E)`$ is the photo-absorption cross-section. (Actually, the appropriate molecular weight to use is that per nucleus, but in hot star winds, there is no neutral gas, and so every nucleus is an ion.) The ratio $`n_\mathrm{j}/n_\mathrm{i}`$ defines the relative abundance by number for atomic species j, where $`n_\mathrm{i}`$ is the number density of ions. The X-ray absorption (apart from prominent K-shell edges) can be represented roughly by a power-law in energy: $`\kappa _\mathrm{w}(E)\kappa _0E^\gamma `$, with $`\gamma `$ in the range 2–3, depending on ionisation structure and chemical composition (see, e.g. Hillier et al. 1993; Cohen et al. 1996). The parameter $`\kappa _0`$ is a constant that depends sensitively on the abundances.

To proceed, we now develop a phenomenological model in which the wind inhomogeneities develop at random times near some inner boundary $`R_0`$. Then they propagate radially according to a monotonically increasing velocity law $`v(r)`$. Let us introduce dimensionless notations for velocity $`u=v(r)/v_{\mathrm{}}`$, and for distance $`x=r/R_0`$. Further, we normalise time to the characteristic time $`𝒯_\mathrm{f}=R_0/v_{\mathrm{}}`$, which is equal to the flow time in the specific case $`R_0=R_{}`$. Then for unitless time we have $`\mathrm{t}=𝒯/𝒯_\mathrm{f}`$.

To avoid confusion we would like to clarify here, that further out, we shall use two descriptions for the distribution of the hot gas. In the exospheric approximation described in Sect. 3, we assume that the hot material is present in the form of spherical shells. The exospheric approach with spherical shells are simplifying assumptions, but the model does allow us to gain some insight into the nature of X-ray variability from wind shocks. In Sect. 4, we allow for randomly distributed wind shocks of arbitrary geometrical form (as long as the hot zones are not exceedingly “large”, to be defined later), and we consider the full radiative transfer problem for the emergence of X-rays from the wind flow.

For the motion of the hot zones (e.g., shells), we consider an ensemble of zones filled by hot gas propagating in the radial direction with constant velocity $`u=1`$. The different zones in our model are labelled by their times of appearance $`\mathrm{t}`$ at the distance $`x_0=1`$. The first zone is at distance $`x_0=1`$ at the moment of time $`\mathrm{t}_0=0`$. At time $`\mathrm{t}_k`$ when shell number $`k`$ is crossing $`x_0=1`$, the first zone is at distance $`x^1(\mathrm{t}_k)=x_0+\mathrm{t}_k`$ with $`\mathrm{t}_k=_{i=1}^k\delta \mathrm{t}_i`$ and zone number $`j`$ is at distance $`x^j(\mathrm{t}_k)=x_0+(\mathrm{t}_k\mathrm{t}_j)`$. Here, we choose to use the time intervals $`\delta \mathrm{t}_i`$ as drawn from uniformly random numbers in the range $`[0.5,\mathrm{\hspace{0.17em}1.5}]`$ multiplied by a dimensionless parameter $`T`$:

$$\delta \mathrm{t}_i=𝒫([0.5,1.5])T,$$
(6)
and so $`\mathrm{t}_i=T`$. The parameter $`T`$ is in fact an average separation in time between the appearance of subsequent zones and is measured in flow times $`𝒯_\mathrm{f}`$. Further we will explore the influence of the average time separation $`T`$ on the character of the variability.

The parameter $`T`$ is introduced because detailed models of X-ray production demand the presence of an initial perturbation of the stellar wind. For example, in models which depict X-ray production mechanisms in O stars (Feldmeier et al. 1997b), the existence of large scale turbulence near the wind base is postulated. As it will become clear below, the average time separation $`T`$ is the only essential parameter of the models under consideration. We thus wish to emphasis that in principle, one can expect $`T`$ to be an observable value, which could be inferred from the X-ray emission light curve analysis as the average separation between the peaks of intensity.

The emission measure of each zone is $`\delta \mathrm{EM}=n_\mathrm{e}^\mathrm{X}n_\mathrm{i}^\mathrm{X}\delta V`$. Using Eq. (3), the squared number density inversely scales with the fourth power of distance: $`n_\mathrm{e}n_\mathrm{i}=n_0^2x^4`$ , where

$$n_0^2=D_\mathrm{s}^2\left(\frac{\dot{M}}{4\pi v_{\mathrm{}}R_0^2}\right)^2\frac{1}{\mu _\mathrm{e}^\mathrm{x}\mu _\mathrm{i}^\mathrm{x}m_\mathrm{H}^2}.$$
(7)
On the other hand, the volume of a zone $`\delta V`$ with fixed radial thickness and solid angle grows with the square of distance $`\delta Vx^2`$. This implies that the contribution of a zone to the total X-ray luminosity determined by its emission measure $`\delta \mathrm{EM}`$ has an overall decline that goes as $`x^2`$. This is a plausible treatment when mass conservation for the hot shocked material is assumed so that the continuity equation Eq. (3) can be applied.

## 3 Exospheric model for the time dependent X-ray emission of single stars

### 3.1 Modelling time dependent X-ray emission

Let us assume that the observed X-ray emission arising from the hot gas emerges only from radii exterior to the optical depth unity surface of radius $`r_1`$ (exosphere), with X-rays at smaller radii completely attenuated. In this approximation the X-ray luminosity has the proper scaling with mass loss and other key wind parameters. The exospheric radius $`r_1`$ of optical depth unity is (Owocki & Cohen 1999)

$$r_1(E)=\frac{\dot{M}}{4\pi v_{\mathrm{}}}\kappa _\mathrm{w}(E).$$
(8)
The emission measure of the cool material outside the exosphere which is assumed to be transparent to the X-rays is given by

$$\mathrm{EM}_\mathrm{w}4\pi _{r_1}^{\mathrm{}}n_\mathrm{i}^\mathrm{w}n_\mathrm{e}(1\mathrm{W}(r_1))r^2𝑑r,$$
(9)
where $`\mathrm{W}(r_1)`$ is a dilution factor defined to be $`\omega _{r_1}/4\pi `$ where $`\omega _{r_1}`$ is the solid angle subtended by the exosphere and the parenthetical term accounts for geometric occultation by the exospheric surface of radius $`r_1`$.

Let us make the simplifying assumption that X-ray emitting material is present in the form of spherical shells. All shells have the same small fractional thickness $`=\mathrm{\Delta }R_{\mathrm{shell}}/R_0`$ at a given radius in the envelope. Recall that we assume isothermal X-ray emitting zones and thus the parameter $``$ is not to be misconstrued as a cooling length but rather the geometrical fractional thicknessof a shell filled with the hot gas.

In our phenomenological exospheric model, X-ray emitting shells are seen only beyond radius $`r_1=x_1R_0`$, so that effectively $`x_0x_1`$. The shells start their motion outwards in the atmosphere at quasi-random times $`\mathrm{t}`$k. At that moment, the X-ray flux ($`\delta F_\mathrm{X}`$) is increased by the new source. Then owing to the decrease in emission measure of the zone with distance ($`\delta \mathrm{EM}_\mathrm{X}x^2`$ for constant shell thickness), the incremental flux $`\delta F_\mathrm{X}`$ decreases. Thus, the total X-ray output of an ensemble of X-ray emitting zones will be observed to vary with time.

The radius $`x_1`$ has a strong dependence on mass-loss rate, chemical composition and energy (see Fig. 1). This means that in the model under consideration, the X-ray emitting zones will cross the surface of optical depth unity at different distances for different energies. Therefore their emission measures at the moment when they appear for an external observer will be different. Obviously then, one may expect the character of variability to change as a function of wavelength and for stars of different abundances.

The emission measure of an individual shell $`j`$ at time $`\mathrm{t}_k`$ is

$$\delta \mathrm{EM}_j\mathrm{\hspace{0.17em}4}\pi R_0^3n_\mathrm{e}^\mathrm{X}n_\mathrm{i}^\mathrm{X}\left[x_1(E)+(\mathrm{t}_k\mathrm{t}_j)\right]^2.$$
(10)
Combining Eqs. (7) and (10), the emission measure of the ensemble of hot shells becomes includegraphicsincludegraphics

$$\mathrm{EM}_\mathrm{X}(E,\mathrm{t}_{k_{\mathrm{max}}})\mathrm{\hspace{0.17em}4}\pi R_0^3n_0^2\underset{k=1}{\overset{k_{\mathrm{max}}}{}}\frac{1}{[x_1(E)+\mathrm{t}_k]^2}.$$
(11)
Here the integral over the volume is replaced by summation over all X-ray emitting shells outside the sphere of $`x_1(E)`$, and $`k_{\mathrm{max}}`$ is the maximum number of shells under consideration in the wind. The value of $`k_{\mathrm{max}}`$ is chosen so that the most distant shells make a negligibly small contribution to the emission measure (i.e., $`\delta EM_\mathrm{X}(x_{\mathrm{max}})0`$ for $`xx_1`$). We choose $`k_{\mathrm{max}}10^5`$ which allows us to account for the emission of a shell till it reaches a distance $`x_{\mathrm{max}}=k_{\mathrm{max}}T10^5`$.

The X-ray luminosity with time is $`L_\mathrm{X}(E,t)=\mathrm{EM}_\mathrm{X}(E,\mathrm{t})\mathrm{\Lambda }_\nu (T_\mathrm{X})`$. Given that hot gas zones are a function of radius only, and not of latitude or azimuth, a relatively simple expression for the relative variability, $`\sigma _\mathrm{L}/L_\mathrm{X}`$, can be derived. Recognising that the X-ray emission from any given zone drops rather steeply with radius, both the emission and variability will be predominantly governed by zones that appear at the $`x_1`$ radius. Hence the fractional variability can be estimated based on Poisson statistics. The number density of zones above a given radius $`x`$ scales as $`1/x_1^2`$; the volume in that vicinity scales as $`x_1^3`$; so the number of zones near $`x_1`$ is $`Nx_1`$. Consequently the relative variability scales as $`\sigma _\mathrm{L}/L_\mathrm{X}1/\sqrt{x_1}`$. Given that $`x_1\kappa _\nu (E)`$, we can expect the relative variability to be a function of energy. For purposes of illustration, if $`\kappa _\nu (E)E^\gamma `$, then the relative variability will be $`\sigma _\mathrm{L}/L_\mathrm{X}E^{\gamma /2}`$, which is an increasing function of energy.

Figure 2 shows changes in the emission measure of the shell ensemble, computed using Eq. (11). We took the values of $`x_1`$ corresponding to 1.5 keV for O, WN and WC stars. As is clear from the figure, the character of variability depends on the chemical composition. As can be seen from Fig. 1, the radius of optical depth equal unity near 1 keV for a WC star is located about an order of magnitude further away in the wind than is the case for a WN star. This leads to an overall reduction of the X-ray emission, and some suppression of the variability in WC stars at this energy. In Sect. 4 we revisit these simulations using the formal solution of radiative transfer with angle dependent optical depth, instead of the simplified exospheric approximation. But first we make some comments on model constraints obtainable from filling factor data.

### 3.2 Using filling factors to constrain model parameters

In Ignace et al. (2000), lower limits to the filling factors of most putatively single Galactic WR stars were determined from an analysis of Rosat All-Sky Survey (RASS) observations. A natural question arises: Is it possible to extract information about the spatial distribution of hot, X-ray emitting gas by knowing $`f_\mathrm{X}`$ from RASS passband observations in the framework of the model under consideration? Let us again assume that the fractional thickness of the X-rayemitting layers $``$ does not depend on distance. By definition of Eq. (1) and using Eq. (9) and Eq. (11), we have:

$$f_\mathrm{X}=D_\mathrm{s}^2\left(\frac{\mu _\mathrm{e}^\mathrm{w}}{\mu _\mathrm{e}^\mathrm{X}}\right)^2\frac{_{k=1}^N_E(x_1+\mathrm{t}_k)^2𝑑E}{_Ex_1^2𝑑E}.$$
(12)
Here, we neglect the occultation term in Eq. (9). Note that $`x_1=x_1(E)`$ is implicit. One infers from Eq. (12) that in fact there are two dependent parameters in the problem: $`T`$ and $``$. Let’s assume that it is possible to infer $`T`$ from observations of the variability of X-ray emission. Then making assumptions about the density ratio $`D_\mathrm{s}`$, one may determine $``$ if $`L_\mathrm{X}`$ (equivalently $`f_\mathrm{X}`$) is known. Therefore from Eq. (12), we may estimate $``$ from known values of $`f_\mathrm{X}`$. In Table 1 we give typical values of $``$ using $`D_\mathrm{s}=4`$ (the same value as in Baum et al. 1992 and Hillier et al. 1993) for WN and WC stars.

So far, we have used the exospheric approximation to demonstrate that the X-ray variability is strongly influenced by the opaque cool material of the stellar wind and therefore (a) is a function of energy and (b) depends on the chemical composition of the stellar wind (see Fig. 2). Also, it was shown that lower limits to filling factors can be used to infer lower limits to the thickness of spherical X-ray emitting shells (as illustrated in Tab. 1), if the average time interval between subsequent shocks $`T`$ is known from observations.

## 4 Formal solution for X-ray emitting zones

While it appears that the majority of WR and O stars display basically spherical winds on the large scale, there is compelling observational evidence that the winds are clumped on small scales (e.g. Lépine & Moffat 1999). In Lépine & Moffat (1999) a phenomenological model was introduced, wherein WR winds are depicted as consisting of a large number of randomly distributed, radially propagating, discrete wind emission elements. This model was used to simulate line profile variability patterns in emission lines from a clumped wind. We adopt this approximation for the X-ray producing region and consider the X-ray emission as arising from analogous optically thin zones with temperature $`T_\mathrm{X}`$, which as before are embedded in the cooler bulk wind flow. However, the zones are no longer treated as spherical shells.

We treat these hot zones as piecewise elements, presumably produced by shocks. They expand adiabatically moving radially according to a $`\beta `$-law for the velocity. Zones start their motion at randomly distributed moments of time from random but uniformly distributed positions at an imaginary spherical surface of radius $`R_0`$ – the minimal predicted radius for shocks to occur.

The main difference with the exospheric model from the previous section is that the optical depth of cool material now depends not only on radius but also on the spherical polar angle as seen by an observer. Further we now compute a proper line-of-sight optical depth to each hot emitting zone, in contrast to the exospheric approximation for which hot zones are either completely attenuated, or completely unattenuated. To simplify the numerical simulation, we do require that the emission measure of every zone $`\delta \mathrm{EM}`$ is the same at a given radius and that the spatial size of each zone is small enough to assume the wind attenuation to every point in the zone is approximately constant (e.g., the zones cannot be hemispherical shells, since the wind attenuation would strongly vary across such a structure).

In this section we adopt the usual cylindrical coordinate system (e.g. Lamers et al. 1987). The coordinates are: distance from the star $`x`$, impact parameter for the line-of-sight $`p`$, and distance along the line of sight $`z`$, all normalised to the stellar radius. The angle between the line of sight and the radial vector is $`\theta =\mathrm{arccos}\mu `$, with $`\mu =1`$ if the radial vector points to the observer. For any point in the wind, $`x`$, $`p`$ and $`z`$ are related by

$$z=\mu x=\mu (p^2+z^2)^{1/2}$$
(13)
With our assumptions, the emission measure changes with distance from the star as $`\delta \mathrm{EM}_\mathrm{X}=\delta \mathrm{EM}_0(1/x)^2`$, with $`\delta \mathrm{EM}_0`$ the emission measure of a zone at distance $`x_0`$. The luminosity of each zone for which the optical depth of the ambient wind is $`\tau _\mathrm{w}`$ is

$$\delta L_\nu =\mathrm{\Lambda }_\nu (T_\mathrm{X})\delta \mathrm{EM}_\mathrm{X}e^{\tau _\mathrm{w}}$$
(14)
and the luminosity of a hot gas zone occulted by the stellar core is assumed to equal zero. Applying the equation of continuity Eq. (3) with Eq. (5), we can derive the wind optical depth to be

$$\tau _\mathrm{w}=\frac{\sigma (E)\dot{M}}{4\pi R_0v_{\mathrm{}}\mu _\mathrm{i}m_\mathrm{H}}_z^{\mathrm{}}\frac{dz}{x^2u(x)}\tau _0(E)_z^{\mathrm{}}\frac{dz}{x^2u(x)},$$
(15)
where $`x=\sqrt{p^2+z^2}`$. Now the optical depth for each zone is angle dependent. For soft (0.4–2.4 keV) X-ray energies at radii where the wind acceleration takes place, the optical depth is generally quite large (see Tab. 2).

Although our numerical calculations can account for the X-ray emission from zones at small radii, Table 2 indicates that their contribution to the total X-ray flux will be quite negligible owing to the rather extreme attenuation. So, we simplify our calculations by considering only those zones which have already reached the constant velocity region in their motion throughout the wind, hence $`u=1`$. Using Eq. (15), it is possible to obtain an analytic expression for the optical depth (MacFarlane et al. 1991):

$$\tau _\mathrm{w}=\tau _0(E)_z^{\mathrm{}}\frac{dz}{x^2}=\frac{\tau _0(E)}{x}\frac{\theta }{\mathrm{sin}\theta }.$$
(16)
So, the optical depth depends on the location of a given zone in the wind as well as on the energy.

In Fig. 3 the change of X-ray flux for a single X-ray emitting zone along with its motion through the flow is shown. The luminosity grows rapidly with distance owing to the decrease in optical depth by the overlying wind. In general the zones which are propagating in directions close to $`\mathrm{cos}\theta =1`$ have higher luminosity due to the fact that the optical depth is smaller in this direction. The emission rises to a peak followed by a decline. At this point, the attenuation is negligible, and the decreasing emission is a consequence of the inverse square fall-off of the emission measure with distance. The distance $`x_{\mathrm{peak}}`$ at which an X-ray emitting zone yields peak emission depends strongly on the angle $`\theta `$. Further, noting that $`\delta \mathrm{EM}_\mathrm{X}x^2`$ and $`\tau _\mathrm{w}x^1`$ it can be shown (from an analysis of Eq. (14)) that the maximum of the curve occurs when $`\tau _\mathrm{w}(E,\theta )=2`$.

Let us assume that a stellar wind contains numerous X-ray emitting zones. Clearly an ensemble of such zones will also produce time variable X-ray emission. The rise in luminosity of a single zone is much sharper than its decrease after peaking. Suppose that at a given time and a given energy, the X-ray luminosity of only one zone could reach its maximum. Subsequently, the emission from this zone will be decreasing, and until the next zone reaches its maximum, the emission of the whole ensemble will drop. Thus, in the regime of constant expansion the separation between maxima of emission reflects the average time separation $`T`$. We may expect larger levels of variability on a short time scale (order of flow times) for larger values of $``$. So, we may conclude that the observed lack of variability suggests values of $`T1`$ and therefore values of $`10^2`$, which are marginally consistent with estimations of cooling lengths (e.g. Hillier et al. 1993). It is necessary to point out here that the available observations so far have not been capable of detecting such small fluctuations of X-ray flux.

To proceed further, let us recall that the total X-ray luminosity depends on the random variable $`\mathrm{t}_i`$. For the angular distribution of hot zones, we select $`\mathrm{cos}\theta `$ as a uniform random variable in interval $`[1,1]`$, and for the azimuth, $`\varphi _k`$ is uniform random in the interval $`[0,2\pi ]`$. It should be noted that although $`\mathrm{cos}\theta _k`$ is uniform, $`\tau _{\mathrm{w},k}`$ is not. The total emission for an ensemble of zones is a summation over all contributors with the appropriate attenuation:

$$L_\mathrm{X}^{\mathrm{tot}}(E)=\mathrm{\Lambda }_\nu (T_\mathrm{X})\delta \mathrm{EM}_0\underset{k}{\overset{k_{\mathrm{max}}}{}}\frac{1}{(1+\mathrm{t}_\mathrm{k})^2}e^{\tau _{\mathrm{w},k}(E)}.$$
(17)
The wind will be more opaque to X-rays at soft energies than to those at hard energies. Therefore, we expect the phenomenological picture wherein an X-ray emitting zone has peak flux at higher energy earlier than it has peak flux at softer energy.

Using $`\delta \mathrm{EM}_0`$ as a parameter of the model, we avoid direct references to the density of hot material or volume of the zones filled by this material. In our model at each moment of time, only one X-ray emitting zone is at a given distance from the inner boundary $`R_0`$. This is similar to the propagation of subsequent spherical shells which was described in the previous section in the framework of the exospheric approximation. In Section 3.2 we have shown that the fractional thickness of a spherical layer $``$ and $`T`$ are not independent parameters but coupled with each other (see Eq. 12). To place constraints on $`\delta \mathrm{EM}_0`$, let us assume now that $`\delta \mathrm{EM}_0`$ is the same as the emission measure of a spherical shell with thickness $``$ located at distance $`R_0`$ with corresponding density $`n_0`$ as in Eq. (7).

The number of hot gas zones and their initial emission measure $`\delta \mathrm{EM}_0`$ is set by the average separation in time $`T`$. The assumption of the constant filling factor leads to the correlation between average time interval and initial emission measures of hot zones. This means that small hot zones with small separation in time should produce small changes in X-ray flux. On the other hand, large sized hot zones produce changes in X-ray flux with a more significant amplitude. To complete the picture, it is necessary to bear in mind that the degree of variability also depends on chemical composition.

Fig. 4 presents numerical simulations of short-scale time variability of the monochromatic X-ray luminosity at 1 keV for two stars with different cool wind opacities, namely O stars and WN stars. Clearly, the character of variability strongly depends on the average time separation $`T`$ and the spectral type of the star. The mass loss rates of WN stars are much higher then those of O stars; therefore, we attribute the difference in attenuation between WR stars and O stars as the main cause for different character of variability shown in Fig. 4 owing to the fact that $`x_1`$ is much bigger for WN stars than for O stars.

The changes in luminosity seen in Fig. 4 are, in fact, due to the superposition of light curves for many zones propagating in different directions similar to those shown in Fig. 3. The average separation between subpeaks on the curves from Fig. 4 reflects the average time between two subsequent zones having peaks in their luminosity. However, when $`T`$ is small, the flux fluctuations are negligible (order $``$ few per cent) and simply cannot be resolved. The interesting question to address is what is the plausible range for the parameter $`T`$. Observed fluctuations in X-ray emission for $`\zeta `$ Ori (O4f) (increase in the count rate of $``$ 30% for 2 days) and $`\zeta `$ Pup (O9.5Ia) (modulations with a period of 16.7 hours and amplitude $``$ 10%, Berghöfer & Schmitt 1994, Berghöfer et al. 1996) suggest rather small values of $`T`$ ($`1`$ flow time), assuming that the X-rays form in wind shocks.

As shown by our numerical simulations, on the time scale of several hours, the variability of X-ray luminosity might be negligible for stars of spectral type WN. The reasons for this apparent lack of variability are that in addition to small size and short time separation of hot zones, the total optical depth for the indicated energies is quite large in the WN star winds. As can be seen from Eq. (17), the exponential term $`e^{\tau _{\mathrm{w},\mathrm{i}}}`$ suppresses the diference in values of $`L_\mathrm{X}`$ . Although not shown in Fig. 4, our modelling reveals that the X-ray luminosity of WC stars is almost constant due to their quite opaque stellar winds. As obvious from the figure, even the transparent winds of O stars may not demonstrate detectable levels of variability in X-rays on short time scales of about 1 hour in the case of small $`T`$.

Although on a time-scale of hours, the X-ray variability may be quite small, the dynamical flow time $`r/v(r)`$ for distances $`r100R_{}`$ is some $`10^5`$ seconds, of order a day. So we have performed simulations of changes in flux over longer time intervals. The results of this calculation are shown in Fig. 5. Surprisingly, even for $`T=1`$, the variability becomes significant ($``$ 10%) on time scales of thousands of hours. In this case, the separation between peaks in the light curves do not reflect the time separation between X-ray emitting zones. To understand the character of this variability let us consider Eq. (17). The time dependent luminosity is governed by two uniformly distributed random variables $`\delta \mathrm{t}_i`$, and $`\mathrm{cos}\theta _i`$. Let us assume, for simplicity, that zones are launched within constant time intervals equal to unity and concentrate on random variations in $`\tau _\mathrm{w}`$. Then $`\delta \mathrm{t}_i=1`$ and from Eq. (17)

$$L_\mathrm{X}^{\mathrm{tot}}\underset{k}{\overset{\mathrm{}}{}}\frac{1}{(1+k)^2}e^{\tau _{w,k}}.$$
(18)
From Eq. (16), $`\tau _{w,k}\tau _0\theta _k/k\mathrm{sin}\theta _k`$ and $`\mathrm{cos}\theta _k`$ is a random variable from $`[1,1]`$. Thus, when the number of realisations $`k`$ is large enough, even for rather opaque winds with $`\tau _01`$, the variance of the flux is significant.

The likelihood of detecting variability of X-ray flux increases drastically at hard energies. This is due to the strong dependence of the optical depth on energy. To illustrate this, if $`\kappa _\mathrm{w}(E)E^\gamma `$ then using Eq. (18) and neglecting the second order terms, the standard deviation of the X-ray luminosity becomes $`\sigma _\mathrm{L}E^{\frac{\gamma }{2}}(\dot{M}/v_{\mathrm{}})^{\frac{1}{2}}`$. This is an increasing function of energy and decreasing function of the stellar wind density. As numerical simulations show, in the case of random times $`\mathrm{t}`$i, the general trend of the dependence of variance on energy and density is an increasing function of energy. The only condition is that the variance of emission measures of the hot zones should be smaller than the variance of the optical depths associated with them.

Fig. 6 represents the ratio of the standard deviation $`\sigma _\mathrm{L}`$ to the average value of the X-ray luminosity as a function of energy, computed using Eq. (17). Apart from ionisation edges, Fig. 6 clearly shows that the relative variability is an increasing function of energy. Fig. 6 confirms the basic scaling of variability with parameter $`x_1`$ that was obtained under the exospheric approximation with $`\sigma _\mathrm{L}/L_\mathrm{X}1/\sqrt{x_1}`$. Presumably, this will fail for large values of $`Tx_1`$, because then there are relatively few hot zones in the outflow. At very low $`T`$ one hits another regime where the flow approaches homogeneity.

As seen in Fig. 6 for $`T=1`$, the amplitude of variability for O stars can reach as much as 80% near 1.5 keV but only 10% at softer energies around 0.5 keV. At the same time, the relative variability is only a few percent near 1 keV for the WN stars. For WC stars, the level of X-ray variability is negligible. Therefore, it is not surprising that with such low sensitivity X-ray emission detectors such as Rosat PSPC (0.2–2.4 keV) or Einstein IPC (0.2–4.0 keV), the detected X-ray flux of WR stars appeared to be fairly constant X-ray sources.

## 5 Summary and conclusions

We present simulations of the expected X-ray variability for early-type stars in the framework of a shock model for X-ray production. We assumed that the optically thin hot X-ray emitting material is embedded in a cool X-ray absorbing stellar wind, described by the standard model. For such a medium we used the concept of filling factor, which is the ratio of emission measures of X-ray emitting and absorbing material. Assuming an isothermal X-ray emitting hot gas, the filling factor is proportional to the ratio of volumes filled by hot and cool components and to the square of the ratio of densities of hot and cool material. To examine the basic processes of X-ray emission and absorption we considered two different models.

To begin, we employed the exospheric approximation with angle independent optical depth. In this approximation we used a model of the stellar wind where hot gas is present in the form of spherical shells propagating with constant velocity. The emission of such an outflow is time dependent. This allowed us to derive a expression restricting the fractional thickness of a spherical shell emitting X-rays by using filling factors. The lower limits on filling factors for isothermal X-ray emitting material were derived in Ignace et al. (2000). Using these data, we found lower limits to the thickness of the spherical shells filled by hot gas to produce the observed level of X-ray emission. It was shown that the thickness of an isothermal spherical shell can be inferred from the analysis of the X-ray light curve and the observed filling factor.

Further, the model of an envelope consisting of a number of individual radially propagating X-ray emitting zones was developed by analogy with stochastic wind models. The total X-ray luminosity of such an atmosphere was obtained by using the formal solution of the radiation transfer equation for angle dependent optical depth. We did not place any constraints on the form, volume or density of the zones of hot gas except to demand that the total emission measure of each zone is the same as the emission measure of a spherical shell at the same distance from the stellar core, and that the zone must be relatively “small”. This allows us to restrict the number of independent parameters when performing the numerical simulations. The rationale for doing this was that broadband X-ray luminosities are known for most hot stars. Therefore, any combination of parameters describing the distribution of hot gas should provide the observed level of X-ray emission for each particular star as its first priority.

The two major simplifying assumptions of our modelling are those of isothermal X-ray sources and smooth cool material. We may speculate that including radiatively cooling sources will lead to a decrease in the level of variability. On the other hand, taking into account the clumpy structure of the background wind will lead to an increase of the variability of X-ray output. A detailed consideration of these effects are a matter of future work.

Despite its appoximative nature, the analysis described here revealed the basic properties of variability of X-ray emission for early type stars. Some of the main points are:

1. The apparent lack of short time-scale variability (order of an hour) of X-ray emission cannot be considered as a deficiency of shock models for X-ray production.
2. We may expect stochastic variability on long time-scales (thousands of flow times), especially for optically thick winds.
3. The level of X-ray variability depends on chemical composition and density of the stellar wind and differs for stars of different spectral classes even for similar mechanisms of X-ray production. It is governed by the opacity of the cool material and is substantially lower for the more opaque winds of WR stars and particularly WC stars.
4. The dependence of the wind opacity with energy means that the X-ray emission may be highly variable in the part of the spectrum, where the cool material is optically thin and practically constant in optically thick parts. Therefore the hard energies are specially apt for detection of variability of X-ray emission. Recall that existing claims of X-ray variability show differences between soft and hard passbands (e.g., possibly detected for $`\zeta `$ Ori by Berghöfer & Schmitt 1994).
5. Whether X-ray variability is detected or not for a given energy provides valuable information about the spatial distribution and properties of X-ray emitting material, and if variability is detected, its dependence on energy would be especially telling of the wind structure.

###### Acknowledgements.

This research was supported by a PPARC grant (LMO, JCB) and NASA grants NAG5-9964 (RI) and NAG5-9226 (JPC). We acknowledge discussions with K. Gayley, M. Hendry and S. Rauzy, and we thank the referee L. Koesterke for comments beneficial to this paper.


format: text

text:
# Mixed Lattice and Electronic States in High-Temperature Superconductors

\[

## Abstract

Inelastic neutron scattering measurements are presented which show the abrupt development of new oxygen lattice vibrations near the doping-induced metal-insulator transition in La2-xSrxCuO4. A direct correlation is established between these lattice modes and the electronic susceptibility (as measured by photoemission) inferring that such modes mix strongly with charge fluctuations. This electron-lattice coupling can be characterized as a localized one-dimensional response of the lattice to short-ranged metallic charge fluctuations.

\]

High-temperature superconductors are based on antiferromagnetic insulating materials caused by strong electronic correlations. Doping charge carriers (holes) into this system creates a two-dimensional correlated metallic state in the CuO2 plane that becomes superconducting at low temperatures. Much effort has been focused on the role of antiferromagnetic spin fluctuations in transport properties and superconductivity in the metallic state, mainly because it is believed that superconducting transition temperatures are too high to arise solely from electron-phonon coupling. Recently, this correlated metallic state has been discussed in terms of localized, atomic-scale charge density fluctuations. As hole doping affects mainly the hybridized Cu 3d$`_{x^2y^2}`$ and O 2p(x,y) anti-bonding electronic states, the lattice should couple strongly to a localized matrix of slow charge fluctuations, perhaps decisively influencing charge transport and superconducting properties. Inelastic neutron scattering measurements of the lattice dynamics do show evidence of strong and unusual electron-lattice coupling in high-Tc compounds with various crystal structures and methods of doping. However, many issues surrounding this anomalous coupling remain unclear. Of course, the ultimate issue concerns its role in superconductivity. Here, we report the systematic development of these anomalous modes on hole concentration and the abrupt formation of new oxygen lattice vibrations near the metal-insulator transition (MIT) in La2-xSrxCuO4. These new lattice modes may be associated with the bosons which interact with electronic states, causing a kink in the electronic dispersions observed by photoemission. These results are the best experimental evidence yet that the lattice is strongly mixed with charge dynamics in the high-temperature superconductors.

We performed inelastic neutron scattering measurements of phonon densities of states (DOS) in La2-xSrxCuO4 (LSCO) spanning hole concentrations from the undoped insulator to the optimally doped superconductor ($`0x0.15`$). The measurements were performed at $`T=10K`$ on the LRMECS spectrometer at the Intense Pulsed Neutron Source at Argonne National Laboratory. Sample preparation details and extraction of the DOS from the raw data are similar to those previously reported for La2-xSrxNiO4. The results of the present measurements are shown in Fig. 1. Of primary importance is the abrupt development of new lattice modes near $``$70 meV for hole concentrations between 6-8% ($`x=0.060.08`$). Because the DOS is not an analytic function and cannot be reliably fit to a series of peaks, this development is characterized by calculating the curvature of the DOS at 72 meV (inset of Fig. 1) which changes sign in this concentration range.

The abrupt change in the DOS is not due to the tetragonal-orthorhombic structural phase transition, as all samples remain in orthorhombic at $`T=10K`$. The new band is also not caused by electrostatic impurity effects from Sr substitution because identical measurements on the isostructural nickelate compounds show no such band formation at these Sr concentrations. Rather, the abrupt development of the 70 meV band must be related to the doping-induced MIT near $`x0.05`$, where trapped holes begin to become mobile. The onset of superconductivity and dynamic incommensurate spin fluctuations (signatures of the stripe correlations) also occurs in this critical concentration range. The nature of this MIT is unclear at present, however the transition does occur without a corresponding change in lattice symmetry. Thus, the 70 meV band signifies strong and unusual electron-lattice coupling in the metallic state of high temperature superconductors that is not present in the insulating state. The MIT also affects lower energy phonons, especially near 30 meV. Unfortunately, a plethora of phonon modes existing at this energy make analysis of the low energy DOS features difficult, and they are not discussed further.

The new lattice modes at 70 meV consist, at least partially, of half-breathing-like oxygen phonon modes that propagate in the CuO2 plane. Phonon dispersion measurements of LSCO by inelastic neutron scattering show that the half-breathing oxygen modes between $`\stackrel{}{q}=(\pi /2,0,0)`$ and $`(\pi ,0,0)`$ soften anomalously from 80 meV ($`x=0`$) to 70 meV ($`x=0.15`$) with doping.  Figure 2(a) demonstrates clearly that this flat half-breathing branch contributes to the DOS as a van Hove singularity for La1.85Sr0.15CuO4 near 70 meV. The half-breathing polarization is shown in Fig. 2(b).

To learn more about the nature of the electron-lattice coupling of the half-breathing modes, we use a simple lattice dynamics model to isolate the effective interaction being mediated by the mobile holes. We are only interested in the highest energy modes propagating within CuO2 plane. Also, for these vibrational energies (above 60 meV) the phonons are comprised entirely of in-plane polarized oxygen modes. Therefore, a two-dimensional harmonic force constant (ball and spring) model of the CuO2 plane is used to reproduce the high-frequency DOS. The initial set of pairwise force constants are chosen to reproduce the measured in-plane phonon dispersion and DOS of undoped La2CuO4. Doping is introduced by varying the magnitude of the force constants while maintaining the full periodicity of the CuO2 plane. Models studied with spatially varying force constants (superlattice models) produced only small variations in the DOS. The model reproduces well the subtle changes observed in the high frequency oxygen phonons within the insulating phase ($`x<0.06`$), such as the gradual weak softening of the $``$88 meV phonon band, by reducing the Cu-O nearest-neighbor force constant. This is shown in figure 3(a).

The model has difficulty reproducing the large changes of the oxygen phonons in the metallic state at $`x=0.08`$. One is required to introduce a repulsive force constant between next-nearest neighbor oxygens with a bridging copper to produce a band near 70 meV originating from the half-breathing mode. More importantly, this repulsion must be introduced anisotropically between only one oxygen pair in a given CuO4 plaquette (not the orthogonal pair), thereby breaking the crystal symmetry (although perhaps only locally and over phonon time scales). Introducing the repulsion over both pairs strongly softens oxygen breathing modes (near 85 meV) in addition to the half-breathing modes, inconsistent with experimental results. The calculated and measured differences in the $`x=0.06`$ and $`0.08`$ DOS are shown in figure 3(b) and results of the full DOS for $`x=0`$, $`0.06`$, and $`0.08`$ are shown in figure 3(c). This simple model suggests that the electron-lattice coupling is one-dimensional in nature and the coupled hole states are oriented along the Cu-O-Cu bond direction.

However, differences between the metallic and insulating DOS cannot be explained fully with the simple model. The intensity of the 70 meV band is relatively large implying that many phonons ($``$ 15% of all possible oxygen modes) are affected. The metallic phonon bands are also narrower in energy than the insulating bands, which suggests phonon localization. The sharpness and intensity of the metallic bands are consistent with phonon dispersion measurements in LSCO where the half-breathing-like modes are observed to be dispersionless in a large region of the Brillouin zone between $`\pi /2a<|q_x|<\pi /a`$, $`|q_y|<0.3\pi /a`$, and independent of $`q_z`$. This entire Brillouin zone pocket centered at $`(\pi ,0,0)`$ forms the rather intense van Hove singularity in the DOS at 70 meV. These observations cannot arise from the usual electron-phonon coupling in a metal which causes Kohn anomalies in the phonon dispersion at specific wavevectors $`\stackrel{}{q}=2\stackrel{}{k}_F`$ (where $`\stackrel{}{k}_F`$ is the Fermi surface wavevector).

Strong electron-lattice coupling and the one-dimensionality of the charge states are corroborated by various angle-resolved photoemission results. Bogdanov et al. show that the electronic band dispersion has a kink due mixing of holes with other bosonic excitations in the range of oxygen optical phonon energies. Lanzara et al.  conclusively identify these excitations as phonons. For LSCO, Lanzara et al. also show that strong electron-lattice coupling occurs with phonons in the 70 meV energy range, i.e. the half-breathing modes. Ino et al.  have studied the doping dependence of the photoemission in La2-xSrxCuO4. These results reveal that the $`(\pi ,0,0)`$ electronic saddlepoint sits $``$500 meV below the Fermi level in the insulating phase. Doping moves the saddlepoint close to the Fermi level, eventually producing a MIT at $`x0.05`$ where charge dynamics are characterized as one-dimensional. The maximum spectral weight of the saddlepoint band crosses the half-breathing phonon energy (70-80 meV) near $`x=0.07`$ consistent with new phonon modes being observed somewhat above the critical hole concentration of the MIT as measured by transport. This behavior is unclear, but may arise from the disparity in the time scales of inelastic neutron scattering and transport measurements.

It is possible that the one-dimensional charge fluctuations originate from states in the extended electronic saddlepoint near $`\stackrel{}{q}=(\pi ,0,0)`$. Then, half-breathing-like phonons and saddlepoint holes form localized states (or ”wave packets” made up of many plane waves). The stripe scenario supports this conjecture by assuming an inhomogeneous and localized charge distribution in the metallic state . Castro Neto has shown that strong phonon and electronic responses at $`(\pi ,0,0)`$ occur in the stripe model from the scattering of hole pairs at stripe domain walls. Another source of charge inhomogeneity originates from large Peierls-type (phonon induced) charge transfer fluctuations. The introduction of charge fluctuations into electronic calculations within the local-density approximation is known to produce strong screening of the half-breathing phonon modes in La2CuO4. Results obtained from exact diagonalization of the one-dimensional Peierls-Hubbard model also indicate that charge transfers caused by half-breathing modes are enhanced when strong correlations are included. These results suggest that the phonon softening in the metallic state is likely due to the strong local interaction of phonons with charge fluctuations above the metal-insulator transition. This important interaction bears consideration in any theory of high temperature superconductivity.

###### Acknowledgements.

RJM would like to thank T. Egami, A. Castro Neto and A. Ramirez for helpful discussions. This work was supported by the U.S. Department of Energy under contract number W-7405-Eng-36 with the University of California. This work has benefited from the use of the Intense Pulsed Neutron Source at Argonne National Laboratory. This facility is funded by the U. S. Department of Energy, BES-Materials Science, under contract W-31-109-Eng-38.


text:
# 1 Typical temperature dependence of the transmitted light intensity in the range of the orientational transition (thickness 𝑙=36.9 𝜇m)

Temperature Dependence of Anchoring Energy of MBBA on $`SiO`$-Evaporated Substrate

M.V. KHAZIMULLIN, A.P. KREKHOV and Y.A. LEBEDEV

Institute of Molecule and Crystal Physics, Russian Academy of Sciences, 450025 Ufa, Russia

The temperature behaviour of the anchoring energy of nematic liquid crystal MBBA at the substrates with oblique $`SiO`$ evaporation and director orientation at the interface have been studied by means of magnetic Fréedericksz transition technique. The temperature dependence of the coefficients in the phenomenological model of surface anchoring energy has been determined.

Keywords: nematic liquid crystal, anchoring energy, orientational transition

INTRODUCTION

The equilibrium director orientation in a nematic liquid crystal (NLC) layer strongly depends on the anchoring properties of the bounding substrates. In the absence of an external field the confining surfaces impose the orientation in the bulk of nematic layer. If an external field is present the director distribution results from the competition between the surface and bulk torques. Study of the orientational behaviour in an external field provides the information about anchoring properties of NLC at the bounding substrates . Besides, a useful tool for the investigation of interfacial phenomena is the temperature-induced orientational transitions which have been found at some interfaces .

Recently the temperature-induced orientational transitions planar $``$ tilted $``$ homeotropic close to the clearing point (nematic - isotropic transition) have been found at the interface between MBBA and a $`SiO`$-evaporated substrate . To describe this orientational behaviour the phenomenological model of surface energy similar to  taking into account the substrate micro relief has been suggested.

In this paper we present the results of experimental investigations of the orientational behaviour of the nematic liquid crystal MBBA in a magnetic field on the same interface as in  in the temperature range where the orientational transitions occurred. The director tilt angle at the substrate as a function of temperature has been measured and compared with theoretical predictions. The parameters of the phenomenological model are determinated and the temperature dependence of surface energy is found.

BASIC EQUATIONS

We consider a nematic liquid crystal layer of thickness $`l`$ confined between two identical plates at $`z=0`$ and $`z=l`$. Weak anchoring is provided by the substrates with micro relief (grooved surface). The substrate with micro relief is characterised by the grooves axis (along $`x`$ axis) and normal vector. Then the phenomenological expression for the surface anchoring energy $`f_s`$ for such substrate can be obtained using the expansion in terms of the nematic tensor order parameter $`Q_{ij}=S(T)[n_in_j1/3\delta _{ij}]`$ up to second order in the scalar order parameter $`S(T)`$ (see  for details)

$$f_s=C(S)\mathrm{cos}^2\theta _0+D(S)\mathrm{cos}^4\theta _0,$$
(1)
with

$$C(S)=aS+bS^2,D(S)=dS^2,$$
(2)
where $`\theta _0`$ is the angle between the director at the substrate and the $`x`$ axis and $`a`$, $`b`$, $`d`$ are the temperature independent model parameters which depend on the micro relief properties and the nematic - substrate interaction. If there is no external field the director distribution in the nematic layer is defined by the surface angle $`\theta _0`$, which can be found by minimisation of the surface energy $`f_s`$ with respect to $`\theta _0`$. Depending on the model parameters ($`a`$, $`b`$ and $`d`$) one could have different scenarios for the temperature behaviour of the surface angle $`\theta _0`$. For example, $`a=0`$, $`b<0`$, $`d>0`$, $`b+2d>0`$ describes a temperature independent pretilt, or, if $`a>0`$, $`b<0`$, $`d>0`$, $`a+b<0`$, $`a+b+2d>0`$, one has a tilted $``$ homeotropic transition similar to that found in . Here we focus on the case of $`a>0`$, $`b<0`$, $`d>0`$, $`a+b<0`$, $`a+b+2d<0`$ when the temperature-induced orientational transition planar $``$ tilted $``$ homeotropic takes place. Minimisation of surface energy (1) gives the temperature behaviour of surface angle

$$\mathrm{cos}^2\theta _0=\frac{S_p}{S_pS_h}(1\frac{S_h}{S}),$$
(3)
with transition points planar $``$ tilted ($`S_p`$) and tilted $``$ homeotropic ($`S_h`$)

$$S_p=\frac{a}{b+2d},S_h=\frac{a}{b}.$$
(4)
The temperature dependence of the surface angle $`\theta _0`$ can be found from Eq.(3) using the experimental data for the scalar order parameter $`S(T)`$ .

In the case when a magnetic field $`𝐇`$ is applied along the $`z`$ axis, the director distribution can be obtained by minimisation of the total free energy (per unit area)

$$F=\frac{1}{2}_0^l\left[K(\theta )\left(\frac{d\theta }{dz}\right)^2\chi _aH^2\mathrm{sin}^2\theta \right]𝑑z+f_{s1}+f_{s2},$$
(5)
where $`K(\theta )=K_1\mathrm{cos}^2\theta +K_3\mathrm{sin}^2\theta `$, $`K_1`$ and $`K_3`$ are elastic constants for “splay” and “bend” deformations respectively, $`\chi _a`$ is the diamagnetic anisotropy and $`f_{si}`$ are the surface energies for the lower ($`i=1`$) and upper ($`i=2`$) substrates. After the standard procedure one gets the equation for the director profile (first integral of Euler-Lagrange equation)

$$\frac{d\theta }{dz}=H\sqrt{\chi _a\frac{\mathrm{cos}^2\theta \mathrm{cos}^2\theta _m}{K(\theta )}}$$
(6)
and boundary condition (in the case of identical substrates $`f_{s1}=f_{s2}=f_s`$)

$$\frac{df_s}{d\theta _0}=H\sqrt{\chi _aK(\theta _0)(\mathrm{cos}^2\theta _0\mathrm{cos}^2\theta _m)},$$
(7)
where $`\theta _m`$ is the angle in the midplane of the nematic layer and the symmetry of the solution with respect to the midplane has been used ($`d\theta /dz=0`$ at $`z=l/2`$). Integrating (6) from $`\theta _0(z=0)`$ to $`\theta _m(z=l/2)`$ gives the relation between $`\theta _m`$ and $`\theta _0`$ by

$$\frac{\pi }{2}\frac{H}{H_F}=_{\theta _0}^{\theta _m}P(\theta )𝑑\theta ,P(\theta )=\sqrt{\frac{1+\eta \mathrm{cos}^2\theta }{\mathrm{cos}^2\theta \mathrm{cos}^2\theta _m}},$$
(8)
where $`H_F=(\pi /l)\sqrt{K_1/\chi _a}`$ is the Fréedericksz transition field for strong planar anchoring and $`\eta =K_3/K_11`$.

Using expressions (1), (2), (4) and (7) the temperature behaviour of the surface angle $`\theta _0`$ in the presence of a magnetic field can be derived in the following form

$`\mathrm{sin}(2\theta _0)bS\{(SS_h)(1{\displaystyle \frac{S_h}{S_p}})S\mathrm{cos}^2\theta _0\}`$
$`=H\sqrt{\chi _aK(\theta _0)(\mathrm{cos}^2\theta _0\mathrm{cos}^2\theta _m)}.`$ (9)
The values of $`S_p`$, $`S_h`$ can be easily found from the experimental data on the planar $``$ tilted $``$ homeotropic transition in zero magnetic field. Measuring the temperature dependence of the surface angle $`\theta _0`$ at different values of the magnetic field one can verify the phenomenological model (1), (2) and determine the model coefficients (which should be independent of the magnetic field).

EXPERIMENTAL

The experimental cell consists of two glass plates with mylar spacers of thickness $`35`$ $`\mu `$m filled by nematic liquid crystal MBBA. The inner surface of the confining plates was covered by thin layer of $`SiO`$, which was vacuum evaporated under the angle $`60^{}`$ with respect to the surface normal. At room temperature one obtains homogeneous planar orientation (along $`x`$ axis). The cell was mounted in the hot stage and demonstrated the temperature-induced orientational transition planar $``$ tilted $``$ homeotropic under heating. The accuracy of the temperature measurements was better than $`10`$ mK. To study the orientational behaviour of NLC in a magnetic field the experimental setup consisting of an electromagnet with maximum field $`12`$ kOe and $`12`$ mm gap between the poles has been designed. The stability of the magnetic field was $`\pm 10`$ Oe. The intensity of transmitted light as a function of temperature $`I(T)`$ at the fixed values of magnetic field has been measured by a photometer MPM-100 (Zeiss) (light source: He-Ne laser, $`\lambda =632`$ nm). All measured signals (magnetic field, temperature, and transmitted intensity) were processed by a computer ADC-card.

For the case of cross polars and $`x`$ axis at $`45^{}`$ to the polarisers one has for the transmitted light intensity

$$I=I_0\mathrm{sin}^2\frac{\delta }{2},\delta =\frac{2\pi }{\lambda }_0^lR(\theta )𝑑z$$
(10)
with

$$R(\theta )=\frac{n_on_e}{\sqrt{n_o^2\mathrm{cos}^2\theta +n_e^2\mathrm{sin}^2\theta }}n_o.$$
(11)
Here $`\delta `$ is the phase difference; $`\lambda `$ is wavelength of light; $`n_o`$, $`n_e`$ are the ordinary and extraordinary refractive indices. Since $`\delta =\delta (T)`$ is monotonically decreasing to zero during planar $``$ tilted $``$ homeotropic transition, it can be easily obtained from the experimental data on the temperature dependence of the transmitted light intensity. Then, taking into account (6), (10) and (11) one obtains the relation between the phase difference and surface angle

$$\delta =\frac{4l}{\lambda }\frac{H_F}{H}_{\theta _m}^{\theta _0}R(\theta )P(\theta )𝑑\theta .$$
(12)
Using the temperature dependence of the material parameters $`K_1`$, $`K_3`$ and $`\chi _a`$ of MBBA , expressions (12), (8) allow to calculate the temperature dependence of the surface angle $`\theta _0(T)`$ from the experimental data on the phase difference $`\delta (T)`$ for the different values of magnetic field.

RESULTS

In Fig. 1 the typical dependence of the transmitted light intensity on the reduced temperature $`\tau =T/T_c`$ for zero magnetic field close to the clearing point is shown. The temperature dependence of the surface angle $`\theta _0`$ derived from these experimental data is plotted in Fig. 2 (triangles). The corresponding transition points planar $``$ tilted ($`\tau _p`$) and tilted $``$ homeotropic ($`\tau _h`$) are shown at Fig. 1 by arrows. The dependence $`\theta _0(\tau )`$ for zero magnetic field has been fitted by Eq.(3) and the transition points have been determined: $`\tau _p=0.9961`$ and $`\tau _h=0.9984`$.

In the presence of a magnetic field the surface angle deviates from the initial value $`\theta _0=0`$ and then monotonically increases up to $`\theta _0=\pi /2`$ with increasing temperature (Fig. 2). Note, that the temperature of the tilted $``$ homeotropic transition is decreased with increasing magnetic field. The solid lines in Fig. 2 show the behaviour of the surface angle calculated from the theoretical model for corresponding magnetic fields. Using in Eq.(S0.Ex1) the values of $`\tau _p`$ and $`\tau _h`$ found for zero magnetic field, one obtains $`b=0.055`$ erg/cm2 by fitting to the experimental data for $`\theta _0(\tau )`$ for different values of magnetic field. In a wide range of magnetic field in the experiments the coefficient $`b`$ was found to be independent of $`H`$.

Finally, in Fig. 3 the temperature dependence of the coefficients $`C(\tau )`$, $`D(\tau )`$ and $`C(\tau )+2D(\tau )`$ are presented together with the surface energy $`f_s`$ \[Eq.(1)\]. One sees that $`C(\tau )`$ changes sign at the tilted $``$ homeotropic transition point ($`\tau _h`$) and $`C(\tau )+2D(\tau )`$ at the planar $``$ tilted transition point. The surface energy $`f_s`$ is quite small and decreases to zero at the clearing point. At temperature $`T=26`$ C one gets for the surface energy $`810^3`$ erg/cm2.

It should be noted, that MBBA is a rather unstable substance, which can have an influence on the determination of model coefficients. For estimation of possible errors the coefficients $`C`$, $`D`$ were calculated taking into account the uncertainties in the experimental data and material parameters. Varying the material parameters in the range of $`\pm 10\%`$ results in $`20\%`$ uncertainty in the anchoring energy. The influence of inaccuracies in the measured quantities (magnetic field, temperature, thickness etc.) was essentially smaller ($`3÷5\%`$).

Thus, the temperature-induced orientational transitions planar $``$ tilted $``$ homeotropic at the MBBA - $`SiO`$-evaporated substrate has been investigated. From the experimental data on the transmitted light intensity in the presence of a magnetic field the temperature dependence of the surface angle has been found and the parameters of the phenomenological model of surface energy have been determined. The results demonstrate a good quantitative agreement between the experimental data and the suggested surface energy model.

Acknowledgments

We thank O. Khasanov for the help in setup construction and substrates preparation and L. Kramer for fruitful discussions and critical reading of the manuscript. A.K. is grateful to the Alexander von Humboldt Stiftung for equipment donation. Financial support from INTAS Grant No. 96-498 is gratefully acknowledged.


text:
# BES R measurements and 𝐽/𝜓 decays

## 1 BES R Measurements

R is defined as,

$$R=\frac{\sigma (e^+e^{}hadrons)}{\sigma (e^+e^{}\mu ^+\mu ^{})}=3\underset{q}{}Q_q^2$$

A precision measurement of R will reduce the uncertainty on $`\alpha (M_z)`$, narrow the predicted mass window of Higgs boson $`m_H`$ and reduce the uncertainty on the anomalous magnetic moment of muon $`a_\mu `$. R has been measured in the energy region from hadron production threshold to $`Z^0`$ pole. Previous R measurements below 5 Gev were carried out almost 20 years ago and the largest uncertainty on R comes from the measurement below 5 GeV energy region. Therefore, it is important to measure R value in the 2 – 5 GeV energy region. BES collaboration has carried out two runs of R scan. 6 energy points were scanned in the 1998’s run. Data samples at 2.6 GeV and 3.55 GeV were used to tune LUND parameters. In the second run of 1999, 85 energy points were scanned, from which 24 points with separated-beam data and 7 points with sigle-beam data. In the experiment, R is determined by the following realtion:

$$R=\frac{\sigma ^0had}{\sigma _{\mu \mu }^0}=\frac{N_{had}^{obs}N_{bg}_lN_{ll}N_{\gamma \gamma }}{\sigma _{\mu \mu }^0\epsilon _{had}\epsilon _{trg}(1+\delta )},$$

where $`N_{had}^{obs}`$ represents the observed hadronic events, $`N_{bg}`$ represents the beam-associated backgrounds, $`N_{ll}`$ represents the lepton pair backgrounds, $`N_{\gamma \gamma }`$ represents two photon process events, $`\epsilon _{had}`$ stands for the detection efficiency, $`\epsilon _{trg}`$ is the trigger efficiency, $`(1+\delta )`$ is the radiative correction, $``$ is the integrated luminosity and $`\sigma _{\mu \mu }^0`$ is the lowest cross section of $`e^+e^{}\mu ^+\mu ^{}`$ which is given by $`\sigma _{\mu \mu }^0(s)=4\pi \alpha ^3/3s`$. One of the most important work in R measurements is hadronic event selection. It is done in three steps. First remove apparent Bhabha events. Then in track level, select track to remove cosmic rays, beam-gas. Finally, in event level, remove beam-gas, two-photon, $`ee`$, $`\mu \mu `$, $`\mathrm{}`$. The main backgrounds come from cosmic rays, lepton pair production, two-photon processes and beam associated background. The most serious background is beam associated background. It is studied in two independent methods. The first method is to use separated-beam data to calculate the event number of beam associated background, the second method is to use z component of event vertex to calculate the event number. The results given by two different method are consistent. Detection efficiency is given by Monte Carlo simulation, which is model dependent. However, results given by different models are consistent in 3 percent level. The integrated luminosity is determined by $`L=\frac{N_{obs}}{\sigma \epsilon \epsilon _{trg}},`$ with $`N_{obs}`$ the number of events (Bhabha, dimuon and $`\gamma \gamma `$), $`\sigma `$ the production cross section, $`\epsilon _{trg}`$ the trigger efficiency and $`\epsilon `$ the detection efficiency. Initial state radiative correction is studied in four different schemes and they are consistent within $`1\%`$ in continuum region. Finally, Crystall Ball scheme is selected to calculate the radiative correction. Uncertainties of R values are reduced down to $`610\%`$. The average uncertainties of R is about $`6.6\%`$. As an example, the error source of R value at 3.0 GeV energy point is shown in Table 1. R values in 1 – 5 GeV energy region is shown in Fig.1.

## 2 Search for $`\sigma `$ in $`J/\psi `$ Hadronic Decay

The early analysis of $`\pi \pi `$ and $`\pi K`$ scattering data shows no pole at the lower mass region. Up to now, the existence of $`\sigma `$ and $`\kappa `$ as resonant particles has not been widely accepted. However, recent re-analysis of the $`\pi \pi `$ and $`\pi K`$ scattering data shows an evidence for existence of the $`\sigma `$ and $`\kappa `$ particle with comparatively light mass. If $`\sigma `$ is a s-channel resonance, it should appear in production process. Therefore, it is important to search for $`\sigma `$ and $`\kappa `$ in production process.

One found that there is a low mass enchancement in the $`\pi ^+\pi ^{}`$ invariant mass spectrum in $`J/\psi \omega \pi \pi `$ (Fig.2). For the sake of complicity, we call this low mass enhancement the $`\sigma `$-particle. In Daliz plot, a clear band which corresponds to the $`\sigma `$-particle can be seen(Fig.2). There are some backgrounds in the above data sample (such as $`J/\psi \rho 3\pi `$, $`\mathrm{}`$). Some of the backgrounds are shown in the $`\omega `$ side-band, which can be removed through side-band subtraction. $`\sigma `$-particle can still be observed in the $`\pi ^+\pi ^{}`$ invariant mass spectrum after side-band subtraction (Fig.2). Threshold effects and phase space effects can also be seen in this figure. Apparently, the first peak does not originate from threshold effect and phase space effect. We also performed detailed Monte Carlo study and found that the background events from other $`J/\psi `$ decay channels are very less. Therefore, the first peak is likely a s-channel resonance. Then, a Partial Wave Analysis(PWA) on the $`\pi ^+\pi ^{}`$ invariant mass spectrum is applied to study the structure of the first peak. Our results strongly favor that the spin-parity is $`0^{++}`$ and the statistical significance of $`\sigma `$ is about 18 $`\sigma `$. Its mass and width are determined through mass and width scan and the results are listed in Table 2. Our scan curve is shown in Fig.3 with the cross representing real data and the histogram being fit curve. The final fit on the global mass spectrum is also shown in Fig.3. $`f_0(980)`$ also appears in this channel. Our results on $`f_0(980)`$ is listed in Table 2. The branching ratios listed there are only for $`J/\psi \omega X\omega \pi ^+\pi ^{}`$.

$`\sigma `$, as a s-channel resonance, is not found in the $`J/\psi \gamma \pi ^+\pi ^{}`$ It means that it does not look like a scalar glueball. According to quark model, the ordinary $`q\overline{q}`$ scalar meson nonet with lowest mass is $`1^3P_0`$ states. Its orbital angular momentum is excited to P wave. The mass of them can not be as low as 390 MeV. $`\sigma `$ can not be filled into the quark model for ordinary $`q\overline{q}`$ mesons.


text:
# On a certain analogy between hydrodynamic flow in porous media and heat conductance in solids

## I Introduction

Imposing an inhomogeneous pressure/temperature field on a (porous) material causes in general mechanical deformations and stresses. These ‘induced’ elastic deformations in turn can and do influence the further mass/heat transport. The problem is described as a coupled field problem consisting of (at least) two equations (a) an elastic equation describing mechanical displacements due to a given source field (pressure/temperature gradient field) and (b) at least one transport equation for mass/heat with a source term originating from the mechanical deformation (divergence of the displacement field).

The entire description bases on continuity equations and constitutive laws which are discussed in textbooks in great detail. Perhaps less known is Biot’s theory of hydrodynamic flow in porous elastic media. We will address this problem below in more detail. In particular we wish to draw some attention on similarities for the flow and the heat problem.

## II General considerations

We give the basic outline of Biot’s theory for the particular case of an incompressible pore fluid. This considerably simplifies the governing equations and their discussion, because all material derivaties can be eliminated. First we discuss the basic elastic equations. This is followed by a detailed discussion of the flow equation for incompressible pore fluids. Thereafter we discuss some properties of the coupled equations, and finally establish their general integral following the spirit of Galerkin’s method. We specialize the obtained general integral for the case of irrotational displacement fields and make the connection to the problem of thermoelasticity.

### A Preliminaries

In the following bold-face letters always indicate vectors/tensors.

#### 1 Lamé equation

The continuity equation for momentum density,

$$div\sigma +𝐟=\rho _s\ddot{𝐮},$$
(1)
represents the mechanical equilibrium condition, where $`\sigma =(\sigma _{ij})`$ is the stress tensor, $`𝐟`$ any external body force, $`\rho _s`$ the solid’s mass density and $`𝐮=(u_i)`$ the mechanical displacement. We will neglect the inertial terms in the following and will consider Biot’s theory of porous media for the case of symmetric elasticity, i.e. $`\sigma _{ij}=\sigma _{ji}`$. Deformations are characterized by the strains $`2u_{ij}=(u_i/x_j+u_j/u_i)`$. Employing Hooke’s linear relation between stresses and strains for isotropic materials one obtains Lamé’s equation for the displacements,

$$(12\nu )\mathrm{\Delta }𝐮+graddiv𝐮+\frac{2}{E}(1+\nu )(12\nu )𝐟=\mathrm{𝟎},$$
(2)
with Poisson contraction number $`\nu `$ and Young modulus E. There are a few quite general results about this equation, which will be useful in the following. Taking the divergence of Eq.(2) one has,

$$\mathrm{\Delta }div𝐮=\frac{(1+\nu )(12\nu )}{(1\nu )E}div𝐟,$$
(3)
i.e., $`div𝐮`$ is the solution of a Poisson equation and for vanishing body forces $`𝐟`$ it is even a potential function. The physical significance of $`div𝐮`$ is that it measures the local relative volume change due to the deformation which is an invariant. Pure shear deformations are unreflected in Eq.(3) because they preserve the volume. Of course all information on vorticity has been lost by passing from Eq.(2) to Eq.(3). On the other hand taking the Laplacian on Eq.(2) and inserting Eq.(3) into yields,

$$\mathrm{\Delta }\mathrm{\Delta }𝐮=\frac{2(1+\nu )}{E}\left(\frac{1}{2(1\nu )}graddiv𝐟\mathrm{\Delta }𝐟\right).$$
(4)
This equation is of fourth order while Eq.(2) is of second order; hence their solutions will be in general not identical. In his seminal work Galerkin resolved this problem. He introduced another vector function $`𝐅`$ satisfying the (inhomogeneous) biharmonic equation,

$$\mathrm{\Delta }\mathrm{\Delta }𝐅=\frac{2}{E}(1+\nu )𝐟,$$
(5)
with $`𝐟`$ being the same body force as in Eq.(2). He showed that the integral of Eq.(2) can be expressed in form of derivatives of $`𝐅`$,

$$𝐮=\mathrm{\Delta }𝐅\frac{1}{2(1\nu )}graddiv𝐅.$$
(6)
We would like to note that the foregoing considerations made no particular assumptions on the body force field $`𝐟`$ in particular they apply to conservative fields. Typically $`𝐟`$ is explicitly known like for gravitational forces, or it is itself a solution of another independent equation. The latter is the case for thermo-elastic problems, in which

$$𝐟=\alpha KgradT$$
(7)
where $`\alpha `$ is the heat expansion coefficient, the $`K=\frac{E}{3(12\nu )}`$ accounts for the material’s compressibility, and $`T`$ is the temperature field. Typically Eqs.(2) and (7) are supplemented by a heat conductance equation,

$$\chi \mathrm{\Delta }T=\frac{}{t}T,$$
(8)
with $`\chi `$ being the coefficient of heat conductance divided by the specific heat capacity. Equation (5) then becomes a biharmonic equation with time dependent –but known inhomogeneitis– which do not dependent on the deformation themselves. As will be shown further below does flow in porous materials imply in general ‘body forces’ being deformation dependent themselves. The deformation acts back on the flow. We will come back to this point in more detail further below.

#### 2 Flow equation

Lets consider fluid flow in a porous continuum. The porosity (local volumetric fluid fraction) can be described by a scalar field $`\varphi (𝐫,t)`$. It is common to assume that the relative fluid-solid velocity $`\varphi (𝐯\dot{𝐮})`$ can be expressed by Darcy’s relation,

$$𝐟=gradp=\varphi \frac{\kappa }{\mu }(𝐯\dot{𝐮}),$$
(9)
where $`𝐯`$ and $`\dot{𝐮}`$ are the fluid and solid velocities, $`\kappa `$ the mechanical permeability, $`\mu `$ the fluid viscosity, and $`p`$ the hydrostatic pressure field.

The continuity equations for solid and fluid mass are,

$$\frac{}{t}\left((1\varphi )\rho _s\right)+div\left((1\varphi )\rho _s\dot{𝐮}\right)=0,$$
(10)
and

$$\frac{}{t}\left(\varphi \rho _l\right)+div\left(\varphi \rho _l𝐯\right)=0,$$
(11)
respectively, where $`\rho _s`$ and $`\rho _l`$ denote the solid and fluid mass densities. In the following we make the two important assumptions that the mass densities are constants ($`\rho _s=const.`$ and $`\rho _l=const.`$).

From Eqs.(9), (10) and (11) and the assumption of constant mass densities follows a) that the local volume is conserved because the inflow of fluid volume into a volume element equals the solid volume outflow, and b) that the flow equation takes the following simple form ($`\kappa /\mu =const`$),

$$\frac{\kappa }{\mu }\mathrm{\Delta }p=div\dot{𝐮},$$
(12)
i.e. a Poisson equation for the pressure with the rate of the relative elastic volume change as source term. The porosity $`\varphi `$ does not enter explicitly into the flow equation because of constant mass densities. As soon as $`\dot{𝐮}`$ has been determined from the coupled Eqs.(2) and (12) $`\varphi `$ can be computed from Eq.(10).

### B The coupled system

#### 1 The Galerkin-Biot potential

After the above preparations we will focus on the coupled Eqs.(2) and (12) describing elastic–flow interactions for incompressible pore fluids. Note that

$$𝐟=gradp,$$
(13)
in Eq.(2), comp. Eq.(9). In general Lamé’s vector equation (2) is much harder to solve than the scalar Poisson Equation (12). However, for the very particular case of irrotational displacements Lamé’s equation can be integrated easily, see Sec. III A. In general (and also for most boundary conditions of practical interest) displacements and velocities will have a non-zero rotational part. Therefore $`div𝐮`$ or its time derivative will not contain the full information about displacements and strains, i.e. all elastic volume preserving deformations (shears) are unreflected in it.

It has been convincingly demonstrated that for heat-conductance in solids the relative volume changes, $`div\dot{𝐮}`$, actually do also appear in Eq.(8). However, the corresponding coefficient is for most solids so small that this term can be safely disregarded simplifying the theoretical treatment considerable. As an analogy to heat-conductance in solids one can think of a strongly compressible fluid, i.e. a gas in a porous medium. One can neglect the back-influence of the pressure induced elastic distortion on the pore pressure field and is basically left with the equations decribing a corresponding thermo-elastic problem. On the other hand does flow of an nearly incompressible fluid in a porous medium represent the ‘opposite’ case: in order flow to happen elastic deformation must take place. Without elastic deformation there is no flow and vice versa.

We reconsider the flow problem Eq.(5) and make some small refinements. Of course the pore pressure field is unknown in Eq.(5) since we have not incorporated any information from Eq.(12). This needs to be done in order to obtain a fully decoupled equation for displacements. Starting from Eq.(12) we have,

$$𝐟=gradp=\frac{\mu }{\kappa }grad_V^{}\frac{div\dot{𝐮}(𝐫^{})}{4\pi |𝐫𝐫^{}|}𝑑V^{}.$$
(14)
This equation still needs to be complemented by boundary conditions for the pressure field. Inserting Eq.(6) into Eq.(14) and integrating employing Helmholtz’s theorem one obtains a relation between $`𝐟`$ and $`𝐅`$ which is in turn inserted in Eq.(5),

$$\mathrm{\Delta }\mathrm{\Delta }𝐅=\frac{1}{D}graddiv\dot{𝐅},$$
(15)
with the abbreviation

$$D=\frac{\kappa }{\mu }\frac{E(1\nu )}{(1+\nu )(12\nu )}.$$
(16)
Equation (15) is the main general result in this section. It needs to be complemented by boundary and initial conditions for displacements and pressure. It is not difficult to see from Eqs.(6) and (15) that if $`curlu0`$ the displacements are characterized by a vector-diffusion equation. However, the general case with non-vanishing curl is much more complicated (it is described by an integro-differential equation in $`𝐮`$).

#### 2 Dispersion relation

While theoretical simplifications within models are usually appreciated a remaining question is whether the simplified model assumptions do lead to a consistent physical picture concerning certain criteria. One such criterion is stability. For example, is it obvious that the model equations are stable if one disregards inertial forces for the elastic continuum phase? This problem is most convenient considered investigating the dispersion relation $`\omega (𝐤)`$ obtained from an expansion after eigenfunctions. For brevity we only present the basic result in the following.

Assuming certain mathematical integrability conditions on displacements $`𝐮(𝐫,t)`$ and pressure $`p(𝐫,t)`$ both fields can be represented as Fourier-Integrals,

$$𝐮(𝐫,t)=𝑑𝐤\stackrel{~}{𝐮}(𝐤)e^{i(\omega t\mathrm{𝐤𝐱})},$$
(17)
and

$$p(𝐫,t)=𝑑𝐤\stackrel{~}{p}(𝐤)e^{i(\omega t\mathrm{𝐤𝐱})}.$$
(18)
Inserting this two representations into the basic Eqs.(2) and (12) one obtains the dispersion relation,

$$\omega (𝐤)=iD𝐤^2,$$
(19)
with $`D`$ being already defined in Eq.(16) The dispersion relation has the form of a diffusion only problem, i.e., all solutions are damped in time and therefore the equations are stable. The damping is controlled by a constant of diffusion $`D`$ in a rough estimate of order $`10^1\frac{m^2}{s}`$ which is orders of magnitude larger than for molecular diffusion (we estimated $`D`$ from $`\kappa 10^{15}m^2`$, $`E10^{11}Pa`$, $`\mu 10^3Pas`$ and $`\nu =0.25`$).

## III The point source <br>hydraulic loading)

In the following we consider a particular simple case in three dimension: at time $`t_0`$ the material is loaded by a pressure peak $`\delta (t)\delta (𝐫)`$. Dimensionless quantities for the length scale, stress and time scale will be useful (denoted by a star);

$$𝐫^{}=\frac{𝐫}{\kappa ^{1/2}};p^{}=\frac{p}{G};t^{}=t\frac{G}{\mu },$$
(20)
with $`\kappa `$, $`2G=E/(1+\nu )`$ and $`\mu `$ being the permeability, elastic modulus and fluid viscosity respectively.

### A Pressures and Displacements

Because the chosen symmetry of boundary conditions guarantuees the existence of an irrotational solution Lamé’s equation Eq.(2) simplifies considerable,

$$div^{}𝐮^{}=a_{}^{}{}_{}{}^{2}(p^{}p_0^{}),$$
(21)
where $`p_0^{}`$ denotes the equilibrium pressure and $`2a_{}^{}{}_{}{}^{2}=(12\nu )/(1\nu )`$ a material constant. Note that the solution for $`𝐮^{}`$ can be obtained by direct integration of Eq.(21),

$$𝐮^{}(𝐫^{})=\frac{a_{}^{}{}_{}{}^{2}}{4\pi }grad^{}_V^{}𝑑V_{}^{}{}_{}{}^{}\frac{p^{}(𝐫_{}^{}{}_{}{}^{})p_0^{}}{|𝐫^{}𝐫_{}^{}{}_{}{}^{}|}.$$
(22)
Equation (21) can be inserted into the pressure equation (12) giving an (inhomogeneous) heat conductance equation,

$$\mathrm{\Delta }^{}p^{}a_{}^{}{}_{}{}^{2}\frac{p^{}}{t^{}}=\delta (t^{})\delta (𝐫^{}),$$
(23)
with three dimensional solution (Green’s function),

$$p^{}(𝐫^{},t^{})p_0^{}=\frac{1}{a_{}^{}{}_{}{}^{2}}(A/\pi )^{3/2}e^{Ar_{}^{}{}_{}{}^{2}},$$
(24)
where,

$$A=\frac{a_{}^{}{}_{}{}^{2}}{4t^{}}=\frac{(12\nu )\mu }{8(1\nu )Gt}.$$
(25)
Thus incompressible flow in a porous medium whose displacement field is irrotational corresponds to a heat conductance problem in a fiktive non-porous medium. One can insert the pressure field Eq.(24) into Eq.(22) and evaluate the integral,

$$𝐮^{}(𝐫^{})=\frac{1}{2\pi ^{3/2}}\left(\frac{1}{𝐫_{}^{}{}_{}{}^{2}}_0^{\sqrt{A}|𝐫^{}|}e^{z^2}𝑑z\frac{\sqrt{A}}{|𝐫^{}|}e^{A𝐫_{}^{}{}_{}{}^{2}}\right)𝐞_r.$$
(26)
In Fig.1 we show the dimensionless radial displacements, $`u_r^{}`$, for two different times.

The displacements being zero at the symmetry center raise to a maximum value and drop for larger radii towards zero again. As can be seen from Eq.(26) do the displacements follow a scaling form. The scaling variable $`\delta =\sqrt{A}𝐫^{}`$ may be anticipated for a diffusion-related problem.

### B Deformations

It is obvious that the deformations can be written in scaling form,

$$u_{rr}^{}=\left(\frac{A}{\pi }\right)^{3/2}H(\delta ),$$
(27)
with $`H(\delta )`$ being a scaling function of the form,

$$H(\delta )=e^{\delta ^2}+\frac{e^{\delta ^2}}{\delta ^2}\frac{1}{\delta ^3}_0^\delta 𝑑ze^{z^2},$$
(28)
being plotted in Fig. 2.

Let $`\delta _0`$ denote the root of the forgoing equation. We find,

$$\delta _0=\sqrt{A}r_0^{}1,$$
(29)
within three percent error. The volumes exhibiting positive and negative stresses are therefore separated by a surface of radius $`r_0^{}`$ at time $`A`$ around the center, or equivalently $`r_0=2\sqrt{Dt}`$ with diffusity $`D`$ defined in Eq. (16), eg. an order of magnitude estimate gives $`r_020cm`$ for $`t=0,1sec`$ for the neutral zone. Another interesting property is the asymptotic behaviour of the deformation for $`\delta \mathrm{}`$. The asymptotics is given as,

$$\underset{\delta \mathrm{}}{lim}\frac{u^{}}{r^{}}=\frac{1}{2\pi r_{}^{}{}_{}{}^{3}}=\frac{1}{2\pi }\left(\frac{\kappa ^{1/2}}{r}\right)^3.$$
(30)
It is interesting to note that neither elastic constants nor fluid properties enter this relation. Only the permeability enters this compressional far-field decay. This relation should be compared to the simple case of a cavity of radius $`R`$ bearing pressure $`p`$ within an infinite non-Biot medium with vanishing pressure at infinity,

$$u_{rr}=p\frac{(1+\nu )}{E}\left(\frac{R}{r}\right)^3=\frac{p^{}}{2}\left(\frac{R}{r}\right)^3,$$
(31)
when one measures the bearing pressure in units of the shear modul $`G=\frac{E}{2(1+\nu )}`$. Both relations are essentially the same, where the typical length-scale $`\kappa ^{1/2}`$ in Eq.(30) given by the solids permeability acts similar to the cavity-size in non-Biot elasticity. In principle this should allow to measure the permeability experimentally for $`\delta \mathrm{}`$, which means in practice for distances $`r>>2\sqrt{Dt}`$. This observation leads us to the conclusion that the more complicated ’Biot elasticity’ reduces for length scales much larger than the scale of diffusion towards convential linear elasticity.

The maxima of the tensile deformations, $`\delta 0`$, are readily calculated using $`lim_{\delta 0}H(\delta )=1/3`$ ,

$$\underset{\delta 0}{lim}\frac{u^{}}{r^{}}=\frac{A^{3/2}}{3\pi ^{3/2}}.$$
(32)
This value does not depend on the material’s permeability.

Before turning towards the mechanical stresses we quickly summarize the results for the polar deformations $`u_{\varphi \varphi }^{}`$. Fig. 3 shows the polar deformations. in a rescaled form.

One has,

$$u_{\varphi \varphi }^{}=\frac{1}{2}\left(\frac{A}{\pi }\right)^{3/2}\left(e^{\delta ^2}H(\delta )\right),$$
(33)
with $`H(\delta )`$ and $`\delta `$ as defined in Eq. (28). The limiting cases are,

$$\underset{\delta 0}{lim}u_{\varphi \varphi }^{}=\underset{\delta 0}{lim}u_{rr}^{},$$
(34)
and

$$\underset{\delta \mathrm{}}{lim}u_{\varphi \varphi }^{}=\frac{1}{2}\underset{\delta \mathrm{}}{lim}u_{rr}^{}.$$
(35)
### C Stresses

We now proceed to calculate the radial stresses $`\sigma _{rr}^{}`$ which we measure in units of the shear modul $`G`$,

$$\sigma _{rr}^{}=\frac{2}{(12\nu )}\left((1\nu )u_{rr}^{}+2\nu u_{\varphi \varphi }^{}\right).$$
(36)
which is equivalent to

$$\frac{12\nu }{2}\left(\frac{A}{\pi }\right)^{3/2}\sigma _{rr}^{}=(1\nu )H(\delta )+\nu (e^{\delta ^2}H(\delta )).$$
(37)
We have plotted the rescaled stresses Eq. (37) in Fig. 4 for Poisson numbers $`\nu =0.2`$ and $`\nu =0.4`$.

The characteristic points of the rescaled stress $`\mathrm{\Sigma }_{rr}^{}`$ do now depend on $`\nu `$ albeit not very strongly, i.e., for lower $`\nu `$ one finds lower maximum tensile stresses, lower roots $`\delta _0`$ and higher compressive stresses. Generally spoken do the polar deformations $`u_{\varphi \varphi }^{}`$ strengthen the tensile stresses and weaken the compressive region. Higher contraction numbers $`\nu `$ do stronger weight polar deformations on radial stresses. The two limiting cases for the radial stresses are,

$$\underset{\delta 0}{lim}\sigma _{rr}^{}=\frac{2(1+\nu )}{3(12\nu )}\left(\frac{A}{\pi }\right)^{3/2},$$
(38)
and

$$\underset{\delta \mathrm{}}{lim}\sigma _{rr}^{}=\frac{1}{\pi (12\nu )r_{}^{}{}_{}{}^{3}}.$$
(39)
## IV Conclusion

We have considered a simplistic model for the fluid flow within a porous elastic solid. The presented description is given by linear equations. Essential simplifications are: Vanishing inertial forces of the solid phase and vanishing fluid compressibility. The first assumption relieves us of the problem of wave propagation in Biot media. The second assumption allows us to construct a linear description. We tried to obtain some general results for the linear problem. It appears that the general integrals for Biot’s and the thermo-elastic problem are different even within a linear description. While in both cases the body force field is governed by scalar potentials, e.g. pressure and temperature respectively, Biot’s problem is a fully cross-coupled field problem. This changes the structure of the governing Galerkin-Biot vector-potential, Eq. (15), substantially in such a way that mixed space-time derivatives do appear. However, we were able to show that an irrotational displacement field garantuees an equivalence of both problems. Though irrotational displacement fields are rare in geometries of practical interest, we studied the point symmetric case in three dimensions in order to exemplify the fluids action on the elastic strain/stress distribution. We quantitatively discussed the situation that arises when a pressure pulse within a ’borhole’ drives fluid through the solids pore-space. We found that radial tensile stresses arise in a region where non-permeable solids exhibit only compressive stresses. This surely is of interest to selected problem in soil meachanics and fracture mechanics, i.e. hydraulic fracturing.

## Acknowledgments

F. T. would like to acknowledge financial support from CEC under grant number ERBFMBICT950009 and thank for hospitality at NTNU Trondheim. Many thanks also go to Alex Hansen, Magnus Wangen and Thomas Ihle for stimulating discussions.