text:
# The Radiative Stress Tensor

## 1 Introduction

The presence of a radiation field will subject an otherwise perfect fluid to dissipation. In the limit of short photon mean free paths, dissipation is expressed as diffusion of momentum and energy and is regulated by the transport coefficients of viscosity and conductivity. These were obtained by Thomas (1930) who started from the relativistic transfer equation and developed the solutions in an expansion in photon mean free path. Since Thomas retained only terms of first order, his was a diffusion theory whose results have been extended to include the effects of general relativity (Hämeen-Antilla and Suhonen, 1971) Thomson scattering (Masaki, 1971; Hsieh and Spiegel, 1976) and Compton scattering (Masaki, 1981).

Unfortunately, a diffusion approximation does not provide an adequate description of the effects of radiation on the medium and its dynamics when vigorous radiative fluid dynamics produces large fluctuations in physical conditions. In such conditions, transparent regions may form and the photon mean free path can become relatively long so that better approximations are needed. One approach is to continue using the equations for the first two moments of the radiative intensity, energy density and flux, and closing the system by providing an improved approximation for the pressure tensor. A step in this direction is to use the Eddington factor methods (Mihalas and Mihalas, 1984) as has been done by a sequence of authors (Castor, 1972; Fu, 1987; Minerbo, 1983; Dominguez, 1995, 1996, 1997). As Castor (1972; see also (Mihalas, 1983)) has reported, the scheme seems to be usable for long photon mean free paths, but it does call for the introduction of additional information. Typically the variable Eddington factor is to be found from numerical solutions of the transfer equation in related conditions. It is therefore desirable to seek a self-contained approximation for the pressure tensor and we shall provide one here.

We shall here use an approach that originates, in part, with a method used for closure in the study of thermal effects of radiation (Unno and Spiegel, 1966). In the treatment of the thermal effects alone, the Eddington approximation gives rather good results. This is so because the Eddington closure relation for the radiative stress-energy tensor is the result of resummation of terms of all orders in an expansion in photon mean free path (Unno and Spiegel, 1966), and it captures aspects of both optically thin and thick structures. When we use the Eddington approximation we get a radiative heat equation (for a static medium) that is schematically of the form (Unno and Spiegel, 1966):

$$_tS=\kappa _1𝒟^2S,𝒟=\left(1\kappa _2^2\right)^1,$$
(1)
where $`\kappa _1`$ and $`\kappa _2`$ are suitable “diffusion” coefficients. The Laplacian, which is the usual diffusion operator, is here augmented by the additional smoothing operator $`𝒟`$. Of course, this is the bare essential of such a result, whose simplicity depends on a number of idealizations, most notably that the mean free path of photons is a constant, and on the neglect of retardation effects. Nevertheless, $`𝒟`$ carries much of the qualitative character of the radiative smoothing process. If we develop the smoothing operator in equation (1) for large scales (small $`^2`$, qualitatively speaking) we obtain $`𝒟^2^2+\kappa _2^4`$. This illustrates why we may recover the full expression for $`𝒟`$ from a Padé approximation based on only the first two terms in the expansion.

On the other hand, for the case of viscous stresses, the Eddington closure approximation fails to produce a physically acceptable result (Anderson and Spiegel, 1972) and an improved version is needed. The aim of the present paper is to devise a closure for the second order moments of the radiation field that may be used to deal with both the dynamical and thermal effects of the radiation field for both optically thick and thin structures. Similar techniques have recently been used in kinetic theory where an improved viscous stress tensor has been obtained by resumming terms up to only second order in an expansion in mean free path (Rosenau, 1989, 1993; Slemrod, 1997). The resummed forms, when expanded, contain terms of all orders, as in the earlier radiative thermal study.

Starting from the manifestly covariant form of the transfer equation for a moving medium, we develop the solution of the transfer equation in photon mean free path. Then we group terms within this expansion according to angular complexity. This leads to multiple sums and we rearrange the series so as to be able to compare terms of like character. From each order, we take the dominant term in power of the mean free path. Then we choose the leading portion at each order in terms of angular complexity, in what is really a Galerkin truncation. The resulting infinite sum is then converted to finite form by a resummation, or Padé, technique (Baker, 1975) to provide a compact form for the radiative viscosity tensor.

## 2 Basic Formulas

### 2.1 Transfer Equation

To include the effects of differential motion of a medium on the transfer of radiation, Thomas (1930) used the transfer equation in relativistic form. In this way, he readily included effects such as aberration and Doppler shift and was able to obtain the correct radiative viscosity where others had failed (Mihalas 1983). Even if one is not concerned with relativistic fluid velocities, it is convenient adopt the relativistic formulation and also to follow Hazelhurst and Sargent (1959) in using the covariant form of the equation of radiative transfer for the case of moving media. This approach is especially effective if one makes use of the one-particle distribution function for photons, $`f(p,x)`$, where $`p`$ (or $`p^\mu `$) is the four-momentum of a photon and $`x`$ (or $`x^\mu `$) is its location in space-time. According to the way one proceeds, $`f`$ is a scalar (as here) or a scalar density.

If we neglect gravitational bending of light and refraction by the medium, we may write the transfer equation in the general form (Anderson and Spiegel, 1972)

$$p^\mu _\mu f=\stackrel{~}{\rho }(\alpha \beta f)$$
(2)
where $`_\mu `$ stands for $`/x^\mu `$ and summation over any repeated index is understood. Here $`\alpha `$ is the rate of injection of photons with momentum $`p`$ by emission and $`\beta f`$ is the rate at which photons are removed by absorption; both rates are per unit mass. The quantity $`\stackrel{~}{\rho }`$ is the number density of absorbers times the rest mass per absorber.

Unlike $`f`$, the specific intensity, $``$, is not a scalar. Its transformation properties, as given by Thomas (1930), may be deduced from the definition

$$=\frac{h^4\nu ^3f}{c^3},$$
(3)
where $`h`$ is Planck’s constant, $`c`$ is the speed of light and $`\nu `$ is the frequency. We shall use units in which $`c`$ and $`h`$ are unity.

The transformation properties of the other quantities in the usual transfer equation are obtained by comparing that equation in the local rest frame of the matter with equation (2). For our purposes, there are two important reference frames, a basic inertial frame (such as that of a stationary star) and a frame locally comoving with the matter. The four-velocity of the matter, $`u^\mu `$, has components $`(1,0,0,0)`$ in the locally comoving frame.

The invariant

$$\stackrel{~}{\nu }=u^\mu p_\mu .$$
(4)
is the photon frequency in the local rest frame of the matter. By comparing the transfer equation in a general frame with that in the comoving frame, Thomas (1930) found $`\beta =\stackrel{~}{\nu }\kappa (\stackrel{~}{\nu })`$, where $`\kappa `$ is the opacity, and $`\alpha =\stackrel{~}{\nu }^2j(\stackrel{~}{\nu })`$, where $`j`$ is the emissivity. Since $`\alpha `$, $`\beta `$ and $`\stackrel{~}{\nu }`$ are scalars, we then know the transformation rules for $`j`$ and $`\kappa `$ once we introduce the (relativistic) Doppler relation between general frequency $`\nu `$ and $`\stackrel{~}{\nu }`$.

### 2.2 Radiative Stress Tensor

A basic object of radiative fluid dynamics is the radiative stress tensor, $`T^{\mu \nu }`$. It may be defined as an integral of the distribution function over momentum space:

$$T^{\mu \nu }=p^\mu p^\nu f𝑑P,$$
(5)
where $`dP`$ is the covariant volume element in momentum space. With $`\stackrel{~}{\nu }`$ as the radial coordinate in spherical polar coordinates in momentum space the element of volume in momentum space is $`dP=\stackrel{~}{\nu }d\stackrel{~}{\nu }d\mathrm{\Omega }`$, where $`d\mathrm{\Omega }`$ is the element of solid angle (Landau and Lifshitz, 1984). We may then factor the frequency from $`p^\mu `$ and define the null vector

$$n^\mu =p^\mu /\stackrel{~}{\nu },$$
(6)
so that the stress tensor is

$$T^{\mu \nu }=𝑑\mathrm{\Omega }n^\mu n^\nu I$$
(7)
where we have introduced the frequency-integrated specific intensity

$$I=_0^{\mathrm{}}\stackrel{~}{\nu }^3f𝑑\stackrel{~}{\nu }.$$
(8)
The projection operator

$$h^{\mu \nu }=\eta ^{\mu \nu }u^\mu u^\nu $$
(9)
picks out space-like components in the frame of the matter, where $`\eta ^{\mu \nu }`$ is the Minkowski metric, diag$`(1,1,1,1)`$. Then we may introduce the spacelike unit four-vector,

$$l^\mu =h_\nu ^\mu n^\nu ,$$
(10)
so that

$$n^\mu =u^\mu +l^\mu .$$
(11)
We see that $`l^\mu `$ must be perpendicular to $`u^\mu `$ and that $`l_\mu l^\mu =1`$. Then with the introduction of (11) we arrive at a standard way to write the stress tensor, namely

$$T^{\mu \nu }=P^{\mu \nu }+F^\mu u^\nu +u^\mu F^\nu +Eu^\mu u^\nu $$
(12)
with

$$E=I𝑑\mathrm{\Omega },F^\mu =l^\mu I𝑑\mathrm{\Omega },P^{\mu \nu }=l^\mu l^\nu I𝑑\mathrm{\Omega },$$
(13)
where these integrals are over all solid angle.

By taking moments of the transfer equation itself, we derive equations for moments like those in (12). This procedure leads to an infinite hierarchy of moment equations that may be truncated by introducing a relation amongst the moments, $`E`$, $`F^\mu `$ and $`P^{\mu \nu }`$. Such closure relations can be found by assuming some particular forms for the angular dependence of $`I`$, usually expressed in terms of the $`l^\mu `$ or $`n^\mu `$ and this requires angular integrals of their products. Though we do not follow this procedure, we do need such integrals and so we give a list of formulas for the angular integrals in Appendix A. In the simplest case, when the isotropic approximation is introduced for the intensity, we need the integral

$$l^\mu l^\nu 𝑑\mathrm{\Omega }=\frac{4\pi }{3}h^{\mu \nu }$$
(14)
for the evaluation of $`P^{\mu \nu }`$ that is used in the Eddington closure approximation (Anderson and Spiegel, 1972).

### 2.3 The Matter and the Total Stress Tensors

To describe the matter macroscopically, we adopt, with Weinberg (1971), the procedures given by Eckart (1940), who noted that a vector $`V^\mu `$ may be written as $`V^\mu =(u_\alpha V^\alpha )u^\mu +h_\alpha ^\mu V^\alpha `$, as in equation (11). A similar decomposition is possible for tensors and, in particular, for the radiative stress tensor, we see that

$$E=u_\alpha u_\beta T^{\alpha \beta },F^\mu =u_\alpha T^{\alpha \beta }h_\beta ^\mu ,P^{\mu \nu }=h_\alpha ^\mu h_\beta ^\nu T^{\alpha \beta }.$$
(15)
As in Eckart (1940)’s approach, we may use the form of (15) to define the three important moments of kinetic theory in terms of the matter stress tensor. If we denote the stress-energy tensor of the matter as $`𝒯^{\mu \nu }`$, the mass-energy density, the matter flux and the pressure tensor are given (as in (15)) by

$$=u_\alpha u_\beta 𝒯^{\alpha \beta },^\mu =u_\alpha 𝒯^{\alpha \beta }h_\beta ^\mu ,𝒫^{\mu \nu }=h_\alpha ^\mu h_\beta ^\nu 𝒯^{\alpha \beta }.$$
(16)
To concentrate on radiative effects, we assume that, in the absence of interaction with radiation, the matter would be a perfect fluid, for which $`^\mu =0`$ and $`𝒫^{\mu \nu }=𝒫h^{\mu \nu }`$, where $`𝒫`$ is the partial pressure of the matter. The matter stress tensor in that case would be

$$𝒯^{\mu \nu }=u^\mu u^\nu 𝒫h^{\mu \nu }.$$
(17)
The total stress tensor (matter plus radiation) is

$$\mathrm{\Theta }^{\mu \nu }=T^{\mu \nu }+𝒯^{\mu \nu }$$
(18)
and the equations of motion of the mixture are expressed as

$$\mathrm{\Theta }_{,\nu }^{\mu \nu }=0.$$
(19)
## 3 Diffusion Theory

For a discussion of the extension of Thomas’ theory, a brief restatement of his results is useful for bringing out some of the issues to be faced. To provide one here, we introduce the source function, $`𝒮=\alpha /\beta `$, and the local photon mean free path, $`\epsilon =(\stackrel{~}{\rho }\kappa )^1`$ where $`\kappa `$ is generally a function of $`\stackrel{~}{\nu }`$ as well as of the state variables of the medium. The mean free path appears explicitly in the transfer equation (2), which we rewrite as

$$f=𝒮\epsilon n^\mu f_{,\mu }$$
(20)
where $`n^\mu `$ is defined in (6). This way of writing the transfer equation is meant to take advantage of the smallness of $`\epsilon `$, which is assumed in Thomas’ theory. The fact that $`\epsilon `$ appears in front of the derivative means that the perturbation theory based on its smallness will be singular. Generally speaking, singular perturbations expansions give rise to (divergent) asymptotic series but these may be made useful by appropriate (re)summation techniques (Bender, 1978) such as we shall use here.

In first approximation, we have $`f𝒮`$. If we iterate on this once in (20) we get Thomas’ approximation,

$$f=𝒮\epsilon n^\mu 𝒮_{,\mu }.$$
(21)
To compute the stress tensor, we first perform the integration over frequency after multiplying by $`\stackrel{~}{\nu }^3`$ to obtain an approximation for the frequency-integrated intensity,

$$I=Sn^\mu _0^{\mathrm{}}\stackrel{~}{\nu }^3\epsilon 𝒮_{,\mu }𝑑\stackrel{~}{\nu },$$
(22)
where we have introduced the integrated source function

$$S(T)=_0^{\mathrm{}}\stackrel{~}{\nu }^3𝒮(\stackrel{~}{\nu },T)𝑑\stackrel{~}{\nu }.$$
(23)
The integral in formula (22) is a simple example of the kind of integrals that the higher order development entails, so it is worth discussing its evaluation.

In Thomas theory, $`\stackrel{~}{\nu }^3𝒮`$ is the Planck distribution and it depends explicitly on $`\stackrel{~}{\nu }`$ and temperature, $`T`$, from which it derives an implicit dependence on position and time. Therefore we may write

$$𝒮_{,\mu }=\stackrel{~}{\nu }_{,\mu }_{\stackrel{~}{\nu }}𝒮+T_{,\mu }_T𝒮.$$
(24)
On recalling that $`\stackrel{~}{\nu }=u_\mu p^\mu `$ and that $`p^\mu `$ does not depend on coordinates, we obtain $`\stackrel{~}{\nu }_{,\mu }=p_\nu u_{,\mu }^\nu `$. Then, since $`u_\nu u_{,\mu }^\nu =0`$, we have that

$$\stackrel{~}{\nu }_{,\mu }=\stackrel{~}{\nu }n_\nu u_{,\mu }^\nu .$$
(25)
If, as in local thermodynamic equilibrium, the source function depends on $`\stackrel{~}{\nu }`$ and $`T`$ only through the combination $`\stackrel{~}{\nu }/T`$, we have

$$\stackrel{~}{\nu }_{\stackrel{~}{\nu }}=T_T$$
(26)
and so (24) can be rewritten as

$$𝒮_{,\mu }=\left(T_{,\mu }Tn^\rho u_{\rho ,\mu }\right)_T𝒮.$$
(27)
The integral in (22) then becomes

$$_0^{\mathrm{}}\stackrel{~}{\nu }^3\epsilon 𝒮_{,\mu }𝑑\stackrel{~}{\nu }=(T_{,\mu }Tn^\nu u_{\nu ,\mu })_0^{\mathrm{}}\stackrel{~}{\nu }^3\epsilon _T𝒮d\stackrel{~}{\nu }.$$
(28)
If we pull the $`T`$ derivative out of the integral, we get an extra term involving the $`T`$ derivative of $`\epsilon `$. On the other hand, we can simply introduce a weighted frequency integral through the relation

$$\overline{\epsilon }S_{,T}=_0^{\mathrm{}}\stackrel{~}{\nu }^3𝒮_{,T}(\stackrel{~}{\nu },T)\epsilon 𝑑\stackrel{~}{\nu }.$$
(29)
The quantity $`\overline{\epsilon }`$ is a mean of the inverse of the absorption coefficient with respect to $`𝒮_{,T}`$ in the spirit of the Rosseland mean absorption coefficient. Then we find that, for $`TS_{,T}=4S`$,

$$I=S\overline{\epsilon }n^\mu S_{,\mu }+4\overline{\epsilon }n^\mu n^\nu Su_{\mu ,\nu }.$$
(30)
In the higher order theory to be discussed below, we can use similar tricks by introducing a variety of mean absorption coefficients, as is done sometimes in the theory of stellar atmospheres. More often than not, all these are assumed to be equal in the final analysis. To avoid such purely formal maneuvers, we shall assume in the presentation of higher order theory that the medium is gray.

Then, to this approximation, partial integrations such as we have just described give us the known forms of moments of the intensity in the Thomas approximation:

$$E=4\pi S4\pi \overline{\epsilon }\left[u^\alpha S_{,\alpha }+\frac{4}{3}Su_{,\alpha }^\alpha \right],$$
(31)
$$F^\mu =\frac{4\pi }{3}\overline{\epsilon }h^{\mu \alpha }\left[S_{,\alpha }4Su^\beta u_{\alpha ,\beta }\right],$$
(32)
$`P^{\mu \nu }={\displaystyle \frac{1}{3}}h^{\mu \nu }E+{\displaystyle \frac{16\pi }{15}}\overline{\epsilon }Su_{\rho ,\sigma }\tau ^{\mu \nu \rho \sigma },`$ (33)
where we have used the relevant formulae from Appendix A and have introduced

$$\tau ^{\mu \nu \rho \sigma }=h^{\mu \rho }h^{\nu \sigma }+h^{\mu \sigma }h^{\nu \rho }\frac{2}{3}h^{\mu \nu }h^{\rho \sigma }.$$
(34)
From the expression for the flux we see that, since $`S`$ is proportional to $`T^4`$, there is the usual conduction term proportional to the (four) gradient of $`T`$ with a conductivity coefficient $`16\pi \overline{\epsilon }S/(3T)`$. This is the diffusion limit and it is much less good than the Eddington approximation, for which the flux is proportional to the gradient of $`E`$. With the Eddington approximation, thermal times are well approximated over all ranges of optical thickness (Unno and Spiegel, 1966). The second term in (32) has been considered to be a purely relativistic term (Mihalas, 1983), but a classical analogue exists in the kinetic theory of rarefied gases (Chen, PhD thesis; Chen et al., 2000).

In the expression for the pressure tensor, the first term, to the order given, does mimic the Eddington approximation and we have in addition a viscous term. From the standard form of the pressure tensor, we can read off the viscosity coefficient given by Thomas, $`16ϵ\pi S/15`$. Like the Thomas expression for the flux, this is also a diffusive approximation good only for short mean free paths. Our aim here is to improve on the non-diagonal viscous tensor in the same way that the Eddington approximation improves on the conduction term. In a subsequent paper, we shall describe an improved approximation for the diagonal terms for the case where scattering is important.

## 4 The Expansion of $`f`$

To go beyond the Thomas approximation, we could continue the iteration process described in the previous section, but we find it more effective to proceed by formal expansions. Our small parameter will be a typical value of $`\epsilon `$, call it $`ϵ`$. Then we form the expansion

$$f=\underset{m=0}{\overset{\mathrm{}}{}}f_{(m)}ϵ^m,$$
(35)
which we introduce into the transfer equation. On demanding that the expanded equation is satisfied term by term and prescribing that $`\epsilon `$ is everywhere of order $`ϵ`$, we get

$$f_{(0)}=𝒮,$$
(36)
as in comoving local thermodynamic equilibrium and, for the higher orders, we obtain

$$f_{(m+1)}=\frac{\epsilon }{ϵ}n^\mu f_{(m),\mu }.$$
(37)
In the first order, as in diffusion theory, the complication caused by having $`\epsilon `$ in the formulas is slight, but in higher orders, it makes the development very cumbersome especially on account of the presence of derivatives of $`\epsilon `$ in the expressions of $`f_{(m)}`$. These may be converted to derivatives with respect to rest frequency and state variables as for the source function but still, in the performance of the integrals over frequency, it soon becomes necessary to introduce a plethora of mean absorption coefficients, and even tensorial mean absorption coefficients (Anderson and Spiegel, 1972). Though we do not report on this here, we have retained variation of $`\epsilon `$ in calcualtions out to the third order in $`\epsilon `$ in (Chen, PhD thesis). This is far enough for one aspect of our approach to be carried out, though we shall not go into this here.

More generally, it is possible to make allowances for the dependence of $`\epsilon `$ on $`x^\mu `$ by a formal device like used in stellar atmosphere theory. We may introduce new coordinates $`y^\mu `$ through the transformation $`ϵdy^\mu =\overline{\epsilon }dx^\mu `$ where $`\overline{\epsilon }`$ is a mean (over frequency) absorption coefficient. This formal simplification is analogous to a transformation to the local rest frame of the matter. But since both $`\overline{\epsilon }`$ and the $`u^\mu `$ are dependent on the radiation field, such implicit definitions of coordinates necessitate, at the end of the calculations, some rather hard inversions. Moreover, we would need to take frequency means to use such transformations effectively. We shall avoid all these complications here and will stick to inertial coordinates, requiring that $`\epsilon `$ is constant. In that case, we obtain

$$f_{(m+1)}=n^\mu f_{(m),\mu },$$
(38)
or, in terms of $`𝒮`$,

$$f_{(m)}=(1)^m\frac{D^m𝒮}{Ds^m},$$
(39)
where

$$\frac{D}{Ds}=n^\mu _\mu .$$
(40)
To perform the indicated differentiations, it is useful to have a formula for the derivative of the null vector, $`n^\mu `$. On differentiating $`p^\mu =n^\mu \stackrel{~}{\nu }`$ and making use of (25), we find that

$$n_{,\nu }^\mu =n^\mu n^\rho u_{\rho ,\nu }.$$
(41)
We can then develop explicit formulas for the $`f_{(m)}`$ in terms of the source function and these together lead to a series providing $`f`$ as a functional of $`𝒮`$. This is tantamount to what Chandrasekhar (1960) has called the formal solution of the transfer equation. (Its integral expression in terms of optical distance is helpful in treating the case of variable $`\epsilon `$.) The idea of the formal solution is that the source function and the absorption coefficient are material properties of the medium and are to be considered as given for the purposes of solving the transfer problem, though the situation is reversed when we try to solve equations for the properties of the medium. So we need to decide how the source function is to be specified.

In the simplest case mathematically, we would be given $`𝒮`$ explicitly as a function of space and time. From this, we would compute frequency integrals to get the integrated intensity through formula (8). But it is more representative of the general situation to assume, as we did for diffusion theory, that $`𝒮`$ is expressible in terms of $`T`$ and $`\stackrel{~}{\nu }`$. Though the source function does not depend on direction, we do have effects of angular distribution entering in through the transfer equation itself. To make allowances for that also, we write

$$_\mu =T_{,\mu }_T+p^\rho u_{\rho ,\mu }_{\stackrel{~}{\nu }}+n_{,\mu }^\sigma _{n^\sigma }$$
(42)
and we are led to a generalization of the differentiation formula implicit in (30):

$$_\mu =\left(T_{,\mu }Tn^\rho u_{\rho ,\mu }\right)_Tn^\sigma n^\rho u_{\rho ,\mu }_{n^\sigma }.$$
(43)
On noting that the quantity $`f_{(0)}=𝒮`$ does not depend on direction ($`n^\mu `$) we can evaluate $`f_{(1)}`$. We find using (43) that

$$f_{(1)}=𝒮_{,T}T_{,\mu }n^\mu +T𝒮_{,T}u_{\nu ,\mu }n^\mu n^\nu .$$
(44)
If we use this formula to compute the moments, we recover the results of the previous section with $`ϵ`$ in place of $`\overline{\epsilon }`$.

From (39), we see on using (41), that

$$f_{(2)}=n^\mu n^\nu 𝒮_{,\mu \nu }n^\mu n^\nu n^\rho u_{\rho ,\mu }𝒮_{,\nu }.$$
(45)
We then have to cope with the higher derivatives of the source function $`𝒮`$. If we continue to assume that $`𝒮`$ depends only on $`\stackrel{~}{\nu }/T`$, since $`𝒮`$ does not depend on $`n^\mu `$, (43) gives us

$$𝒮_{,\mu }=T𝒮_{,T}[(lnT)_{,\mu }(ln\stackrel{~}{\nu })_{,\mu }].$$
(46)
A similar, but lengthier, calculation leads to

$`𝒮_{,\mu \nu }`$ $`=`$ $`T𝒮_{,T}\{[(lnT)_{,\mu }(ln\stackrel{~}{\nu })_{,\mu }](lnT)_{,\nu }+(lnT)_{,\mu \nu }(ln\stackrel{~}{\nu })_{,\mu \nu }\}`$
$`+`$ $`[T^2𝒮_{,TT}(lnT)_{,\nu }(T𝒮_{,T}+T^2𝒮_{,TT})(ln\stackrel{~}{\nu })_{,\nu }]\left[(lnT)_{,\mu }(ln\stackrel{~}{\nu })_{,\mu }\right]`$
where we have used (25) and (41).

We then find that (45) takes the form

$$f_{(2)}=n^\mu n^\nu _{\mu \nu }^{(22)}n^\mu n^\nu n^\rho _{\mu \nu \rho }^{(23)}+n^\mu n^\nu n^\rho n^\sigma _{\mu \nu \rho \sigma }^{(24)}$$
(48)
with

$`_{\mu \nu }^{(22)}`$ $`=`$ $`𝒮_{,T}T_{,\mu \nu }+𝒮_{,TT}T_{,\mu }T_{,\nu }`$
$`_{\mu \nu \rho }^{(23)}`$ $`=`$ $`T𝒮_{,T}u_{\mu ,\nu \rho }+(3𝒮_{,T}+2T𝒮_{,TT})u_{\mu ,\nu }T_{,\rho }`$
$`_{\mu \nu \rho \sigma }^{(24)}`$ $`=`$ $`(3T𝒮_{,T}+T^2𝒮_{,TT})u_{\mu ,\nu }u_{\rho ,\sigma }.`$ (49)
This way of writing $`f_{(2)}`$ is useful in computing terms in the next order and from (39) and (40) and the forms of $`f_{(1)}`$ and $`f_{(2)}`$, we may surmise that

$$f_{(m)}=\underset{k=m}{\overset{2m}{}}(1)^k_{\mu _1\mathrm{}\mu _k}^{(mk)}n^{\mu _1}\mathrm{}n^{\mu _k}$$
(50)
where the quantities $``$ are functions of $`𝒮`$ and its derivatives with respect to $`T`$; that is, they are functionals of $`𝒮`$. (If we had not assumed constant photon mean free path, they would be functionals of $`\mathrm{log}\epsilon `$ as well.) The calculation of the distribution function at any order then boils down to the calculation of these functionals, which is a straightforward, if demanding, task. For $`m=3`$, we are led to the results

$$_{\mu \nu \alpha }^{(33)}=𝒮_{,TTT}T_{,\mu }T_{,\nu }T_{,\alpha }+3𝒮_{,TT}T_{,\mu \nu }T_{,\alpha }+𝒮_{,T}T_{,\mu \alpha \nu }$$
(51)
$`_{\mu \nu \alpha \beta }^{(34)}`$ $`=`$ $`T𝒮_{,T}u_{\alpha ,\beta \mu \nu }+(4𝒮_{,T}+3T𝒮_{,TT})u_{\alpha ,\beta \mu }T_{,\nu }+(6𝒮_{,T}+3T𝒮_{,TT})u_{\alpha ,\beta }T_{,\mu \nu }`$ (52)
$`+(8𝒮_{,TT}+3T𝒮_{,TTT})u_{\alpha ,\beta }T_{,\mu }T_{,\nu }`$
$`_{\mu \nu \alpha \beta \rho }^{(35)}`$ $`=`$ $`u_{\alpha ,\beta }[(10T𝒮_{,T}+3T^2𝒮_{,TT})u_{\mu ,\nu \rho }`$ (53)
$`+(15𝒮_{,T}+18T𝒮_{,TT}+3T^2𝒮_{,TTT})u_{\mu ,\nu }T_{,\rho }]`$
$$_{\mu \nu \rho \sigma \alpha \beta }^{(36)}=[15T𝒮_{,T}+9T^2𝒮_{,TT}+T^3𝒮_{,TTT}]u_{\mu ,\nu }u_{\rho ,\sigma }u_{\alpha ,\beta }$$
(54)
We thus have all the terms up to third order.

In this way we see that we may order the terms by increasing angular complexity. In that case, if we go up to four factors of $`n^\mu `$, we have an expansion of $`f`$ in the form

$`f`$ $`=`$ $`𝒮ϵ_\mu ^{(11)}n^\mu +ϵ[_{\mu \nu }^{(12)}+ϵ_{\mu \nu }^{(22)}]n^\mu n^\nu ϵ^2[_{\mu \nu \rho }^{(23)}+ϵ_{\mu \nu \rho }^{(33)}]n^\mu n^\nu n^\rho `$ (55)
$`ϵ^2[_{\mu \nu \rho \sigma }^{(24)}+ϵ_{\mu \nu \rho \sigma }^{(34)}+ϵ^2_{\mu \nu \rho \sigma }^{(44)}]n^\mu n^\nu n^\rho n^\sigma +\mathrm{}.`$
We may compute the coefficients $``$ at any order using this procedure.

## 5 The Intensity Expansion

An intermediate step on the way to computing the moments that make up the radiative stress tensor, is the computation of the (frequency) integrated intensity defined in (8). Its expansion is

$$I=\underset{m=0}{\overset{\mathrm{}}{}}ϵ^mI_{(m)},$$
(56)
and when we introduce (50) into (8) we obtain for the coefficients in (56),

$$I_{(m)}=\underset{k=m}{\overset{2m}{}}(1)^kY_{\mu _1\mathrm{}\mu _k}^{(m,k)}n^{\mu _1}\mathrm{}n^{\mu _k}$$
(57)
where

$$Y_{\mu _1\mathrm{}\mu _k}^{(m,k)}=\stackrel{~}{\nu }^3_{\mu _1\mathrm{}\mu _k}^{(mk)}𝑑\stackrel{~}{\nu }.$$
(58)
With the expressions for the $`_{\mu _1\mathrm{}\mu _k}^{(mk)}`$ given in the previous section, we may develop explicit expressions for $`Y_{\mu _1\mathrm{}\mu _k}^{(m,k)}`$. We do this by working out a recursion formula for these coefficients. First, we introduce the expression (50) for $`f_{(m)}`$ into (38) to obtain

$$f_{(m+1)}=n^\mu \underset{k=m}{\overset{2m}{}}(1)^k\left[_{\mu _1\mathrm{}\mu _k,\mu }^{(mk)}n^{\mu _1}\mathrm{}n^{\mu _k}+_{\mu _1\mathrm{}\mu _k}^{(mk)}(n^{\mu _1}\mathrm{}n^{\mu _k})_{,\mu }\right].$$
(59)
Since from (41), we have

$$(n^{\mu _1}\mathrm{}n^{\mu _k})_{,\mu }=kn^\rho u_{\rho ,\mu }(n^{\mu _1}\mathrm{}n^{\mu _k}),$$
(60)
the only real complication comes from the derivatives of the $``$.

Next, we reexamine the derivation of the explicit formulae for the $``$, and see that their dependence on position and time comes about through their dependence on $`\stackrel{~}{\nu }`$, on $`T`$ and on $`u^\mu `$, as well as the derivatives of these quantities with respect to $`x^\mu `$. So each $``$ depends on $`\stackrel{~}{\nu }`$ and on $`N`$ other variables that we shall call $`V^A`$, where $`A=1,2,\mathrm{},N`$. The $`V^A`$ are $`T_{,\mu }`$, $`T_{,\mu \nu }`$, …, $`u_{\mu ,\nu }`$, and so on. Let

$$\mathrm{\Lambda }_\mu ^{(m,k)}=V_{,\mu }^A_{V^A}$$
(61)
where summation over repeated $`A`$ is understood. Then

$$_\mu =\stackrel{~}{\nu }_{,\mu }_{\stackrel{~}{\nu }}+\mathrm{\Lambda }_\mu ^{(m,k)}$$
(62)
On recalling that $`\stackrel{~}{\nu }_{,\mu }=\stackrel{~}{\nu }n^\rho u_{\rho ,\mu }`$ we can rewrite (59) as

$`f_{(m+1)}`$ $`=`$ $`{\displaystyle \underset{k=m}{\overset{2m}{}}}(1)^{k+1}\stackrel{~}{\nu }n^\rho n^\mu u_{\rho ,\mu }_{\stackrel{~}{\nu }}_{\mu _1\mathrm{}\mu _k}^{(mk)}n^{\mu _1}\mathrm{}n^{\mu _k}`$ (63)
$`+`$ $`{\displaystyle \underset{k=m}{\overset{2m}{}}}(1)^{k+1}\left[n^\mu \mathrm{\Lambda }_\mu ^{(mk)}kn^\rho n^\mu u_{\rho ,\mu }\right]_{\mu _1\mathrm{}\mu _k}^{(mk)}n^{\mu _1}\mathrm{}n^{\mu _k}`$
To compute the $`I_{(m)}`$, we carry out the frequency integral indicated in (8) and integrate by parts in the term with $`_{\stackrel{~}{\nu }}`$ to obtain

$$I_{(m+1)}=\underset{k=m}{\overset{2m}{}}(1)^k\left[(k+4)n^\rho n^\mu u_{\rho ,\mu }n^\mu \mathrm{\Lambda }_\mu ^{(m,k)}\right]Y_{\mu _1\mathrm{}\mu _k}^{(m,k)}n^{\mu _1}\mathrm{}n^{\mu _k}.$$
(64)
Now $`\stackrel{~}{\nu }`$ is no longer in the problem, and so (62) simplifies to $`\mathrm{\Lambda }_\mu ^{(m,k)}=_\mu `$; we shall make this replacement in the rest of the development.

Formula (64) has two terms, each having a different angular structure. To keep to our general attempt to order things by angular structure as well as by mean free path, we make a further rearrangement. So far, in the expressions for the $``$, for different pairs $`(mk)`$, we have had need only of values for which $`mk2m`$. Therefore we are at liberty to impose that

$$_{\mu _1\mathrm{}\mu _k}^{(mk)}=0\mathrm{for}k<m\mathrm{and}k>2m.$$
(65)
Accordingly

$$Y_{\mu _1\mathrm{}\mu _k}^{(m,k)}=0\mathrm{for}k<m\mathrm{and}k>2m.$$
(66)
Moreover, for either $`m<0`$ or $`k<0`$, we define $`Y_{\mu _1\mathrm{}\mu _k}^{(m,k)}`$ to be zero.

Now we can adjust the indices in the summation (64) by replacing $`k`$ in the first term with $`k2`$ and $`k`$ in the second term with $`k1`$. We then have

$$I_{(m+1)}=\underset{k=m+1}{\overset{2(m+1)}{}}(1)^k\left[(k+2)Y_{\mu _1\mathrm{}\mu _{k2}}^{(m,k2)}u_{\mu _{k1},\mu _k}+Y_{\mu _1\mathrm{}\mu _{k1},\mu _k}^{(m,k1)}\right]n^{\mu _1}\mathrm{}n^{\mu _k}$$
(67)
and, on comparing with (57), we observe that

$$Y_{\mu _1\mathrm{}\mu _k}^{(m+1,k)}=Y_{\mu _1\mathrm{}\mu _{k1},\mu _k}^{(m,k1)}+(k+2)Y_{\mu _1\mathrm{}\mu _{k2}}^{(m,k2)}u_{\mu _{k1},\mu _k},$$
(68)
which will be used in the sequel. The $`Y`$ sequence is to be built up from $`Y^{(0,0)}`$. As we see from (58), this is the integrated source function, that is

$$Y^{(0,0)}=S.$$
(69)
Hence, as seen above, $`I`$ is a functional of $`S`$ as well as of $`u^\mu `$ and the angular variables.

## 6 The Viscosity Tensor

Once we have calculated the integrated intensity, as in the previous section, we can evaluate the physically important angular moments,

$$E=\underset{m=0}{\overset{\mathrm{}}{}}ϵ^mE_{(m)},F^\mu =\underset{m=0}{\overset{\mathrm{}}{}}ϵ^mF_{(m)}^\mu ,P^{\mu \nu }=\underset{m=0}{\overset{\mathrm{}}{}}ϵ^mP_{(m)}^{\mu \nu }.$$
(70)
Since the $`Y`$ defined in (58) are independent of angle, the terms in these expansions are

$$E_{(m)}=\underset{k=m}{\overset{2m}{}}(1)^kY_{\mu _1\mathrm{}\mu _k}^{(m,k)}n^{\mu _1}n^{\mu _2}\mathrm{}n^{\mu _k}𝑑\mathrm{\Omega }$$
(71)
$$F_{(m)}^\mu =\underset{k=m}{\overset{2m}{}}(1)^kY_{\mu _1\mathrm{}\mu _k}^{(m,k)}l^\mu n^{\mu _1}n^{\mu _2}\mathrm{}n^{\mu _k}𝑑\mathrm{\Omega }$$
(72)
$$P_{(m)}^{\mu \nu }=\underset{k=m}{\overset{2m}{}}(1)^kY_{\mu _1\mathrm{}\mu _k}^{(m,k)}l^\mu l^\nu n^{\mu _1}n^{\mu _2}\mathrm{}n^{\mu _k}𝑑\mathrm{\Omega }.$$
(73)
Though we shall quote some results about the approximations to $`E`$ and $`F^\mu `$, our main interest is in providing an improved treatment of the viscous effects for transparent regions of the system under study. So we shall detail only the calculation of an improved approximation for the pressure tensor. If we introduce the Eddington closure for the diagonal terms in the stress tensor, we may write

$$P^{\mu \nu }=\frac{1}{3}Eh^{\mu \nu }+\mathrm{\Xi }^{\mu \nu },$$
(74)
where $`\mathrm{\Xi }^{\mu \nu }`$ is the viscous stress tensor, whose expression the main goal of this work. As mentioned, in a later paper, we shall show how the first term in (74) can be improved considerably.

To seek a representation for $`\mathrm{\Xi }^{\mu \nu }`$ that is suitable for both long and short photon mean free paths, we begin by expanding the viscosity tensor as

$$\mathrm{\Xi }^{\mu \nu }=\underset{m=0}{\overset{\mathrm{}}{}}ϵ^m\mathrm{\Xi }_{(m)}^{\mu \nu }.$$
(75)
From (71) and (73) we find the coefficients in this expansion and so can write

$$\mathrm{\Xi }_{(m)}^{\mu \nu }=\underset{k=m}{\overset{2m}{}}(1)^kY_{\mu _1\mathrm{}\mu _k}^{(m,k)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _k},$$
(76)
where

$$\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _k}(l^\mu l^\nu +\frac{1}{3}h^{\mu \nu })n^{\mu _1}n^{\mu _2}\mathrm{}n^{\mu _k}𝑑\mathrm{\Omega }.$$
(77)
Since the trace of $`l^\mu l^\nu +\frac{1}{3}h^{\mu \nu }`$ is zero, $`\mathrm{\Xi }_{(m)}^{\mu \nu }`$ and $`\mathrm{\Xi }^{\mu \nu }`$, are traceless, symmetric tensors.

For the $`Y`$s we have the recursion formula (68) given at the end of the previous section, while to reexpress the angular integrals usefully, we can refer to the formula of Appendix A. We see from these that

$$\mathrm{\Omega }^{\mu \nu }=0,\mathrm{\Omega }^{\mu \nu \rho }=0\mathrm{and}\mathrm{\Omega }^{\mu \nu \rho \sigma }=\frac{4\pi }{15}\tau ^{\mu \nu \rho \sigma },$$
(78)
where $`\tau ^{\mu \nu \rho \sigma }`$ is defined in (34). To extend the results of Thomas we next generalize $`\tau ^{\mu \nu \rho \sigma }`$ into a series of orthogonal polynomials that represent the expansion of $`\mathrm{\Omega }`$.

If we use (A6), we can develop the products of the $`n^\mu `$ and express $`\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _k}`$ in terms of the moments of $`l^\mu `$. Then we have

$$\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _k}=\underset{i=0}{\overset{k}{}}\left(\begin{array}{c}i\\ k\end{array}\right)W^{\mu \nu (\mu _1\mathrm{}\mu _i}u^{\mu _{i+1}}\mathrm{}u^{\mu _k)}$$
(79)
where

$$W^{\mu \nu \mu _1\mathrm{}\mu _{\mathrm{}}}=M^{\mu \nu \mu _1\mathrm{}\mu _{\mathrm{}}}+\frac{1}{3}h^{\mu \nu }M^{\mu _1\mathrm{}\mu _{\mathrm{}}}$$
(80)
and

$$M^{\mu _1\mathrm{}\mu _{\mathrm{}}}=l^{\mu _1}l^{\mu _2}\mathrm{}l^\mu _{\mathrm{}}𝑑\mathrm{\Omega }.$$
(81)
The integral for an odd number of factors $`l^\mu `$ is zero so, for odd $`k`$, $`M^{\mu _1\mathrm{}\mu _k}0`$, and hence $`W^{\mu \nu \mu _1\mathrm{}\mu _k}0`$. Therefore we may rewrite (79) in a way that will be convenient later, namely,

$$\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _k}=\underset{i=0}{\overset{[\frac{k}{2}]}{}}\left(\begin{array}{c}2i\\ k\end{array}\right)W^{\mu \nu (\mu _1\mathrm{}\mu _{2i}}u^{\mu _{2i+1}}\mathrm{}u^{\mu _k)}$$
(82)
where $`[\frac{k}{2}]`$ means the integer part of $`\frac{k}{2}`$.

When $`k=1`$, the upper limit of the sum is $`i=0`$, and $`\mathrm{\Omega }^{\mu \nu \rho }=W^{\mu \nu }u^\rho `$. Because $`M=4\pi `$ and, according to (14), $`M^{\mu \nu }=\frac{4\pi }{3}h^{\mu \nu }`$, we have $`W^{\mu \nu }=0`$ and $`\mathrm{\Omega }^{\mu \nu \rho }=0`$ as just noted. In the same way, at the next level, from (A2) we see that

$$\tau ^{\mu \nu \rho \sigma }=W^{\mu \nu \rho \sigma }$$
(83)
which is equivalent to the third expression in (78). We then define a generalized $`\tau `$ as

$$\tau ^{\mu \nu \mu _1\mathrm{}\mu _2\mathrm{}}=W^{\mu \nu \mu _1\mathrm{}\mu _2\mathrm{}}+\frac{1}{3}W^{\mu \nu \mu _1\mathrm{}\mu _{2\mathrm{}3}\mu _{2\mathrm{}2}}h^{\mu _{2\mathrm{}1}\mu _2\mathrm{}}.$$
(84)
To introduce this into (82) we rewrite (84) as

$$W^{\mu \nu \mu _1\mathrm{}\mu _2\mathrm{}}=\tau ^{\mu \nu \mu _1\mathrm{}\mu _2\mathrm{}}\frac{1}{3}W^{\mu \nu \mu _1\mathrm{}\mu _{2\mathrm{}3}\mu _{2\mathrm{}2}}h^{\mu _{2\mathrm{}1}\mu _2\mathrm{}},$$
(85)
where it is understood that $`h^{\mu _p\mu _q}`$ is zero if either $`p`$ or $`q`$ is negative. Hence, for $`\mathrm{}=1`$ we recover (83) and we may then proceed recursively to express $`W`$ in terms of $`\tau `$. When $`\mathrm{}=2`$ we have

$$W^{\mu \nu \mu _1\mu _2\mu _3\mu _4}=\tau ^{\mu \nu \mu _1\mu _2\mu _3\mu _4}\frac{1}{3}\tau ^{\mu \nu \mu _1\mu _2}h^{\mu _3\mu _4}.$$
(86)
From this we see that

$$W^{\mu \nu \mu _1\mathrm{}\mu _{2p}}=\underset{q=0}{\overset{p}{}}\left(\frac{1}{3}\right)^{pq}\tau ^{\mu \nu \mu _1\mathrm{}\mu _{2q}}h^{\mu _{2q+1}\mu _{2q+2}}\mathrm{}h^{\mu _{2p1}\mu _{2p}}.$$
(87)
When we put (87) into (82) we obtain the desired expansion of $`\mathrm{\Omega }`$ in terms of the generalized $`\tau `$:

$`\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _k}`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{[\frac{k}{2}]}{}}}\left(\begin{array}{c}2i\\ k\end{array}\right){\displaystyle \underset{j=0}{\overset{i}{}}}\left({\displaystyle \frac{1}{3}}\right)^{ij}`$ (88)
$`\times `$ $`\tau ^{\mu \nu (\mu _1\mathrm{}\mu _{2j}}h^{\mu _{2j+1}}\mathrm{}h^{\mu _{2i1}\mu _{2i}}u^{\mu _{2i+1}}\mathrm{}u^{\mu _k)}.`$
The generalized $`\tau `$ are useful for developing the angular structure of tensors and we find them more suited to this problem than the spherical harmonics of the usual moment schemes.

## 7 Ordering the Summations

In a calculation that is in some ways the predecessor of this one (Unno and Spiegel, 1966), a series of angular moments was rearranged and resummed to provide a closure approximation for the thermal terms. These terms lacked off-diagonal elements and so gave only the Eddington approximation, which is nevertheless a great improvement on the diffusion approximation. The moral is that it is very effective to ensure that all orders in mean free path made themselves felt in the result. However, since the approximation used in that earlier work fails to capture the viscous effects, we here propose one that retains the off-diagonal terms and leads to results good both for long and short photon mean free paths. To bring out the nature of the rearrangement we propose, let us shift the index $`k`$ in (76) to $`k+m`$ so that (75)-(76) becomes

$$\mathrm{\Xi }^{\mu \nu }=\underset{m=0}{\overset{\mathrm{}}{}}\underset{k=0}{\overset{m}{}}\mathrm{{\rm Y}}_{(m,k)}^{\mu \nu }$$
(89)
where

$$\mathrm{{\rm Y}}_{(m,k)}^{\mu \nu }=(1)^{k+m}ϵ^mY_{\mu _1\mathrm{}\mu _{k+m}}^{(m,m+k)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _{k+m}}.$$
(90)
When we write out the $`\mathrm{{\rm Y}}_{m,k}^{\mu \nu }`$ as an array of tensors, with $`m`$ designating the rows and $`k`$ designating the columns, it takes the form

$$\begin{array}{cccc}\mathrm{{\rm Y}}_{(0,0)}^{\mu \nu }& 0& 0& \mathrm{}\\ \mathrm{{\rm Y}}_{(1,0)}^{\mu \nu }& \mathrm{{\rm Y}}_{(1,1)}^{\mu \nu }& 0& \mathrm{}\\ \mathrm{{\rm Y}}_{(2,0)}^{\mu \nu }& \mathrm{{\rm Y}}_{(2,1)}^{\mu \nu }& \mathrm{{\rm Y}}_{(2,2)}^{\mu \nu }& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}$$
(91)
Thus, the $`\mathrm{{\rm Y}}_{(m,k)}^{\mu \nu }`$ form a lower triangular matrix of tensors.

The first sum in (89) is taken over the $`m^{\mathrm{th}}`$ row in the array of $`\mathrm{{\rm Y}}`$s and the second sum is over the columns. In the spirit of the earlier derivation of the Eddington approximation by resummation, we rearrange the terms so that the sum is first taken over columns, then over the rows. That is, we can sum first over the $`m`$. Since, for fixed $`k`$, the sum over $`m`$ begins with the value $`m=k`$, it is then convenient to shift the index $`m`$ by $`k`$. This simple rearrangement amounts to

$$\underset{m=0}{\overset{\mathrm{}}{}}\underset{k=0}{\overset{m}{}}\mathrm{{\rm Y}}_{(m,k)}^{\mu \nu }=\underset{k=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\mathrm{{\rm Y}}_{(m+k,k)}^{\mu \nu }.$$
(92)
This device, which we shall employ again, turns (89) into

$$\mathrm{\Xi }^{\mu \nu }=\underset{k=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\mathrm{{\rm Y}}_{(m+k,k)}^{\mu \nu }$$
(93)
so that the expansion of the viscosity tensor is now

$$\mathrm{\Xi }^{\mu \nu }=\underset{k=0}{\overset{\mathrm{}}{}}ϵ^k\mathrm{\Xi }_{[k]}^{\mu \nu }$$
(94)
with

$$\mathrm{\Xi }_{[k]}^{\mu \nu }=\underset{m=0}{\overset{\mathrm{}}{}}(1)^mϵ^mY_{\mu _1\mathrm{}\mu _{m+2k}}^{(m+k,m+2k)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _{m+2k}}.$$
(95)
The two expansions for $`\mathrm{\Xi }^{\mu \nu }`$, (75) and (94) look similar but, when we compare them carefully, we detect an important difference: only the series for $`\mathrm{\Xi }_{[k]}^{\mu \nu }`$ contains terms of all orders in $`ϵ`$. We see also that (93) involves $`\mathrm{{\rm Y}}_{(m+k,k)}^{\mu \nu }`$, which in turn requires that we be more explicit about the $`Y_{\mu _1\mathrm{}\mu _{m+2k}}^{(m+k,m+2k)}`$ in (95). The relevant development can proceed with the help of the iteration formula (68). We first use (66) in (68) with $`mm1`$ and $`k=m`$ to obtain

$$Y_{\mu _1\mathrm{}\mu _m}^{(m,m)}=Y_{\mu _1\mathrm{}\mu _{m1},\mu _m}^{(m1,m1)}.$$
(96)
Iteration of this formula gives

$$Y_{\mu _1\mathrm{}\mu _m}^{(m,m)}=Y_{,\mu _1\mathrm{}\mu _m}^{(0,0)}.$$
(97)
More generally, for $`k>0`$, we may use $`(\text{68})`$ to find

$$Y_{\mu _1\mathrm{}\mu _{2k}}^{(k,2k)}=(2k+2)Y_{\mu _1\mathrm{}\mu _{2k2}}^{(k1,2k2)}u_{\mu _{2k1},\mu _{2k}}.$$
(98)
Again, using (68), we see that

$$Y_{\mu _1\mathrm{}\mu _{2k+1}}^{(k+1,2k+1)}=Y_{\mu _1\mathrm{}\mu _{2k},\mu _{2k+1}}^{(k,2k)}+(2k+3)Y_{\mu _1\mathrm{}\mu _{2k1}}^{(k,2k1)}u_{\mu _{2k},\mu _{2k+1}},$$
(99)
and so on. For the special case of $`k=0`$, we obtain (97).

Now we can use the rearrangement trick again by noting that if we put (98)-(99) into (95), it takes the form

$$\mathrm{\Xi }_{[k]}^{\mu \nu }=\underset{m=0}{\overset{\mathrm{}}{}}\underset{i=0}{\overset{m}{}}X_{(m+k,k,i)}^{\mu \nu }$$
(100)
where the expressions for the $`X`$ can be obtained by comparing to (95). We do not need their explicit form since we need only to observe that by a repetition of the formula (92) we get

$$\mathrm{\Xi }_{[k]}^{\mu \nu }=\underset{i=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}X_{(m+k+i,k,i)}^{\mu \nu }.$$
(101)
We may then carry out the necessary substitutions, and separate out the exceptional case $`k=0`$, to obtain

$`\mathrm{\Xi }^{\mu \nu }`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(1)^{m+i}ϵ^{m+k+i}(2k+i+2)\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _{2k+m+i}}`$ (102)
$`\times \left[Y_{\mu _1\mathrm{}\mu _{2k+i2}}^{(k+i1,2k+i2)}u_{\mu _{2k+i1},\mu _{2k+i}}\right]_{,\mu _{2k+1+i}\mathrm{}\mu _{2k+m+i}}`$
$`+{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(1)^mϵ^mY_{,\mu _1\mathrm{}\mu _m}^{(0,0)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _m}`$
This completes the formal expansion procedure for the viscosity tensor, which we can now express in terms of the generalized $`\tau `$ and the coefficients of the intensity expansion. For both of these sets of quantities we have recursion formulae. We next begin to make approximations to prepare the way for representing the results in terms of rational approximations to the series.

## 8 Approximating the Orders

So far we have been proceeding with formal expansions in mean free path and making convenient rearrangements of the series obtained. The result is that terms in the series depend on $`ϵ`$ so that terms of similar character can be compared side by side. Now we turn to some approximations to make the series more manageable. What we shall propose in this section is the neglect of those terms in each order of the formal expansion that are clearly (and asymptotically) smaller than their fellows in the same order, as in the earlier study of the radiative heat equation. There it was observed that by retaining a dominant term in every order one is led to a useful outcome (in that case, the Eddington approximation). That approach is what we mean by ‘approximating the orders.’

We note that in (102) we may change the order of summation, since the sums are each taken to infinity in that expression. So we may rewrite (102) more compactly as

$$\mathrm{\Xi }^{\mu \nu }=\underset{m=0}{\overset{\mathrm{}}{}}(ϵ)^m\left\{\underset{k=1}{\overset{\mathrm{}}{}}\underset{i=0}{\overset{\mathrm{}}{}}(1)^iϵ^{k+i}A_{(m,k,i)}^{\mu \nu }+Y_{,\mu _1\mathrm{}\mu _m}^{(0,0)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _m}\right\}$$
(103)
where

$$A_{(m,k,i)}^{\mu \nu }=(2k+i+2)\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _{2k+m+i}}\left[Y_{\mu _1\mathrm{}\mu _{2k+i2}}^{(k+i1,2k+i2)}u_{\mu _{2k+i1},\mu _{2k+i}}\right]_{,\mu _{2k+1+i}\mathrm{}\mu _{2k+m+i}}.$$

(104)
We operate first on the double sum of the summand of (103) by separating out its leading term for $`k=1`$ and $`i=0`$. The double sum can then be rewritten as

$$ϵA_{(m,1,0)}^{\mu \nu }+\underset{k=2}{\overset{\mathrm{}}{}}\underset{i=1}{\overset{\mathrm{}}{}}(1)^iϵ^{k+i}A_{(m,k,i)}^{\mu \nu }+\underset{k=2}{\overset{\mathrm{}}{}}ϵ^kA_{(m,k,0)}^{\mu \nu }+\underset{i=1}{\overset{\mathrm{}}{}}(1)^iϵ^{i+1}A_{(m,1,i)}^{\mu \nu }.$$
(105)
In the double sum of (105), we shift the indices so that $`k`$ becomes $`k+2`$ and $`i`$ goes over into $`i+1`$. Similar shifts can be carried out for the last two terms of sum. Then (105) becomes

$$ϵA_{(m,1,0)}^{\mu \nu }ϵ^3\underset{k=0}{\overset{\mathrm{}}{}}\underset{i=0}{\overset{\mathrm{}}{}}(1)^iϵ^{k+i}A_{(m,k+2,i+1)}^{\mu \nu }+ϵ^2\underset{k=0}{\overset{\mathrm{}}{}}ϵ^k\left[A_{(m,k+2,0)}^{\mu \nu }(1)^kA_{(m,1,k+1)}^{\mu \nu }\right].$$
(106)
Since the last two terms in this expression are down by at least a factor of $`ϵ`$ from the first term, from which they are otherwise not dissimilar, we may neglect them.

Next, let us look at the second term in the summand of (103). A simplification of this term can be based on the fact, seen on inspection, that $`\mathrm{\Omega }^{\mu \nu }=0`$, together with the result given in (78) that $`\mathrm{\Omega }^{\mu \nu \rho }=0`$. With the first two terms gone from this second sum over $`m`$ it becomes

$$\underset{m=0}{\overset{\mathrm{}}{}}(ϵ)^mY_{,\mu _1\mathrm{}\mu _m}^{(0,0)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _m}=\underset{m=2}{\overset{\mathrm{}}{}}(ϵ)^mY_{,\mu _1\mathrm{}\mu _m}^{(0,0)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _m}.$$
(107)
By shifting $`m`$ to $`m+2`$, we can turn (107) into

$$\underset{m=0}{\overset{\mathrm{}}{}}(ϵ)^{m+2}Y_{,\mu _1\mathrm{}\mu _{m+2}}^{(0,0)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _{m+2}},$$
(108)
and so, we can write (103) as

$`\mathrm{\Xi }^{\mu \nu }={\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^m\{ϵA_{(m,1,0)}^{\mu \nu }`$ $``$ $`ϵ^3{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}(1)^iϵ^{k+i}A_{(m,k+2,i+1)}^{\mu \nu }`$ (109)
$`+`$ $`ϵ^2Y_{,\mu _1\mathrm{}\mu _{m+2}}^{(0,0)}\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _{m+2}}\}.`$
As we have stated, we are going to neglect the middle term in this summand and now we see that we may also neglect the last term. Our expression for the radiative viscous stress tensor may then be approximated as

$$\mathrm{\Xi }^{\mu \nu }=4ϵ\underset{m=0}{\overset{\mathrm{}}{}}(ϵ)^m\mathrm{\Omega }^{\mu \nu \rho \sigma \mu _1\mathrm{}\mu _m}[Y^{(0,0)}u_{\rho ,\sigma }]_{,\mu _1\mathrm{}\mu _m}$$
(110)
where each of the neglected terms is $`O(ϵ)`$ compared to its counterparts. We note that the term for $`m=0`$ of this expansion is the Thomas approximation to the stress tensor.

## 9 The Viscous Smoothing Operator

The thermal smoothing operator in (1) is a rational fraction of polynomials in the Laplacian operator. It was obtained by first expanding the original integral operator of radiative smoothing in a series of powers of photon mean free path and approximating each term according to various criteria, especially simplicity of geometrical structure. In attempting a similar approach here for the viscosity problem we encounter a more complicated development. We have grouped terms in our expansion by a combination of order in the mean free path expansion and geometric complexity. This approach has permitted us to compare terms of like character and then to keep in each cluster of terms the leading ones in powers of mean free path. Now we wish to examine the remaining terms and keep only the geometrically simplest terms, as was done in the derivation of $`𝒟`$ of (1). We shall also make some further approximations in terms of $`ϵ`$.

As always, much depends on $`Y^{(0,0)}=S`$ (see (69)), which we now introduce into (110). We then observe that the viscous stress tensor takes the simple form

$$\mathrm{\Xi }^{\mu \nu }=4ϵID^{\mu \nu \rho \sigma }\left[Su_{\rho ,\sigma }\right]$$
(111)
where the smoothing operator, $`ID`$, the centerpiece of this investigation, is given as

$$ID^{\mu \nu \rho \sigma }=\underset{m=0}{\overset{\mathrm{}}{}}(ϵ)^m\mathrm{\Omega }^{\mu \nu \rho \sigma \mu _1\mathrm{}\mu _m}_{\mu _1\mathrm{}\mu _m}.$$
(112)
We now wish to extract the geometrically simplest part from each term in this expansion. To do this, we recall formula (88) for $`\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _m}`$, in which the leading nonvanishing term in the series involves the tensor $`\tau ^{\mu \nu \rho \sigma }`$ that also appears in the standard form of the viscous stress tensor. The origin of the general $`\tau `$ has been discussed in the derivation of (87).

To obtain the usual form of the Eddington approximation, one normally omits higher order angular moments of the intensity field in series such as we have obtained. But the conventionally neglected moments contain projections onto $`\tau ^{\mu \nu \rho \sigma }`$ and so that form of truncation misrepresents the viscosity tensor. That is why we developed the $`\tau `$ tensors — by retaining only the leading angular approximation in these quantities we keep just the angular quantity that is central to the representation of viscosity. Therefore, we here select from (88) only the geometrically most elementary part of the term in each order and write

$$\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _m}=\underset{\mathrm{}=1}{\overset{[\frac{m}{2}]}{}}\left(\begin{array}{c}2\mathrm{}\\ m\end{array}\right)\left(\frac{1}{3}\right)^\mathrm{}1\tau ^{\mu \nu (\mu _1\mu _2}h^{\mu _3\mu _4}\mathrm{}h^{\mu _{2\mathrm{}1}\mu _2\mathrm{}}u^{\mu _{2\mathrm{}1}}\mathrm{}u^{\mu _m)}.$$
(113)
In allowing contributions to each term from only the leading generalized $`\tau `$ we have made what may be called a Galerkin truncation.

Next we make an index shift $`j\mathrm{}1`$, so that (113) becomes

$$\mathrm{\Omega }^{\mu \nu \mu _1\mathrm{}\mu _m}=\underset{j=0}{\overset{[\frac{m2}{2}]}{}}\left(\begin{array}{c}2j+2\\ m\end{array}\right)\left(\frac{1}{3}\right)^j\tau ^{\mu \nu (\mu _1\mu _2}h^{\mu _3\mu _4}\mathrm{}h^{\mu _{2j+1}\mu _{2j+2}}u^{\mu _{2j+3}}\mathrm{}u^{\mu _m)}.$$
(114)
This in turn can be written as

$$\mathrm{\Omega }^{\mu \nu \rho \sigma \mu _1\mathrm{}\mu _m}=\underset{j=0}{\overset{[\frac{m}{2}]}{}}\left(\begin{array}{c}2j+2\\ m+2\end{array}\right)\left(\frac{1}{3}\right)^j\tau ^{\mu \nu (\rho \sigma }h^{\mu _1\mu _2}\mathrm{}h^{\mu _{2j1}\mu _{2j}}u^{\mu _{2j+1}}\mathrm{}u^{\mu _m)}.$$
(115)
With this approximation for $`\mathrm{\Omega }`$, the viscous smoothing operator is

$`ID^{\mu \nu \rho \sigma }`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^m{\displaystyle \underset{j=0}{\overset{[\frac{m}{2}]}{}}}\left(\begin{array}{c}2j+2\\ m+2\end{array}\right)\left({\displaystyle \frac{1}{3}}\right)^j`$ (116)
$`\times `$ $`\tau ^{\mu \nu (\rho \sigma }h^{\mu _1\mu _2}\mathrm{}h^{\mu _{2j1}\mu _{2j}}u^{\mu _{2j+1}}\mathrm{}u^{\mu _m)}_{\mu _1\mathrm{}\mu _m}.`$
Because $`_{\mu _1\mathrm{}\mu _m}`$ is symmetric, the symmetrization of the $`\mu _i`$ is redundant and may be omitted. However, this is not the case for indices $`\rho `$ and $`\sigma `$. If we approximate (116) by keeping only the term with factor $`\tau ^{\mu \nu \rho \sigma }`$ after symmetrization, there is an extra numerical factor $`\frac{2}{(m+2)(m+1)}`$. Taking all these into account and interchanging indices $`j`$ and $`m`$ in the summation, we obtain

$$ID^{\mu \nu \rho \sigma }=\tau ^{\mu \nu \rho \sigma }\underset{j=0}{\overset{\mathrm{}}{}}\underset{m=2j}{\overset{\mathrm{}}{}}\frac{2}{(m+2)(m+1)}\left(\begin{array}{c}2j+2\\ m+2\end{array}\right)(ϵ)^m\left(\frac{1}{3}\right)^jIL_{\{j,m2j\}}$$
(117)
where

$$IL_{\{j,m\}}=h^{\mu _1\mu _2}\mathrm{}h^{\mu _{2j1}\mu _{2j}}u^{\mu _{2j+1}}\mathrm{}u^{\mu _{m+2j}}_{\mu _1\mathrm{}\mu _{m+2j}}.$$
(118)
When we replace the index $`m`$ by $`m+2j`$ (117) becomes

$$ID^{\mu \nu \rho \sigma }=\tau ^{\mu \nu \rho \sigma }\underset{j=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\frac{2}{(m+2)(m+1)}\left(\begin{array}{c}2j+2\\ m+2j+2\end{array}\right)(ϵ)^{m+2j}\left(\frac{1}{3}\right)^jIL_{\{j,m\}}.$$
(119)
We next want to simplify the differential operator $`IL`$ by approximating it with

$$_{\{j,m\}}=h^{\mu _1\mu _2}_{\mu _1\mu _2}\mathrm{}h^{\mu _{2j1}\mu _{2j}}_{\mu _{2j1}\mu _{2j}}u^{\mu _{2j+1}}_{\mu _{2j+1}}\mathrm{}u^{\mu _{m+2j}}_{\mu _{m+2j}}.$$
(120)
Direct evaluation shows that $`IL_{\{0,0\}}=_{\{0,0\}},`$ $`IL_{\{1,0\}}=_{\{1,0\}}`$ and $`IL_{\{0,1\}}=_{\{0,1\}},`$ and this means that we would be making approximations only in the higher order terms in (119) in replacing $`IL`$ by $``$. That is, this replacement leads to the neglect of terms of relative order $`ϵ^2`$ in only the higher terms, as we detail in Appendix B.

The advantage in making the approximate replacement of $`IL`$ is that we may write

$$_{\{j,m\}}=(h^{\alpha \beta }_{\alpha \beta })^j(u^\alpha _\alpha )^m;$$
(121)
this lends itself to the rational approximations that we wish to develop for the viscous smoothing operator, which we may now write as

$$ID^{\mu \nu \rho \sigma }=\tau ^{\mu \nu \rho \sigma }\underset{j=0}{\overset{\mathrm{}}{}}\left(\frac{ϵ^2}{3}h^{\alpha \beta }_{\alpha \beta }\right)^j\underset{m=0}{\overset{\mathrm{}}{}}\frac{2}{(m+2)(m+1)}\left(\begin{array}{c}2j+2\\ m+2j+2\end{array}\right)(ϵu^\alpha _\alpha )^m.$$
(122)
## 10 Rational Approximation for $`ID^{\mu \nu \rho \sigma }`$

An asymptotic expansion such as we have developed can be made more useful by the technique of resummation (Baker, 1975; Bender, 1978), of which there are various forms. The general idea is to replace the series by a rational function of operators in one of several well documented ways. For a given function $`g(ϵ)`$ with series representation

$$g(ϵ)=\underset{i=0}{\overset{\mathrm{}}{}}\gamma _iϵ^i$$
(123)
the Padé approximation of order $`[m,n]`$ is written as the rational function

$$\mathrm{}_{[m,n]}[g(ϵ)]\frac{P_m(ϵ)}{Q_n(ϵ)}$$
(124)
where $`P_m(ϵ)`$ and $`Q_n(ϵ)`$ are polynomials in $`ϵ`$ of degrees $`m`$ and $`n`$, respectively. The coefficients in these polynomials are to be fixed according to the condition that the first $`m+n`$ coefficients of the expansion of $`\mathrm{}_{[m,n]}`$ should match the first $`m+n`$ coefficients in the series for $`g(ϵ)`$.

This procedure, for all its wide acceptance, has not been rigorously justified, but a considerable experience shows that most of the reasonable choices are qualitatively better than a truncated asymptotic expansion. An intuitive explanation is that the rational function can have poles at the same points that the function to be approximated does. Because of this, the range of validity of the expansion may be extended.

In a simple version of resummation, we might try to represent only a few terms in $`ID^{\mu \nu \rho \sigma }`$, say out to order $`ϵ^2`$, as a rational fraction of operators. The first few terms from (122) are

$$ID^{\mu \nu \rho \sigma }=\tau ^{\mu \nu \rho \sigma }\left[1ϵu^\alpha _\alpha +(ϵu^\alpha _\alpha )^2\frac{ϵ^2}{3}h^{\alpha \beta }_{\alpha \beta }+\mathrm{}\right].$$
(125)
The $`[0,1]`$ order Padé approximation is

$$ID^{\mu \nu \rho \sigma }=\frac{\tau ^{\mu \nu \rho \sigma }}{1+ϵu^\alpha _\alpha }$$
(126)
while the order $`[0,2]`$ approximation can be expressed as

$$ID^{\mu \nu \rho \sigma }=\frac{\tau ^{\mu \nu \rho \sigma }}{[1+ϵu^\alpha _\alpha +\frac{1}{3}ϵ^2h^{\kappa \lambda }_{\kappa \lambda }]}.$$
(127)
These representations of the viscous stress tensor should already provide a considerable improvement on the standard diffusive form, but we would like to suggest what we think is yet a better version. In the earlier work on the thermal smoothing problem, it was possible to obtain an operator that had very accurate limits at large and small $`ϵ`$. In the normal case, where the transition between the two situations is smooth, this meant that the approximation would be reasonably good for all values of $`ϵ`$. We would like to achieve a similar accuracy for viscous smoothing and, since (127) is always good for vanishing $`ϵ`$, we need to focus on large $`ϵ`$ in our determination. For this, we propose the following generalization of the prescription for fixing the numerical coefficients in the rational approximation.

We introduce a new rational approximation

$$\mathrm{}_{[m,n]}^{(N)}[g(ϵ)]\frac{P_m(ϵ)}{Q_n(ϵ)}$$
(128)
where, as before, $`P_m(ϵ)`$ and $`Q_n(ϵ)`$ are polynomials in $`ϵ`$ whose orders are to be specified at our convenience. However, this time, the coefficients are to be chosen according to a more general matching prescription. We suppose that the asymptotic series for $`\mathrm{}_{[m,n]}^{(N)}`$ is

$$\mathrm{}_{[m,n]}^{(N)}[g(ϵ)]=\underset{i=0}{\overset{\mathrm{}}{}}R_i^{(N)}ϵ^i$$
(129)
and we determine its coefficients by the conditions that $`R_i^{(N)}=\gamma _i`$ for $`i=0,1,,\mathrm{},m1`$ and for $`i=N+m,\mathrm{},N+n+m1`$ where the orders $`m`$ and $`n`$ are to be specified as before. Thus we leave a gap between the highest and lowest order terms that we choose to match with the asymptotic series of the original function.

By varying $`N`$ we control how much influence the higher terms in the series can have on our rational approximation. On the other hand, nothing has been lost since it is clear that when $`N=0`$ this procedure reduces to the standard one: $`\mathrm{}_{[m,n]}^{(0)}=\mathrm{}_{[m,n]}`$. However, we shall be especially interested in

$$R_{[m,n]}^{(\mathrm{})}[g(ϵ)]=\underset{N\mathrm{}}{lim}R_{[m,n]}^{(N)}[g(ϵ)].$$
(130)
For example, if we represent $`g(ϵ)`$ at the level of approximation where we have

$$\mathrm{}_{[0,1]}^{(N)}[g(ϵ)]=\frac{\gamma _0}{1q_Nϵ},$$
(131)
we obtain

$$q_N=\left(\frac{\gamma _N}{\gamma _0}\right)^{1/N}$$
(132)
and

$$q_N=\frac{\gamma _{N+1}}{\gamma _N}.$$
(133)
We have then, in particular, that

$$\mathrm{}_{[0,1]}^{(\mathrm{})}[g(ϵ)]=\frac{\gamma _0}{1q_{\mathrm{}}ϵ}$$
(134)
with

$$q_{\mathrm{}}=\underset{N\mathrm{}}{lim}\left(\frac{\gamma _N}{\gamma _0}\right)^{1/N}.$$
(135)
Since $`\underset{N\mathrm{}}{lim}(\gamma _0)^{1/N}=1`$ for any finite $`\gamma _0`$, we have

$$q_{\mathrm{}}=\underset{N\mathrm{}}{lim}(\gamma _N)^{1/N}=\underset{N\mathrm{}}{lim}\frac{\gamma _{N+1}}{\gamma _N}.$$
(136)
For the present study, we shall use only the representation for the viscous smoothing operator with $`N=\mathrm{}`$ and $`[m,n]=[0,1]`$. It is therefore convenient to use a simple notation for this case, namely

$$[\mathrm{}]\mathrm{}_{[0,1]}^{(\mathrm{})}[\mathrm{}].$$
(137)
Returning now to the study of the viscous smoothing operator in (122), we apply the resummation procedure to find that

$$\left[\underset{n=0}{\overset{\mathrm{}}{}}\left(\begin{array}{c}2j+2\\ n+2j+2\end{array}\right)(ϵu^\alpha _\alpha )^n\right]=\frac{1}{1+ϵu^\alpha _\alpha }.$$
(138)
The disappearance of $`j`$ from this expression is a result of the large $`N`$ limit.

Then we note that

$$\underset{j=0}{\overset{\mathrm{}}{}}ϵ^{2j}(\frac{1}{3}h^{\alpha \beta }_{\alpha \beta })^j=\frac{1}{1+\frac{1}{3}ϵ^2h^{\alpha \beta }_{\alpha \beta }}.$$
(139)
In terms of these newly defined smoothing operators, we may write

$$ID^{\mu \nu \rho \sigma }=\tau ^{\mu \nu \rho \sigma }𝒟_x𝒟_\tau $$
(140)
where

$$𝒟_x=\frac{1}{1+\frac{1}{3}ϵ^2h^{\alpha \beta }_{\alpha \beta }},𝒟_\tau =\frac{1}{1+ϵu^\alpha _\alpha }.$$
(141)
Then the viscous stress tensor can be written as

$$\mathrm{\Xi }^{\mu \nu }=\frac{16\pi ϵ}{15}\tau ^{\mu \nu \rho \sigma }𝒟_x𝒟_\tau (Su_{\rho ,\sigma })$$
(142)
or, more compactly, within our present level of approximation, as

$$\mathrm{\Xi }^{\mu \nu }=𝒟_x𝒟_\tau [\mu _0\tau ^{\mu \nu \rho \sigma }u_{\rho ,\sigma }]$$
(143)
where $`\mu _0=\frac{16\pi ϵS}{15}`$ is the first order shear viscosity of the Thomas theory.

If we denote the Thomas stress tensor of (33) by

$$\mathrm{\Xi }_T^{\mu \nu }=\mu _0\tau ^{\mu \nu \rho \sigma }u_{\rho ,\sigma },$$
(144)
then, from (143), (138) and (139), we arrive at this linear partial differential equation for $`\mathrm{\Xi }^{\mu \nu }`$:

$$(1+\frac{1}{3}ϵ^2h^{\alpha \beta }_{\alpha \beta })(1+ϵu^\alpha _\alpha )\mathrm{\Xi }^{\mu \nu }=\mathrm{\Xi }_T^{\mu \nu }.$$
(145)
This is an evolution equation or, what meteorologists would call, a prognostic equation. However, the operator $`u^\mu _\mu `$ is effectively a time derivative and, for times longer than the flight times of photons, it is negligible in this formula. So when we are dealing with nonrelativistic motions, (145) reduces to the diagnostic equation

$$\frac{1}{3}ϵ^2h^{\alpha \beta }_{\alpha \beta }\mathrm{\Xi }^{\mu \nu }=\mathrm{\Xi }_T^{\mu \nu }\mathrm{\Xi }_{\mu \nu }.$$
(146)
An illustration of an application of this formula is given elsewhere (Chen, PhD thesis).

In the nonrelativistic case, (146) could be written as

$$\frac{1}{3}ϵ^2\mathrm{\Delta }\mathrm{\Xi }^{\mu \nu }=\mathrm{\Xi }_T^{\mu \nu }\mathrm{\Xi }^{\mu \nu }$$
(147)
where $`\mathrm{\Delta }`$ is the three-dimensional Laplacian operator. However, this last step involves the neglect of retardation times, which is not always a safe practice (Delache and Froeschlé, 1972).

## 11 Transport Coefficients

Our formula for the radiative stress tensor provides a closure relation that may be used in the standard evolution equations for the radiative energy density and flux. These equations, in conjunction with the fluid dynamical equations, provide a two-fluid description of radiative fluid dynamics (Simon, 1963; Hsieh and Spiegel, 1976) that is reasonably complete. However, the simpler one-fluid formulation is more tractable and could be used in certain cases, as when radiative forces are weak.

In the single fluid approximation we may think of the radiation as providing a source of dissipation for the fluid. In that case, we would treat the fluid as imperfect with appropriate dissipative terms and use the transfer theory to provide the corresponding transport coefficients. This version of radiative fluid dynamics has been developed by Weinberg (1972). We add to this theory here by providing some resummed results for the coefficients of viscosity and conductivity.

Since the procedures we have laid down provide formal solutions to the transfer equation, we may compute the radiative quantities in terms of the material properties of the radiating medium. We have described the calculation of the viscous stress tensor since that is the object needed to close the hierarchy of moment equations at a low order. We can similarly solve for the energy density and radiative flux.

On resumming the expansions for the energy density and flux of the radiation we obtain expressions for them like that for the viscosity tensor. We shall not give the calculations for these other moments since they parallel those given above and, in any case, our intention in this section is merely to indicate this other possible direction for treating the dynamics of radiating fluids. So we simply report that

$$E=\frac{4\pi }{3}𝒟_x𝒟_\tau \left[\left(34ϵu_{,\rho }^\rho \right)S\right]$$
(148)
and

$$F^\mu =\frac{4\pi ϵ}{3}𝒟_x𝒟_\tau \left[h^{\mu \rho }\left(S_{,\rho }4Su^\sigma u_{\rho ,\sigma }\right)\right].$$
(149)
These expressions may be introduced into the radiative stress tensor, which we write here as

$$T^{\mu \nu }=T_0^{\mu \nu }+\mathrm{\Delta }T^{\mu \nu }$$
(150)
where the equilibrium part is given by

$$T_0^{\mu \nu }=4\pi Su^\mu u^\nu \frac{4\pi }{3}Sh^{\mu \nu }$$
(151)
and the nonequilibrium contribution is

$$\mathrm{\Delta }T^{\mu \nu }=(E4\pi S)u^\mu u^\nu +F^\mu u^\nu +F^\nu u^\mu \frac{1}{3}(E4\pi S)h^{\mu \nu }+\mathrm{\Xi }^{\mu \nu }.$$
(152)
With the formula for $`E`$ just given, we find that

$$E4\pi S=\frac{4\pi ϵ}{3}𝒟_x𝒟_\tau (3u^\alpha S_{,\alpha }+4Su_{,\rho }^\rho ),$$
(153)
while $`F^\mu `$ is given by (149).

The transport coefficients can be identified from the expression for entropy generation. This is given by Weinberg (1972) as

$$S_{,\mu }^\mu =\mathrm{\Delta }T^{\mu \nu }\left(\frac{u_{\mu ,\nu }}{T}\frac{u_\mu T_{,\nu }}{T^2}\right)$$
(154)
where $`T`$ is temperature, the entropy current is

$$S^\mu =su^\mu T^1u_\nu \mathrm{\Delta }T^{\mu \nu }$$
(155)
and $`s`$ is the total (matter plus radiation) entropy density in a frame locally moving with the matter, as advocated by Mihalas (1984). (We have made some sign modifications from Weinberg to adapt these expressions to the metric signature used here.)

The calculations needed in the evaluation of the entropy generation are not difficult, but they are lengthy and not central to this work, so we omit them here and merely quote the result, which is

$`S_{,\mu }^\mu `$ $`=`$ $`\eta h^{\mu \nu }\left({\displaystyle \frac{T_{,\mu }}{T}}u^\rho u_{\mu ,\rho }\right)\left({\displaystyle \frac{T_{,\nu }}{T}}u^\rho u_{\nu ,\rho }\right)+\xi \theta ^2`$ (156)
$`+`$ $`{\displaystyle \frac{2\mu }{T}}h^{\mu \rho }h^{\nu \sigma }(u_{\mu ,\nu }+u_{\nu ,\mu }{\displaystyle \frac{2}{3}}h_{\mu \nu }\theta )(u_{\rho ,\sigma }+u_{\sigma ,\rho }{\displaystyle \frac{2}{3}}h_{\rho \sigma }\theta )`$
where $`\theta =u_{,\mu }^\mu `$ and the transport coefficients are given by

$$\eta =𝒟_x𝒟_\tau \left(\frac{16\pi ϵS}{3T}\right),$$
(157)
$$\xi =𝒟_x𝒟_\tau \left[\frac{16\pi ϵS}{T}\left(\frac{u^\alpha T_{,\alpha }}{T\theta }+\frac{1}{3}\right)^2\right]$$
(158)
and

$$\mu =𝒟_x𝒟_\tau \left(\frac{4\pi }{15}ϵS\right).$$
(159)
The positivity of the three generalized transport coefficients is a sufficient condition for the entropy generation rate to be positive. This is of interest in that a failure of this rate to be positive would signal a breakdown in the procedures.

The first term in (156) represents entropy generation by thermal processes — heat flux, momentum exchange with the radiation field — and the corresponding transport coefficient, $`\eta `$, is a generalized thermal conductivity. The second term comes from entropy generation due to the bulk volume expansion, and its coefficient, $`\xi `$, is similarly a form of bulk viscosity. Finally, the third term results from dissipation by shearing motions, with $`\mu `$ as a generalization of the Thomas formula. Each of these transport coefficients take the form of a bare coefficient (the coefficient obtained in the diffusive limit) renormalized by the operation of $`𝒟_x𝒟_\tau `$. However, since these formulae are to be used in connection with a stress tensor of the diffusive form, these results are applicable only when the flow does not oscillate rapidly in time or space.

## 12 Conclusion

In the moment method of radiative transfer (Thorne, 1981; Struchtrup, 1997), one replaces the transfer equation by an infinite hierarchy of moment equations. To render this system usable, one truncates it by introducing a supplementary closure relation among the retained moments. That approximation is typically obtained by truncating the expansion for the intensity at an order corresponding to the level of closure desired. When we do this at the lowest reasonable level, only the first two moment equations are retained, and these need to be closed by an expression for the pressure tensor. The standard procedure at this level, the Eddington approximation, fails badly at representing viscous effects in radiative fluid dynamics (Anderson and Spiegel, 1972).

Here, we are proposing to keep the first two moment equations since these contain the main physics of the radiative action on the fluid. But rather than following the procedure of truncating the moment hierarchy with the corresponding Eddington approximation, we seek a closure approximation by another route. Our procedure is to develop an asymptotic expansion for the intensity in the manner of the Hilbert expansion of kinetic theory (Grad, 1963; Caflisch, 1983), but we use this solution only for the purpose of devising an expression for the viscosity tensor; the detailed solution itself does not play a role in the theory. Since the Hilbert expansion is an outer expansion, in the sense of asymptotic theory (Bender, 1978), we extend its validity by the technique of resummation so as to render it useful for situations where the photon mean free path is long

Because of the structure of the transfer equation, our series for the distribution function turns into a development of $`ϵn^\mu _{,\mu }`$, the operator that appears in the transfer equation (20) and that operates iteratively on the source function in the expansions. To extract the desired information from this development we needed to perform integrals over the momentum space of the photons and this required extensive manipulations. Since the expansion parameter $`ϵ`$ appears in front of the highest (indeed the only) derivatives in the problem, we were faced with a problem in singular perturbation theory, a subject where successful approximation techniques are known but whose rigor is rarely, if ever, established.

Even with the help of the asymptotic methods, and resummation, the calculations were demanding and we were led to idealize the problem to the case with constant mean free path, though this approximation could be avoided by a suitable choice of optical coordinates, as we noted. The resummations, in conjunction with approximations, permitted us to obtain a compact expression for the radiative viscosity tensor for use at for all photon free paths. Our form for this tensor can be written as an integral operator on the diffusion tensor of Thomas (1930). Since this integral operator is simple, we could see at once how to convert the result into the linear partial differential equation (145) for $`\mathrm{\Xi }^{\mu \nu }`$. Depending on whether the fluid dynamics is very relativistic or not, this equation can be either prognostic or diagnostic but, in any case, it provides a general closure relation for radiative fluid dynamics.

Asymptotic procedures, when suitably applied, can give very good results, and they have proved their usefulness in kinetic theory (Caflisch, 1983). But of course the results need to be tested, for example, by comparison with the results of other methods. To do this, we must solve specific problems with the methods being examined and to compare results with those from other methods. A first test of this approach would be to compare it with the moment method for this problem since that theory has been under serious development (Thorne, 1981; Struchtrup, 1997) and should provide a good benchmark for new procedures. We know already that to get anything like a reasonable representation of viscous effects in a moment method, we need to truncate at a higher order than in the standard Eddington approximation (Anderson and Spiegel, 1972).

Of course, both these methods give macroscopic descriptions of radiative fluid dynamics and so cannot satisfy the full boundary and initial conditions of transport theory, as is known in kinetic theory (Grad, 1963). But we suggest that this loss is more than compensated by a gain in flexibility and usefulness of our macroscopic description. This we have seen in applying our method to the problem of entropy generation in cosmology (to be reported elsewhere). Those calculations involve relatively smooth flows and are not as demanding of a method as more complicated radiating flows such as may arise in the dynamics of photon bubbles or radiative interaction with vortices and flux tubes in hot media. We believe that the study of such more complicated flows involving strong, possibly turbulent, fluctuations in radiating media will provide the real test of our approach.

Finally, we would suggest that the methods we have used here could be useful in other problems where the mean free path of the basic elements is long compared to scales of variation in the medium. An example is turbulent convection where the local mixing length theory encounters a loss of validity because the mixing length may be comparable to scales of variation in space and time in the medium. The results derived here may point the way to a more useful theory of eddy transport processes.

## Appendix A Useful Angular Integrals

When we introduce the angular expansions into the integrals for the moments of the intensity we naturally encounter a variety of integrals over products of the $`l^\mu `$ and $`n^\mu `$ alone as well as over various mixed products of them. In this appendix we list some formulas for such products that we have had need of throughout this work. Such formulas may be obtained by writing them in terms of linear combinations of $`u^\mu u^\nu `$, $`\eta ^{\mu \nu }`$, and so on, and using special cases to find the coefficients as in previous studies (Anderson and Spiegel, 1972; Thorne, 1981; Struchtrup, 1997, 1998).

First we note that the angular integral of an odd number of factors $`l^\mu `$ vanishes. Therefore, for the general moment over factors of $`l^\mu `$ we consider (81) only for even $`\mathrm{}`$. As can be seen, $`M^{\mu _1\mathrm{}\mu _2\mathrm{}}`$ is orthogonal to $`u^\mu `$ in all suffices and this simplifies the problem.

The special cases for $`\mathrm{}=1`$ and $`\mathrm{}=2`$ are most commonly encountered and so it is worth writing them out explicitly:

$$l^\mu l^\nu 𝑑\mathrm{\Omega }=\frac{4\pi }{3}h^{\mu \nu }$$
(A1)
$$l^\mu l^\nu l^\rho l^\sigma 𝑑\mathrm{\Omega }=\frac{4\pi }{15}(h^{\mu \nu }h^{\rho \sigma }+h^{\mu \rho }h^{\nu \sigma }+h^{\mu \sigma }h^{\nu \rho }).$$
(A2)
The general case is found by induction and we have

$$M^{\mu _1\mathrm{}\mu _2\mathrm{}}=(1)^{\mathrm{}}\frac{4\pi }{2\mathrm{}+1}h^{(\mu _1\mu _2}\mathrm{}h^{\mu _{2\mathrm{}1}\mu _2\mathrm{})},$$
(A3)
where the symmetrization is assumed for the indices inside the parentheses.

We also make use of moments like

$$N^{\mu _1\mu _2\mathrm{}\mu _k}=n^{\mu _1}n^{\nu _2}\mathrm{}n^{\mu _k}𝑑\mathrm{\Omega },$$
(A4)
the simplest of which is

$$n^\mu n^\nu 𝑑\mathrm{\Omega }=\frac{4\pi }{3}(4u^\mu u^\nu \eta ^{\mu \nu })=\frac{4\pi }{3}(3u^\mu u^\nu h^{\mu \nu }).$$
(A5)
The general case is easily obtained using the binomial theorem. On noting that

$$n^{\mu _1}\mathrm{}n^{\mu _k}=\underset{i=0}{\overset{k}{}}\left(\begin{array}{c}i\\ k\end{array}\right)u^{(\mu _1}\mathrm{}u^{\mu _i}l^{\mu _{i+1}}\mathrm{}l^{\mu _k)}.$$
(A6)
where binomial coefficient is $`\left(\begin{array}{c}i\\ k\end{array}\right)=\frac{k!}{i!(ki)!}`$ we have

$$N^{\mu _1\mu _2\mathrm{}\mu _k}=\underset{i=0}{\overset{k}{}}\left(\begin{array}{c}i\\ k\end{array}\right)u^{(\mu _1}u^{\mu _2}\mathrm{}u^{\mu _i}M^{\mu _{i+1}\mathrm{}\mu _k)}.$$
(A7)
## Appendix B The Operators $`IL`$ and $``$

In Section 9 of the text, we have used the following approximation

$$ID\underset{j=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}C_{\{j,m\}}(ϵ)^{m+2j}IL_{\{j,m\}}=\underset{j=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}C_{\{j,m\}}(ϵ)^{m+2j}_{\{j,m\}}$$
(B1)
where $`C_{\{j,m\}}`$ represents the rest coefficient in operator $`ID`$. To justify this approximation, let us define the difference between $`IL`$ and $``$ as

$$L_{\{j,m\}}=IL_{\{j,m\}}_{\{j,m\}}.$$
(B2)
Evidently (as we saw in the text), $`L_{\{0,0\}}=L_{\{1,0\}}=L_{\{0,1\}}=0`$. This leads to

$`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^{m+2j}C_{\{j,m\}}L_{\{j,m\}}`$ $`=`$ $`{\displaystyle \underset{j=2}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^{m+2j}C_{\{j,m\}}L_{\{j,m\}}`$ (B3)
$`+{\displaystyle \underset{m=2}{\overset{\mathrm{}}{}}}(ϵ)^mC_{\{0,m\}}L_{\{0,m\}}`$
$`+{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^{m+2}C_{\{1,m\}}L_{\{1,m\}}.`$
After shifting indices as we did previously, we get

$`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^{m+2j}C_{\{j,m\}}L_{\{j,m\}}={\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^{m+2j}ϵ^2C_{\{j+1,m\}}L_{\{j+1,m\}}`$
$`+{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^mϵ^2(C_{\{0,m+2\}}L_{\{0,m+2\}}+C_{\{1,m\}}L_{\{1,m\}}).`$ (B4)
We can then decompose the summation for $`_{\{j,m\}}`$ as

$$\underset{j=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}(ϵ)^{m+2j}C_{\{j,m\}}_{\{j,m\}}=\underset{j=1}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}(ϵ)^{m+2j}C_{\{j,m\}}_{\{j,m\}}+\underset{m=0}{\overset{\mathrm{}}{}}(ϵ)^mC_{\{0,m\}}_{\{0,m\}}).$$
(B5)
Thus with the decomposition of (B4)-(B5) we have

$`ID`$ $`=`$ $`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}C_{\{j,m\}}(ϵ)^{m+2j}[_{\{j,m\}}+L_{\{j,m\}}]`$
$`=`$ $`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^{m+2j}C_{\{j,m\}}(_{\{j,m\}}+ϵ^2L_{\{j,m\}})`$
$`+`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^m\left[C_{\{0,m\}}_{\{0,m\}}+ϵ^2(C_{\{0,m+2\}}L_{\{0,m+2\}}+C_{\{1,m\}}L_{\{1,m\}})\right].`$
We see that the term by term difference between $`IL`$ and $``$ is of order $`ϵ^2`$. If we omit in each term the higher correction associated with $`L`$, we are left with

$`ID`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^{m+2j}C_{\{j,m\}}_{\{j,m\}}+{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^mC_{\{0,m\}}_{\{0,m\}}`$ (B7)
$`=`$ $`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(ϵ)^{m+2j}C_{\{j,m\}}_{\{j,m\}}.`$
This is the approximation we used in the text.


format: text

text:
# On the pseudogap phase in high-𝑇_𝑐 superconductors

## Abstract

We describe the approach of the superconducting state as a sequence of cross-over phenomena. As the temperature is decreased, uncorrelated pairing of the electrons leads to the opening of a pseudogap at $`T_F^{}`$. Upon further lowering the temperature those electron pairs acquire well behaved itinerant features at $`T_B^{}`$, leading to partial Meissner screening and Drude type behavior of the optical conductivity. Further decrease of the temperature leads to their condensation and superconductivity at $`T_c`$. The analysis is done on the basis of the Boson-Fermion model in the cross-over regime between 2D and 3D.

It is now generally accepted that the onset of superconductivity in high temperature superconductors (HTS) is controlled by phase fluctuations, whilst a finite amplitude of the order parameter is expected to persist up to some $`T^{}`$, which can be well above $`T_c`$. It has also been argued that the shrinking of the pseudogap phase (the interval $`[T_c,T^{}]`$) with increased doping should be related to an increase of the superfluid density $`n_s`$ rather than of the coherence length. In that case, doping of HTS should not induce a cross-over between a Bose-Einstein condensation (BEC) of preformed pairs and a BCS state and hence the opening of the pseudogap should be unrelated to the onset of superconducting fluctuations. These expectations seem to be confirmed by specific heat data, transient Meissner screening and Andreev spectroscopy. The picture evolves according to which a pseudogap opens up, exclusively driven by amplitude fluctuations (dynamical pair-formation, uncorrelated in space and not necessarily related to incipient superconductivity). At some lower temperature short-range and short-time correlations of the electron pairs set in, eventually driving the system into a superconducting state.

From a theoretical point of view these are questions which can be addressed without having to resort to a specific microscopic mechanism for electron pairing and can be studied on such models as the generalized BCS hamiltonian, the negative U Hubbard model or the Boson-Fermion model (BFM). These phenomenological models which explicitly incorporate pair correlations, capture a number of normal state properties of the HTS.

In a classical BCS type superconductor long range phase coherence occurs as soon as short range amplitude correlations (in form of electron pairing) set in. A simple mean field treatment can perfectly describe this situation. On the contrary, in HTS amplitude and phase correlations are separated and a minimal treatment must contain the possibility of allowing for propagating modes of electron-pairs on a time scale of the order of the inverse zero temperature gap. We shall here explore the intricate relations between the onset of pairing of the electrons, their becoming itinerant and ultimately their condensation. In order to separate these various features in approaching the superconducting state we study them in view of a dimensional cross-over such as to capture part of the doping induced changes when going from the under into the overdoped regime. The dimensional cross-over will be controlled by the degree of anisotropy of the electron dispersion. The anisotropy in the electric transport coefficients lends itself to such a picture.

This study is based on the BFM, a phenomenological model which is particularly suited for that purpose. The notion of rather long lived short range pair correlations is introduced explicitly into this model by assuming the presence of localized tightly bound electron pairs which hybridize with pairs of itinerant electrons. We assume a system described by effective sites each of which can alternatively be occupied by either a bound electron-pair or a pair of itinerant electrons, uncorrelated with each other. For the present study we shall assume a homogeneous distribution of both the bound electron-pairs and the itinerant electrons. It is however feasible that the distribution between those constituents occurs in form of a phase-separated phase with dynamically fluctuating regions of predominantly bound electron-pairs and regions with predominantly electrons, leading to “stripe phases”. Previous studies of the BFM, based on conserving diagrammatic approximations let one envisage a well defined finite temperature interval for the pseudogap phase in 1D and 2D, while for 3D such a phase should be restricted to a very narrow temperature regime above $`T_c`$. This suggests that the nature of the pseudogap phase should depend on dimensionality and that the pair-fluctuations which control it have, depending on temperature, different characteristics as far as superconducting fluctuations are concerned. We expect that at least two energy scales should be involved here, identifiable as two temperatures: $`T_F^{}T^{}`$ and $`T_B^{}T_F^{}`$,

* $`T_F^{}`$ corresponding to the opening of the pseudogap in the DOS of the electrons due to the onset of strong local correlations leading to electron pairing, uncorrelated in space.
* $`T_B^{}`$ corresponding to those local pair states becoming well defined itinerant excitations due to the onset of temporal and spatial correlations.

We shall determine the dependence of $`T_F^{}`$ and $`T_B^{}`$ on dimensionality by monitoring the anisotropy of the ratio of the electron mass orthogonal to the basal plane to that within it, $`\alpha =m_{}/m_{}`$ in the bare electron dispersion

$`\epsilon _𝐤={\displaystyle \frac{1}{m_{}}}\left({\displaystyle \frac{|𝐤_{}|}{\pi }}\right)^2+{\displaystyle \frac{1}{m_{}}}\left({\displaystyle \frac{|𝐤_{}|}{\pi }}\right)^2\mu `$ (1)
The BFM is then defined by the following Hamiltonian

$`H={\displaystyle \underset{𝐤,\sigma }{}}\epsilon _𝐤c_{𝐤\sigma }^{}c_{𝐤\sigma }+E_0{\displaystyle \underset{𝐪}{}}b_𝐪^{}b_𝐪+v{\displaystyle \underset{𝐤}{}}[b_𝐪^{}c_{𝐤+𝐪}c_𝐤+h.c.],`$ (2)
with $`E_0=\mathrm{\Delta }_B2\mu `$ and $`v`$ denoting the pair-exchange coupling constant. $`c_{𝐤\sigma }^{()}`$ denote fermionic operators for electrons with spin $`\sigma `$ and wavevector $`𝐤`$ and $`b_𝐪^{()}`$ describe tightly bound electron pairs which will be considered as simple Bosons. The chemical potential $`\mu `$ is common to Fermions and Bosons (up to a factor 2 for the Bosons) in order to guarantee charge conservation. We examine this model within the lowest order self-consistent conserving diagrammatic approximation for which the self-energies for the Fermions and Bosons are given by:

$`\mathrm{𝐼𝑚}\mathrm{\Sigma }_B(𝐪,\omega )={\displaystyle \frac{v^2}{2\pi }}{\displaystyle \frac{d^3k}{(2\pi )^3}𝑑\epsilon \mathrm{𝐼𝑚}G_F(𝐤𝐪,\omega \epsilon )}`$ (3)
$`\mathrm{𝐼𝑚}G_F(𝐤,\epsilon )\left(tanh({\displaystyle \frac{\epsilon }{2k_BT}})tanh({\displaystyle \frac{\epsilon \omega }{2k_BT}})\right)`$ (4)
$`\mathrm{𝐼𝑚}\mathrm{\Sigma }_F(𝐤,\omega )={\displaystyle \frac{v^2}{2\pi }}{\displaystyle \frac{d^3q}{(2\pi )^3}𝑑\epsilon \mathrm{𝐼𝑚}G_F(𝐤𝐪,\epsilon \omega )}`$ (5)
$`\mathrm{𝐼𝑚}G_B(𝐪,\epsilon )\left(tanh({\displaystyle \frac{\epsilon \omega }{2k_BT}})cotanh({\displaystyle \frac{\epsilon }{2k_BT}})\right).`$ (6)
$`G_F(𝐤,\omega )=[\omega \epsilon _𝐤\mathrm{𝑅𝑒}\mathrm{\Sigma }_F(𝐤,\omega )i\mathrm{𝐼𝑚}\mathrm{\Sigma }_F(𝐤,\omega )]^1`$ and $`G_B(𝐪,\epsilon )=[\omega E_0\mathrm{𝑅𝑒}\mathrm{\Sigma }_B(𝐪,\omega )i\mathrm{𝐼𝑚}\mathrm{\Sigma }_B(𝐪,\omega )]^1`$ denote the fully renormalized Fermion and Boson Green’s functions which have to be determined selfconsistently. The above set of equations for the Green’s functions is solved by a standard iterative procedure. The integration over the momenta $`𝐤_{}`$ is carried out by summing over a grid of 1025 points in the interval $`[0,\pi ]`$. The integration over frequencies is carried out slightly above the real axis i.e., taking $`\omega =\mathrm{𝑅𝑒}\omega +i\eta `$ with $`\eta =0.01`$ and by summation over 2048 real frequencies in the interval $`[2,2]`$.Throughout the present work all energies are given in units of the bare basal plane bandwidth $`\stackrel{~}{D}=8t_{}`$.

The pseudogap features which result from the solution of the above set of equations depend strongly on dimensionality. In 1D and 2D the Fermionic self-energy is essentially determined by the term proportional to $`cotanh`$ which leads to the most divergent contribution. For 3D the major contribution to the self-energy comes from the term proportional to $`tanh`$ and does not give rise to a noticeable pseudogap effect. A pseudogap is only seen as we lower the temperature and approach $`T_c`$ where the term proportional to $`cotanh`$ again becomes important.

The various parameters entering our Hamiltonian are determined such as to reproduce certain robust features of HTS. Pinning down the number of Bosons $`n_B=b_i^+b_i`$ in this model is not free of ambiguities, given our poor understanding of the doping process in these materials. For a system such as, for instance, YBCO the number of doping induced bound elelctron-pairs could possibly be given by half the number of dopant ions $`O_2^2(1)`$ in the chains, thus varying between $`0`$ and $`0.5`$ per effective site. The number of fermions $`n_F=_\sigma c_{i\sigma }^+c_{i\sigma }`$ is taken to be equal to 1 if the boson-fermion exchange coupling were absent. We thus obtain a total number of charge carriers $`n_{tot}=n_F+2n_B`$ which should be contained between 1 and 2. As a representative example we choose $`n_{tot}=1.25`$. In order to have $`n_F1`$ we fix the bosonic level as $`\mathrm{\Delta }_B=1.1`$. In order to obtain values for $`T^{}`$ of the order of a few hundred degrees K we choose $`v=0.1`$.

The pseudogap manifests itself as a dip which emerges in the DOS of the electrons close to the chemical potential below a certain temperature $`T_F^{}`$. Tracing the value of the DOS at this energy as a function of temperature permits to identify this characteristic temperature. The appearence of this pseudogap is linked, as we have discussed previously, to a breakdown of well defined single-electron excitations close to the Fermi surface. Concomitantly with this trend, two-electron states (local electron-pair resonances) emerge and upon lowering the temperature, acquire itinerant behaviour below a second characteristic temperature $`T_B^{}`$. $`T_B^{}`$ is defined by the condition (see Fig.1) that the imaginary part divided by the real part of the Boson self-energy is small, i.e., $`\gamma _q^B(T)/\omega _q^B(T)\mathrm{\Gamma }_B(q_{},\omega _B(q_{},T),T)/\omega _B(q_{},T)0.1`$ for small $`q`$ vectors such as $`q_{}=0.1`$. Below $`T_B^{}`$ we find $`\omega _B(q_{})=\mathrm{}^2q_{}^2/2m_B(T)`$ for $`q0.2`$) where $`m_B(T)`$ denotes the temperature dependent mass of the itinerant two-particle excitations.

The two-electron states, becoming well defined as the temperature decreases, goes hand in hand with a significant increase of the distribution function for those states with small wave vectors. It is this feature which leads to a finite value of the amplitude of the order parameter in the normal state whilst at the same time phase coherence is totally absent. The temperature $`T_c`$ at which superconductivity sets in is determined by the condition that the two-electron excitations condense in a macroscopic quantum state. According to the Hugenholtz-Pines theorem, it is determined by $`\mathrm{\Delta }_B2\mu Re\mathrm{\Sigma }_B(0,0)=0`$. We plot in Fig.1 these three characteristic temperatures as a function of $`\alpha `$. Notice the opposite trends of $`T_F^{}`$ and $`T_c`$, approaching each other as the isotropic limit is approached, ($`\alpha =10)`$. This is reminiscent of the experimental situation in cuprate HTS materials, as we go from the underdoped towards the optimally doped regime and might hence partly be related to such a dimensional cross-over.

In order to track such cross-over behavior in the pseudogap phase, leading from uncorrelated fluctuating electron pairs to their itinerant behavior as the temperature is decreased, we now show how such features are related to precursor Drude behavior of the optical conductivity and partial Meissner screening above $`T_c`$. We encounter here physics similar to that of the 2D non-interacting charged Bose gas, where a macroscopic number of itinerant bosonic states with momenta $`𝐤1/\xi `$ ($`\xi `$ denoting the coherence length) act together collectively in a way similar to the macroscopically occupied condensed states with $`𝐤=0`$ in the superconducting phase. The optical conductivity $`\sigma (𝐪,\omega )`$ and the diamagnetic susceptiblity $`\chi (𝐪,\omega )`$ are given by the longitudinal, respectively transverse part of the linear response to an external vector potential $`𝐀(𝐱,t)`$ depending on space and time. The resulting current is given by

$`J_i(𝐱,t)={\displaystyle \frac{e^2}{m_{}^2c}}{\displaystyle d^3x𝑑tK_{ij}(𝐱𝐱^{},tt^{})A_j(𝐱^{},t^{})},`$ (7)
$`K_{ij}(𝐱,t)=i\mathrm{\Theta }(t)[j_i(𝐱,t),j_j(0,0)]2m_{}\delta (𝐱)\delta (t)\delta _{ij}n^F(x)`$ (8)
where $`n^F(x)`$ denotes the density of the Fermions and $`j_i^F(𝐪,t)=_𝐤2k_ic_{𝐤𝐪/\mathrm{𝟐}}^{}c_{𝐤+𝐪/\mathrm{𝟐}}`$ the Fourier transform of their current density. Putting

$`K_{ij}(𝐪,\omega )=\delta _{ij}i\omega \sigma (𝐪,\omega )(q_iq_j\delta _{ij}q^2)\chi (𝐪,\omega )`$ (9)
permits to extract the optical conductivity and diamagnetic susceptibility. The kernel $`K_{ij}(𝐱,t)`$ is decomposed into two contributions. A first one is given by the simple bubble for the current autocorrelation function and neglecting vertex corrections. This is justified because of the strong incoherent contributions of $`G_F`$. The second contribution to this kernel is evaluated in terms of the typical Aslamazov-Larkin diagram which takes into account the contributions of the itinerant Bosons. Due to the intrinsically localized nature of the Bosons there are no direct contributions of them to the electrical current. However since below $`T_B^{}`$ those bosonic two-electron excitations become itinerant they will, according to this Aslamazov-Larkin mechanism, contribute to substantially enhance the conductivity. This is manifest in the dc resistivity (see insert of Fig.2) which clearly shows that upon decreasing the anisotropy ratio $`\alpha `$, the resistivity changes qualitatively from an upturn to a downturn as the temperature is lowered. Similarly, the effect of the pseudogap and the precursor to a macroscopic quantum state of the two-particle excitations can be tracked in the optical conductivity (Fig.2). By inspection of Fig. 2 we notice that for frequencies above $`0.01`$ the optical conductivity drops as the temperature is decreased below $`T=0.03`$ (indicative of a remnance of the pseudogap in the DOS of the single-electron states), while for frequencies below $`0.01`$ the optical conductivity for those temperatures increases (indicative of an emerging precursor “Drude” component as we approach $`T_c`$). The crossover temperature at $`T=0.01`$ is identified as the temperature $`T_B^{}`$ where the electron-pairs become itinerant. As the temperature is lowered, spectral weight is shifted downwards from the frequency regime $`[0.01\omega 0.04]`$ thus enhancing this “Drude” component below $`0.01`$. Such a redistribution of spectral weight is in fact observed. For frequencies below 0.001 ($`240GHz`$, taking $`\stackrel{~}{D}=1eV`$) the optical conductivity increases substantially as one approaches $`T_c`$.

Similarly, the real part of the diamagnetic susceptibility $`\chi ^{}(𝐪,\omega )`$ for small wavevectors ($`q1/\xi `$, $`\xi =1/\sqrt{2m_B(T)\omega _0^B}`$ denoting the coherence length) and frequencies is negative and small for high temperatures. It is given by the usual Landau diamagnetism of the electrons arising from the first contribution to the kernel $`K_{ij}`$. For temperatures below $`T_B^{}`$ (which is $`0.01163`$ for a representative example of $`\alpha =10^3`$) we calculate $`\chi ^{}(q,\omega )`$, not taking into account the phase fluctuations. Close to $`T_c`$ $`(=0.01127`$ for $`\alpha =10^3)`$ $`\chi ^{}(q,0)`$ however is expected to diverge as $`1/\xi _{KTB}^2`$, where $`\xi _{KTB}^2`$ is the Kosterlitz-Thouless-Berezinskii coherence length for phase fluctuations. Our calculated $`\chi ^{}(q,\omega )`$ increases significantly (roughly by a factor 50) over the Landau diamagnetism (see Fig. 3). It is given by the Aslamazov-Larkin contribution of $`K_{ij}`$ i.e., $`\chi ^{}(0,0)av^4(\frac{e^2}{24\pi m_B(T)c^2}\frac{k_BT}{\omega _0^B})`$, $`a`$ being a numerical factor of order unity. The expression in the bracket corresponds to that known from the 2D Bose gas. It is this contribution which leads to almost complete Meissner screening. A characteristic frequency scale thus emerges below which phase uncorrelated amplitude fluctuations of the electron-pairs behave in some respects as condensed electron-pairs in the superconducting state. For frequencies $`\omega 10^5`$ we find the usual Landau diamagnetism but for $`\omega 10^5`$ we notice a saturation of $`\chi ^{}(𝐪,\omega )`$ which decreases as the temperature increases. It is in this low frequency regime that we can extract from the diamagnetic susceptibility an almost complete Meissner screening for strong anisotropies. Its physical meaningfulness arises from the fact that for such a system there is phase locking over a given finite distance which is controlled essentially by the thermal wavelength of the itinerant Bosons.

In this Letter we have examined the effect of pure amplitude fluctuations in the pseudogap phase of HTS as the dimensionality of the system varies between quasi-2D and quasi-3D. We conclude that the doping induced cross-over between the under- and the over-doped regime might at least partly be related to such a dimensional cross-over, clearly indicating a shrinking of the temperature regime for the pseudogap phase as the optimally doped limit is approached. We identified two characteristic temperatures in the pseudogap phase. The first one, $`T_F^{}`$, corresponding to the opening of the pseudogap due to the onset of local uncorrelated electron-pair correlations. The second one, $`T_B^{}`$, corresponding to the onset of itinerant behavior of those electron-pairs, below which we expect partial Meissner screening (very similar to that expected for the 2D non-interacting charged Bose gas, i.e., in the absence of any phase fluctuations) and an optical conductivity showing Drude behavior.

We would like to acknowledge helpful discussions with B. K. Chakraverty, T. Domanski, K. Matho and A. Romano. This research was supported in part by a European network contract ERBCHRXCT940438.