text:
# Permeability of self-affine rough fractures

## I Introduction

The transport of fluids through geological media such as hydrocarbon and water reservoirs involves a combination of flow through the microscopic pore space of the rock itself and flows through macroscopic channels such as fractures. The first case is relatively well understood, at least in principle , by means of models which treat the pore space as a random network, and then use effective medium or percolation concepts for the transport. Fractures are typically modeled as simple straight-sided channels, with a cubic relation between fluid flux and average aperture, and a challenging problem is to understand the dynamics of flow in a macroscopic fracture network . Typically the surface of a single fracture appears fairly smooth, aside from some small-scale superficially random roughness, and Poiseuille flow in a straight channel is the obvious model for fluid flow. However, more careful analysis  shows that common geological fractures in fact have correlated, self-affine fractal surfaces. The roughness exponent, whose precise definition is recalled below, is usually found  to be close to 0.8, surprisingly insensitive to the material and the fracturing process. (Other values may arise from intergranular effects .)

The aim of this paper is to study the permeability of single self-affine two-dimensional fractures. We shall see that it is possible to obtain general analytic expressions in two limiting cases, where the roughness associated with the fracture surface is either small or large compared to the mean aperture. The theoretical predictions will be supported by numerical simulations, using the lattice Boltzmann method, which of course can address the case of intermediate sized roughness as well. In subsequent work, we will consider the fully three dimensional case, and go on to consider more general transport processes in fractures, involving both passive tracers and finite-sized suspended particles. Complimentary experimental studies on laboratory samples of fractured rock are in progress elsewhere .

To the extent that the fluctuations in the height are small compared to the aperture width, the effect on permeability is modest, and we will obtain small corrections to the usual cubic law, but the effects on other transport processes are much more significant. For example, passive tracer dispersion is very sensitive to the presence of any slow zones in the velocity field, and even uncorrelated roughness leads to significant long-time tails . Similarly, the deposition and erosion of solid particles from a surface is sensitive to the local shear stress, which in turn can change dramatically as a surface roughens.

In the opposite extreme of a very narrow fracture, the flow field changes qualitatively as the fluid winds through a highly irregular channel. Aside from the challenge of estimating the permeability of difficult geometry in terms of its statistical characteristics, there are further effects arising from the fact that the two sides of the fracture originate from a single crack. In the fractal case, the spatial correlations between the two sides of a fracture lead to velocity field correlations, which again strongly affect tracer motion . The latter references made use of a very approximate velocity field, obtained in the lubrication limit, and one of the motivations of this paper is to examine the validity of the lubrication approximation for the velocity field.

The organization of this paper is as follows. We first recall some basic facts about self-affine surfaces, and give the algorithm used to generate them numerically, and also explain the lattice Boltzmann method used for the numerical simulations. We then consider the case of fractures with a wide average aperture, obtain a perturbative estimate for the permeability, and test it numerically. We then turn to narrow fractures, first in the case where the two sides are simply displaced normally to the mean fracture plane, and secondly when there is a lateral shift as well. Finally, we summarize the results, and indicate the next steps

## II Preliminaries

### A Self-affine roughness

We briefly review the mathematical characterization of self-affinity, and its implementation in this paper. We consider a rock surface without overhangs, whose height is given by a by a single-valued function $`z(x,y)`$, where the coordinates $`x`$ and $`y`$ lie in the mean plane of the fracture. Self-affine surfaces  display scale invariance with different dilation ratios along different spatial directions, remaining unchanged under the rescaling

$$x\lambda _1xy\lambda _2yz\lambda _3z$$
(1)
Here we consider disordered media, so these scaling laws apply only in an ensemble or spatial average sense. Experiment indicates that for many materials isotropy can be assumed in the mean plane, implying that there is only one non-trivial exponent relating the dilation ratio in the mean plane to the scaling in the perpendicular direction, i.e., $`\lambda _1=\lambda _2\lambda `$, and $`\lambda _3=\lambda ^\zeta `$, with

$$z(x,y)=\lambda ^\zeta z(\lambda x,\lambda y)$$
(2)
where $`\zeta `$ is the roughness or Hurst exponent .

We assume that the process of fracturing the rock is “clean”, in the sense that the rock breaks so as separate along one single valued surface without subsequent deformation and without producing loose interstitial material. We shall emphasize two complimentary limiting situations, in which the two surfaces are either well separated with respect to each other, or very close to each other. If $`L`$ is the lateral size of the fracture in the mean plane, the typical range of the fluctuations (the difference between the maximum and minimum values of $`z`$) scales as $`RL^\zeta `$. If $`h`$ is the mean width of the fracture, with $`hL`$, then the two cases correspond to $`Rh`$ and $`Rh`$, respectively.

The local aperture of the crack is the key parameter determining the fracture permeability, and three different cases will be addressed. For large aperture, the roughness associated with the two sides of the fracture act independently, and it suffices to consider a channel with one rough and one smooth boundary. In the case of a narrow fracture, the correlation between the two sides is a key feature, and we consider two possibilities. We first study fractures where the upper surface has been simply translated a distance $`h`$ normal to the mean plane, so that the local aperture $`a(x,y)`$ equals the constant $`h`$. Alternatively, the two surfaces may have a relative lateral displacement $`d`$ in the mean fracture plane, accompanied by a displacement $`h`$ in the perpendicular direction, so that the two surfaces do not overlap. In this case the local aperture is given by the random variable.

$$a_d(x,y)=z(x+d,y)z(x,y)+h$$
(3)
It turns out  that $`d`$ is the lateral correlation length for fluctuations in the aperture, in the sense that $`a_d(x,y)`$ and $`a_d(x+\mathrm{\Delta }x,y)`$ decorrelate for $`\mathrm{\Delta }xd`$.

In this paper we restrict ourselves to the two dimensional case where the surface is invariant in the $`y`$ direction, $`z(x,y)=z(x)`$, and the flow is forced in the $`x`$ direction by a constant pressure drop. In a subsequent paper we will extend these calculations to fully three dimensional fractures, but it is convenient, both conceptually and in the numerical simulations, to regard the system as having a translationally invariant third dimension.

### B Numerical Methods

Our aim is to study various aspects of the fracture permeability which are sensitive to the fracture roughness. The first ingredient is the height function $`z(x,y)`$, a statistically self-affine surface with periodic boundary conditions. The periodicity is not a physically essential ingredient here, but has some calculational advantages in alleviating finite-size effects. The surface is generated by a Fourier synthesis method, based on power-law filtering of arrays of independent random numbers . The random numbers are generated using a Gaussian distribution, and then modulated by an appropriate power law. If $`\stackrel{~}{Z}(k)`$ is the Fourier transform of the initial Gaussian random array, then the Fourier transform of the surface elevation is chosen to be

$$Z(k)=k^{\zeta 1/2}\stackrel{~}{Z}(k)$$
(4)
(The 1/2 is appropriate for the self-affine curve used here, and would be 1 for a real rock surface.)

In most cases the roughness exponent is chosen as the often-observed value $`\zeta =0.8`$. The amplitude of the roughness can be expressed in terms of variance of the height distribution over the full range,

$$\sigma ^2=(z(x)z(x))^2=\frac{1}{L}_0^L(z(x)z(x))^2𝑑x,$$
(5)
which in this case is $`\sigma ^2=0.49\pm 0.06`$ for $`L=256`$. Alternatively, the amplitude is given indirectly in terms of the range of the fluctuations as $`R=R_1L^\zeta `$, with $`R_1=1.66\pm 0.4`$. In these equations, and in the rest of the paper, the unit of length is that of the spacing between lattice points in the numerical calculations.

Since we consider fluid flow in a highly irregular geometry, and eventually hope to consider dispersion and the evolution of the solid surface due to particle transport, the lattice Boltzmann (LB) method  is particularly convenient. In this algorithm, fictitious particles move between neighboring sites on a lattice, with suitable collision rules, and the boundary of the flow domain is simply a surface of sites where the rule is modified in some way to keep the particles out. We use the Face Centered Hypercubic (FCHC-projected) version of the LB model, with a cubic lattice in 3 dimensions and 19 velocities (D3Q19 in the terminology of ). The collision operator is approximated by a single relaxation time, the Bhatnagar-Gross-Krook model, and the local equilibrium distribution given by  is used. This pseudo-equilibrium distribution locally preserves mass and momentum values, and is formulated specifically to recover the Navier-Stokes equation at large length and time scales. For the no-slip solid boundaries, we use the simplest implementation of particle exclusion – the bounce-back rule, in which a particle incident on the boundary reverses its direction. Periodic boundary conditions are used for the inflow and outflow surfaces. A constant pressure gradient forcing the fluid is added in the $`x`$ direction, while the gap between surfaces extends over $`z`$, and the geometry is taken to be translation invariant in $`y`$.

## III Small surface deviations

In this section we consider the case of a channel with one fractally-rough wall, in the limit where the mean width is large compared to the amplitude of the roughness. We (re-)derive an elegant general result for the permeability, and then an exact variant which gives the permeability of a perturbed pore space. If the pore space perturbation is weak, the leading order correction is easy to evaluate and, when compared to numerical simulation, is seen to provide a reasonably accurate estimate of the permeability.

### A Integral representations

We begin by deriving an integral representation for the permeability, and then considering a perturbation in the boundary shape. Although we suspect that these results are known to many, we are not aware of an earlier publication containing them, although the electrical analog of the integral representation is in the literature .

Consider a rectangular volume of a porous medium, periodic in all directions, and the following surface integral

$$I=_S𝑑\stackrel{}{s}\left[\stackrel{}{u}(\stackrel{}{r})p(\stackrel{}{r})\right]$$
(6)
Here, $`\stackrel{}{u}`$ and $`p`$ are the velocity and pressure fields of a fluid which completely fills the pore space, assumed to satisfy the Stokes equations, and the surface $`S`$ is the complete boundary of the pore space, consisting of the inner grain surfaces and the porous regions of the outer boundary. The pressure and velocity fields are periodic, except that the pressure jumps by a constant amount $`P`$ between two opposite faces in one direction, that of the mean flow. Now the velocity vanishes on the grain boundaries by the no-slip condition, and there is complete cancelation between the opposite faces in the two completely periodic directions, but in the flow direction the remaining two faces of the box combine to give

$$I=P_E𝑑\stackrel{}{s}\stackrel{}{u}(\stackrel{}{r})=PQ=\frac{kP^2S}{\mu L}$$
(7)
Here $`Q`$ is the fluid flux through this end face $`E`$, and the minus sign arises because the outward normal is used in $`I`$. The last equality follows from Darcy’s law, where $`S`$ is the area of end $`E`$, $`k`$ is the permeability, $`\mu `$ the fluid viscosity, and $`L`$ the length of the box.

On the other hand, applying the divergence theorem to Eq. 6, we obtain:

$$I=_V𝑑V(\stackrel{}{u}p)=_V𝑑V(\stackrel{}{u}p)=\mu _V𝑑V(\stackrel{}{u}^2\stackrel{}{u})=\mu _V𝑑V(\stackrel{}{u})^2$$
(8)
where the volume $`V`$ is just the pore space. We have used first the incompressibility of the fluid, second the Stokes equation, and third integration by parts. Comparing the two expressions for $`I`$, we have

$$k=\frac{L}{S}\left(\frac{\mu }{P}\right)^2_V𝑑V\underset{i,j}{}(_iu_j)^2$$
(9)
The analogous formula for the electrical conductivity $`\sigma `$, which follows from Ohm’s law and the Laplace equation , is

$$\sigma =\sigma _0\frac{L}{S}\frac{1}{\mathrm{\Phi }^2}_V𝑑V\underset{i}{}(_i\varphi )^2$$
(10)
Here $`\sigma _0`$ is the conductivity of the pore fluid, $`\varphi (\stackrel{}{r})`$ is the potential and $`\mathrm{\Phi }`$ its difference across the sample. Note that both of these formulae may also be derived by comparing the energy dissipation rate computed microscopically in the pore space to the same quantity computed using average fields.

We next derive an exact variant of the integral expression for the permeability (9) due to Wilkinson , which allows us to implement a perturbation analysis. Suppose we begin with a pore volume $`V_0`$ and known Stokes equation solutions $`p_0`$ and $`\stackrel{}{u}_0`$, and then contract the volume to $`V`$, where the exact (but unknown) solution is $`p=p_0+p_ϵ`$ and $`\stackrel{}{u}=\stackrel{}{u}_0+\stackrel{}{u}_ϵ`$. The result to follow is true even if the volume change and alterations in the fields are not small, but are probably only useful in that limit, hence the suggestive subscript $`ϵ`$. Substituting into the dissipation integral on the right hand of Eq. 9, and suppressing the summation sign for the moment, we obtain

$$_V𝑑V(_iu_j)^2=_V𝑑V\left[(_iu_{0,j})^2+(_iu_{ϵ,j})^2\right]+2_V𝑑V_iu_{0,j}_iu_{ϵ,j}$$
(11)
The last integral on the right hand can be rewritten as

$$_V𝑑V\left[_i(u_{0,j}_iu_{ϵ,j})u_{0,j}_i^2u_{ϵ,j}\right]=_S𝑑S_iu_{0,j}_iu_{ϵ,j}_V𝑑Vu_{0,j}_i^2u_{ϵ,j}.$$
(12)
Noting that $`u_{0,j}=u_{ϵ,j}`$ at the surface $`S`$ due to the no-slip condition, and using $`\mu _i^2u_{ϵ,j}=_jp_ϵ`$, we may rewrite the previous expression as

$$_V𝑑V_iu_{0,j}_iu_{ϵ,j}=_V𝑑V\left[(_iu_{ϵ,j})^2+\mu ^1\stackrel{}{u}p_ϵ\right]$$
(13)
Using incompressibility, the last term in the integrand can be rewritten $`(\stackrel{}{u}p_ϵ)`$ and converted to a surface integral, which vanishes if the pressure is held constant on the ends of the porous medium. Finally, we obtain

$$k=\frac{L}{S}\left(\frac{\mu }{P}\right)^2\underset{i,j}{}_V𝑑V\left[(_iu_{0,j})^2(_iu_{ϵ,j})^2\right]$$
(14)
This result is exact but not useful. However, note that if the relative change in volume is small, $`O(ϵ)`$, the second term is second order in the small parameter, so that to first order the decrease in permeability involves only the integral of the unperturbed velocity over the deleted region, which is easy to calculate.

### B Dependence of permeability on system size

In this section we will determine the decrease in permeability due to the presence of a self-affine perturbation in the lower surface of a straight channel, and in particular its dependence on the size of the channel. To use the previous result, we begin with a straight channel of width $`h_0`$ which completely contains the rough-walled one (see Fig. 1), where we have Poiseuille flow

$$\stackrel{}{u}_0=\frac{1}{2\mu }\frac{P}{L}z(h_0z)\widehat{x}$$
(15)
and an unperturbed permeability $`k_0=h_0^2/12`$. The lower boundary is shifted by $`z(x)>0`$, where we assume that that range $`R`$ of $`z`$ satisfies $`ϵR/h_01`$. The only surviving velocity derivative is

$$_zu_{0,x}=\frac{1}{2\mu }\frac{P}{L}(h_02z)$$
(16)
so that the dissipation integral is simply

$`{\displaystyle _V}𝑑V{\displaystyle \underset{i,j}{}}(_iu_{0,j})^2`$ $`=`$ $`{\displaystyle _0^L}𝑑x{\displaystyle _0^Y}𝑑y{\displaystyle _{z(x)}^{h_0}}𝑑z\left(_zu_{0,x}\right)^2`$ (17)
$`=`$ $`\left({\displaystyle \frac{P}{\mu L}}\right)^2Y{\displaystyle _0^L}𝑑x\left\{{\displaystyle _0^{h_0}}𝑑z{\displaystyle _0^{z(x)}}𝑑z\right\}h_0^2\left({\displaystyle \frac{z}{h_0}}{\displaystyle \frac{1}{2}}\right)^2`$ (18)
Up to first order in $`ϵ`$ this becomes

$`{\displaystyle _V}𝑑V(_iu_{0,j})^2`$ $``$ $`\left({\displaystyle \frac{P}{\mu L}}\right)^2Y\left\{{\displaystyle \frac{h_0^3}{12}}L{\displaystyle \frac{h_0^2}{4}}{\displaystyle _0^L}z(x)𝑑x\right\}`$ (19)
$``$ $`\left({\displaystyle \frac{P}{\mu }}\right)^2{\displaystyle \frac{S}{L}}k_0\left\{1{\displaystyle \frac{3}{h_0}}z\right\}`$ (20)
where $`z`$ is the average of $`z`$ over the channel length $`L`$, and the cross-sectional area is $`S=h_0Y`$. The mean perturbation $`z`$ is computable for a specific profile, but here we wish to use a statistical characterization and relate it to the average properties of an ensemble of self-affine surfaces.

The mean height of the surface is half the range, and therefore has the scaling $`z(x)=\frac{1}{2}R_1L^\zeta `$. Substituting, and replacing the previous equation in Eq. 9 we obtain the first-order perturbative correction to the permeability

$$kk_0\left[1\frac{3}{2}C_1ϵ\right]$$
(21)
We have added an adjustable parameter $`C_1`$, which is expected to lie between $`1<C_1<2`$, to take account of the distinction between open volume and flowing volume. A value $`C_1=1`$ means that the fracture is equivalent to a straight channel of height equal to the mean height of the surface. On the other hand, a value $`C_1=2`$ means that the whole region of fluid below the maximum excursion of the surface is not contributing to the permeability and the system is equivalent to the maximum channel that does not intersect the surface.

Note the lowest order correction to the permeability obtained here results from the reduced volume of a fracture compared to that of a straight channel enclosing it, and happens to coincide with the lubrication approximation. The approach taken here allows us to gain some insight in the drawbacks of this approximation and further improve them. Using the previous equation we may write the first order correction to the flow rate,

$$Q(L)Q_0\left[1\frac{3C_1R_1}{2h_0}L^\zeta \right]$$
(22)
To verify this prediction, we consider a pore space consisting of a cubic lattice with periodic boundaries in $`x`$ and $`y`$, one impermeable wall at $`z=h_0`$, and a self-affine rough boundary lying above $`z=0`$. We generate an ensemble of 10 self-affine surfaces with exponent $`\zeta =0.8`$, as discussed above. In Fig. 2 we display the correction to the flow rate as a function of the size of the system $`L`$. The straight line is a fit to the expected exponent, which does reasonably well except at the smallest and largest values of $`x`$. The origins of the discrepancies are first that the discretization used in the numerical calculation suppresses any roughness less than one unit, which is a significant fraction of the system for small $`L`$. Moreover, the system should be periodic in $`L`$, but the numerical algorithm used to generate it assumes a system size that is a power of 2. As usual, there are computational limitations on the sizes we can simulate, and for lengths between 32 and the accessible maximum of 64, we just truncate the periodic surface. Finally, at the largest system sizes, we are really outside the range of validity of the estimates. In fact, the second term in the flux in Eq. 22 is greater than 1/3 for the parameters here when $`L16`$.

The fitted value $`C_1`$ given that $`\zeta =0.8`$ is $`C_1=1.30\stackrel{+2.1}{1.2}`$, meaning that a portion of the fluid does not contributes to convective transport. Expressed in terms of the surface height variance $`\sigma ^2`$ instead of the span range, the excluded fluid represents $`50\%`$ of the roughness region. The physical consideration absent in the leading order perturbation calculation is the fact that fluid velocity sharply decays inside depressions and corners were resistive eddies are likely to occur . This deficiency is associated with the omitted second term in Eq. 14, i.e., the correction to the unperturbed velocity. In Fig. 1 we magnify two regions close to the lower and upper surfaces, respectively, showing the markedly different way the velocity decays towards the surface.

We have also simulated a system with an alternative roughness exponent $`\zeta =0.95`$, corresponding to a surface with more persistent correlations and in a sense a lower fractal dimension (the dimension of the intersection of the surface with a plane normal to it is $`D=2\zeta `$ ). In this case the prefactor in the range scaling law is $`R_1=0.6\pm 0.12`$. In Fig. 2 we show the numerical results; again, good agreement with the predicted behavior is obtained. The fit gives $`C_1=2.04\stackrel{+3.1}{1.8}`$, suggesting that the low-velocity zones are even more important in this case, presumably because the greater degree of correlation enhances them.

### C Permeability decrease due to zero mean surface fluctuations

We wish to disentangle the permeability decrease due to volume reduction of the fracture from that due to low-velocity zones induced by surface fluctuations. The first step is to consider “zero-mean” self-affine surfaces, whose mean height $`h`$ is constant, but which have a tunable roughness. The previous profile, $`z(x)`$ for $`0<x<L`$, is rescaled so as to have unit variance, and then multiplied by an adjustable amplitude $`A`$, giving adjustable variance $`A^2`$. If we assume that the effect of the boundary perturbation is to give low velocity zones which do not contribute to the flow, and that the volume of these zones contributing to the dissipation integral is proportional to $`A`$, we would write in analogy to Eq. 22

$$Q=Q_0\left[1\frac{3C_2A}{h}\right]$$
(23)
Here $`A`$ is analogous to $`R_1L^\zeta /2`$, and $`C_2`$ an order-1 fitting coefficient. In Fig. 3 we show the results obtained by numerical simulation corresponding to a system with $`h=32`$ and $`L=64`$ (and, to be specific, $`p=6.25\times 10^6`$, $`\mu =0.1`$, $`Y=4`$, which gives $`Q_0=0.683`$). The filled circles are the numerical results, and we see the expected linear decrease of the permeability with $`A`$. In fact, a linear fit to the numerical data gives $`Q_0=0.72\pm 0.03`$ and $`C_2=0.77\pm 0.1`$, which are reasonable numerical values given our assumptions.

Further evidence for our arguments concerning the effects of low-velocity zones can be adduced by showing that the permeability is unchanged if they are deleted. We begin with the family of adjustable-amplitude self-affine surface just discussed, and filter out the smallest wavelengths from the spectrum so as to generate a smother surface. The smooth surface is then shifted upwards so as to enclose the deleted surface fluctuations, in the RMS sense. Quantitatively, if we have $`N_L`$ points on the $`x`$-axis, with $`L=N_L\mathrm{\Delta }`$, and filter all wavelengths smaller than $`\lambda _c=\mathrm{\Delta }N_L/N_C`$, then

$$(\mathrm{\Delta }h)^2=\frac{1}{(N_LN_C)^2}\underset{n=N_C+1}{\overset{N_L}{}}|Z(k_n)|^2k_n=\frac{2\pi n}{\mathrm{\Delta }N_L}$$
(24)
The prefactor on the right arises from the discrete form of Parseval’s theorem , and in practice we choose $`\mathrm{\Delta }=1`$. We then simulate fluid flow in the smoothed fracture, and in Fig. 3 the open squares show the results obtained as a function of roughness amplitude $`A`$. The agreement between the actual and filtered surfaces is clear.

The extension of this analysis to the 3-D case is conceptually straightforward, since the decrease in the permeability in this limit just corresponds to the decrease in pore space volume. The latter would be $`L_x\times L_y\times R`$, due to the self-affine topography of the fracture surface, and the relevant expansion parameter is again $`R/h_0`$. The effects of low-velocity regions is presumably the same as well.

## IV Narrow fractures

We now turn to the complimentary situation of narrow fractures, where the pore space is winding and convoluted. Consider the situation in which a rock is fractured and the two surfaces are simply displaced by a distance $`hL`$, perpendicular to the mean fracture plane, with no lateral shift. The vertical aperture is constant everywhere, but the effective local aperture for fluid flow, i.e., the local width of the channel normal to the mean flow direction, strongly depends on the local angle between the surface and the mean plane. In Fig. 4 we show a fracture of length $`L=64`$, generated by a self-affine surface of exponent $`\zeta =0.8`$, separated a distance $`h=8`$, along with the (lattice Boltzmann) computed velocity field.

The theoretical analysis will be based on a kind of local lubrication approximation, where the channel is divided into a sequence of quasi-linear segments at varying orientation angles. First, we estimate the length $`\xi _{}`$ in the direction of the mean flow over which the channel formed by the two fracture surfaces can be considered approximately straight, which is a typical size over which the fluctuations in the vertical direction are small compared to the effective aperture of the channel. Using the self-affine scaling law for the correlation function,

$$\sigma _z^2(\xi _{})=[z(x,y)z(x+\xi _{},y)]^2=\varphi (\mathrm{})\left(\frac{\xi _{}}{\mathrm{}}\right)^{2\zeta }$$
(25)
where $`\mathrm{}`$ is a microscopic length, say a grain size, and $`\varphi (\mathrm{})`$ is then on the order of a grain size squared, we estimate

$$z(x,y)z(x+\xi _{},y)\mathrm{}\left(\frac{\xi _{}}{\mathrm{}}\right)^\zeta $$
(26)
A segment of length $`\xi _{}`$ is roughly straight when this quantity is a small fraction of $`h`$, which yields

$$\xi _{}\mathrm{}\left(\frac{h}{\mathrm{}}\right)^{1/\zeta }$$
(27)
Returning to the entire fracture, each $`\xi _{}`$-channel is oriented at some angle $`\theta _i`$ with respect to the mean plane, and has effective aperture $`a_i=h\mathrm{cos}\theta _i`$, and length $`\xi _{}^i=\xi _{}/\mathrm{cos}\theta _i`$. Assuming Poiseuille flow in each $`\xi _{}`$-channel, the corresponding local permeability is $`k_i=a_i^2/12`$. We compute the pressure drop across the fracture by adding the the local pressure drops along the straight $`\xi _{}`$-channels along the whole path. Noting that the flow rate $`Q=\overline{U}h=\overline{U}_ia_i`$ is constant,

$`\mathrm{\Delta }P`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{\Delta }P_i={\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{\mu }{k_i}}\overline{U}_i\right)\xi _{}^i=`$ (28)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{12\mu }{h^2\mathrm{cos}^2\theta _i}}\right)\left({\displaystyle \frac{Q}{h\mathrm{cos}\theta _i}}\right)\left({\displaystyle \frac{\xi _{}}{\mathrm{cos}\theta _i}}\right)`$ (29)
$`=`$ $`{\displaystyle \frac{12\mu \xi _{}Q}{h^3}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{cos}^4(\theta _i)={\displaystyle \frac{12\mu \xi _{}\overline{U}}{h^2}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{cos}^4(\theta _i)`$ (30)
Finally, noting that $`N=L/\xi _{}1`$ is the number of channels, we can convert the sum into an average over the distribution of angles, and write for the permeability

$$k=\frac{h^2}{12}\frac{1}{\mathrm{cos}^4(\theta )}$$
(31)
We can give a simple if heuristic estimate of the permeability as follows. Since the exponent $`\zeta <1`$, the channels have small vertical fluctuations, and we can approximate the cosine as

$$\mathrm{cos}\theta =\frac{\xi _{}}{\sqrt{\xi _{}^2+\sigma _z^2(\xi _{})}}1\frac{1}{2}\left(\frac{\sigma _z(\xi _{})}{\xi _{}}\right)^2$$
(32)
Substituting in Eq. 31, we obtain

$$k\frac{h^2}{12}\left[12\left(\frac{\sigma _z(\xi _{})}{\xi _{}}\right)^2\right]$$
(33)
A more convincing evaluation of the angular average makes use of the actual height distribution. Experimental measurements indicate a Gaussian distribution for the spatial correlation in heights, $`Zz(x,y)z(x+\xi _{})`$ , and in fact our numerical procedure for generating self-affine surfaces also gives a Gaussian distribution for $`Z`$. We illustrate this point in Fig. 5, where we plot our generated probability distribution function of $`Z`$, corresponding to different values of $`\xi _{}`$, along with their Gaussian fits

$$p(Z)=\frac{1}{\sqrt{2\pi \sigma _z^2(\xi _{})}}\mathrm{exp}\left(\frac{Z^2}{2\sigma _z^2(\xi _{})}\right)$$
(34)
We also plot the (numerically-obtained) dependence of the mean spatial correlation, $`Z^2=\sigma _z^2(\xi _{})`$, on $`\xi _{}`$ for two choices of $`\zeta `$, confirming the scaling given above in Eq. 25.

The angular average is then given by

$`\mathrm{cos}^4\theta `$ $`=`$ $`{\displaystyle 𝑑Zp(Z)\left[1+2\left(\frac{Z}{\xi _{}}\right)^2+\left(\frac{Z}{\xi _{}}\right)^4\right]}`$ (35)
$`=`$ $`1+2\left({\displaystyle \frac{\sigma _z(\xi _{})}{\xi _{}}}\right)^2+3\left({\displaystyle \frac{\sigma _z(\xi _{})}{\xi _{}}}\right)^4`$ (36)
Using Eqs. 25,31, this result is consistent with the previous estimate Eq. 33 based on simple scaling arguments.

Moreover, aside from the purely numerical coefficients, we can argue that scale invariance requires the angular average to have the form given in (35). The self-affine property of the fracture geometry implies that the probability distribution of heights $`p(Z)`$ should only depend on the rescaled variable $`Z/\xi _{}^\zeta `$ . More specifically,

$$p(Z)=\sigma _z^1(\xi _{})\psi [Z/\sigma _z(\xi _{})]$$
(37)
where the prefactor comes from the normalization. If $`\mathrm{cos}^4\theta `$ is evaluated for this distribution, we obtain a variant of Eq. 35 where the numerical coefficients 2 and 3 are replaced by the moments of the function $`\psi `$, leading to the same scaling form for the permeability.

Finally then, using the leading term in Eq. 35 for the angular average, and Eq. 27 for the value of $`\xi _{}`$, we have the result for the permeability of a narrow two dimensional self-affine fracture,

$`k={\displaystyle \frac{h^2}{12}}\left[12\left({\displaystyle \frac{\varphi (\mathrm{})}{\mathrm{}^2}}\right)\left({\displaystyle \frac{h}{\mathrm{}}}\right)^{\frac{2\zeta 2}{\zeta }}\right]`$ (38)
The principal approximation used in obtaining this result is that the fracture aperture may be regarded as a sequence of almost-straight segments.

To test Eq. 38 numerically, we first recast it in terms of the fluid flux. For a straight channel of height $`h`$, width $`Y`$ and length $`L`$, and pressure drop $`P`$, we have flux $`Q_0=h^3PY/12\mu L`$. Thus

$$Q_0QC_3\frac{PY}{6\mu }\mathrm{}^{\frac{22\zeta }{\zeta }}h^{\frac{5\zeta 2}{\zeta }}$$
(39)
where $`C_3`$ is another adjustable parameter expected to be of order one, and we have taken $`\mathrm{}`$ to be the (unit) lattice spacing.

In Fig. 6 we present the correction to $`Q_0`$ as a function of the distance between the opposite surfaces $`h`$, for both roughness exponents $`\zeta =0.80`$ and $`\zeta =0.95`$. We find the good agreement for the predicted exponent for the case $`\zeta =0.80`$. The only adjustable parameter is the coefficient, which is found to be $`C_3=3.1\pm 1.2`$ in fairly good agreement with the expected value. The fitted exponent is $`2.43\pm 0.08`$, which agrees with the predicted value $`(5\zeta 2)/\zeta =2.5`$. On the other hand, in the case $`\zeta =0.95`$ the numerical results deviates from the asymptotic behavior at small distances between the surfaces. That is because the successive powers given by the expansion of Eq. 38 are very close to each other (the difference between them is $`(2\zeta 2)/\zeta 0.1`$). Again if only the coefficient is fitted, we find $`C_3=0.35\pm 0.21`$. Fitting also the exponent, but using only the numerical simulations with large separation between surfaces, we obtain $`2.71\pm 0.39`$ where $`(5\zeta 2)/\zeta =2.89`$.

In Fig. 7 we compare numerical results corresponding to systems of different sizes with the same roughness exponent $`\zeta =0.80`$. It can be seen that the accuracy in the exponent is improved as the size of the system is larger. We obtain a value $`2.52\pm 0.04`$ for the largest system with size $`L=256`$.

As mentioned, the crossover between small and large surface roughness, (compared to the mean aperture of the fracture), can be addressed numerically. The narrow fracture regime, where $`Q_0Qh^{\frac{5\zeta 2}{\zeta }}`$, is expected to be valid when $`\xi _{}L`$. The critical mean aperture value $`h_c`$ which corresponds to $`\xi _{}L`$ turns out to be, using the fitted value of $`C_3`$, $`h_c100`$. On the other hand, the small surface-deviation scaling behavior can be obtained from Eq. (22), where $`L`$ is the constant size of the system and $`h_0`$ is now the variable mean aperture of the fracture $`h`$,

$$Q_0Q\frac{PY}{12\mu }\frac{3C_1R_1L^\zeta }{2}h^2$$
(40)
In Fig. 8 we display the correction to $`Q_0`$ in the range $`2h128`$ for a system size $`L=64`$. A small deviation from the narrow fractures regime can be observed, starting at $`h=64`$. As the size of the gap becomes larger, the correction to the unperturbed flow rate $`Q_0`$ decreases as a power of $`h`$. Therefore, both larger computational time (starting the simulation from an initially at rest fluid the time interval to reach the steady asymptotic value grows as the square of the gap size ) and higher accuracy in the fluid field close to the solid surface are needed. Our current computational limitations prevents us to simulate systems where $`L\xi _{}`$. Nevertheless, we believe that the observed deviation correspond to a crossover towards a smaller exponent ($`h^2`$) as expected from the previous discussion.

## V Fractures with shifted surfaces

In the field, when a rock is fractured its two sides may be shifted laterally parallel to the mean plane by geological processes. We now consider how our results for permeability are modified by this effect. The situation is somewhat different in two dimensions than in three. In order that a two dimensional fracture remain open it is necessary that the two sides do not touch, whereas in three dimensions in the presence of a lateral shift it is easy for the two surfaces to touch at one or a few points in the plane, while the aperture is still open to flow. Furthermore, in two dimensions when the fracture nearly touches at a point, the permeability will be dominated by the large pressure drop needed to force fluid through this narrow gap. In fact, in this case an estimate of the permeability is simply $`ka_{\mathrm{min}}^2/12`$, where $`a_{\mathrm{min}}`$ is the minimum value of the aperture. In this section we address the complimentary case where the fracture is distinctly open.

As discussed in Section 2, in the presence of a lateral shift $`d`$, the fracture aperture is now the random function $`a_d(x,y)`$ given by Eq. 3 rather than the constant $`h`$. However, there is a residual correlation between the two sides of the fracture, which allows us to relate the permeability to the self-affine statistics. First, we determine a condition on the shift for the fracture to be open. Provided the average properties of the surface are statistically stationary (independent of $`x`$), we have $`a_d(x)=h`$, and the variance of the aperture distribution is

$$\sigma _a^2(d)[a_d(x)a_d(x)]^2=[z(x+d)z(x)]^2=\sigma _z^2(d)=\varphi (\mathrm{})\left(\frac{d}{\mathrm{}}\right)^{2\zeta }$$
(41)
The aperture is surely open when $`\sigma _a(d)h`$, which gives $`d\mathrm{}(h/\mathrm{})^{1/\zeta }`$, but from Eq. 27, this is just $`d\xi _{}`$. In other words, the shift must be small compared to the typical length of a straight segment of unshifted channel, and there is little change in the geometry compared to the unshifted case. We expect, therefore, that the previous result Eq. 39 should apply at least up to a different value of the numerical coefficient $`C_3`$.

In Fig. 9 we show LB simulation results for the correction to the unperturbed flow rate $`Q_0`$, for different values of the shift $`d`$. We see that the behavior is consistent with these arguments for all values of the shift, with deviations occurring when the shift is too large for the mean aperture. Several other shift values were simulated obtaining the same general behaviors as in those shown in Fig. 9. As is often the case with this type of scaling argument, reasoning in terms of asymptotic limits leads to conditions that one quantity must be much larger than another, but in practice the range of validity is much wider.

## VI Conclusions

We have studied the permeability of two-dimensional self-affine fractures, using asymptotic analytic arguments for wide and narrow apertures. Numerical simulations using the lattice Boltzmann method have verified the predictions, and also suggest a smooth crossover between the limits used in the derivations. We have obtained expressions for the permeability in which the usual expression for straight channels is modified by terms related to the Hurst exponent characterizing the fracture surface.

In the wide-gap case, we obtained a perturbative correction related to the roughness exponent and amplitude. The result is formally identical to the lubrication theory prediction, but the corrections are known, in principle, and furthermore we can understand the discrepancies between the derived numerical coefficient and simulation in terms of the effects of low-velocity zones which do not contribute to transport. In the narrow gap case, we also derive an analytic correction related to the roughness exponent and amplitude (or microscopic length $`\mathrm{}`$). The analysis is based in a local lubrication approximation, dividing the whole fracture into a chain of approximately linear channels with varying orientation angles respect to the main flow direction. The results is equivalent to a correction due to the tortuosity of the fracture, where the tortuosity itself depends on the gap size. Finally, we show that in the case when there is a lateral shift between the surfaces, there is no qualitative change compared to the unshifted case. This is basically due to the two dimensional geometry, where in order to have an open channel the shift must be small compared to the size of the typical length of a quasi-linear segment of unshifted fracture.

The two-dimensional nature of the problem imposes restrictions on the range of validity of analytic estimates, but the numerical results are in general agreement over a much wider range.

The extension of this work to fully three-dimensional fractures is in progress. The case of a wide gap does not require any further significant conceptual effort. However, the narrow gap case the fracture geometry is significantly more complicated because it is not feasible to use any quasi-one-dimensional approximation for the flow paths. Numerical simulations are certainly feasible, at least up to size limits imposed by one’s computational hardware, but additional ideas are required for analytic arguments. Future work will consider diffusion, and also convection plus diffusion, in these self-affine fractures, in both two and three dimensions. The LB method is quite simple for numerical simulations, and we will explore the way in which the roughness exponent enters into the quantitative behavior, along the lines given here.

## Acknowledgments

We thank J. P. Hulin, F. Plouraboué and M. Tanksley for discussions. This research was supported by the Geosciences Research Program, Office of Basic Energy Sciences, U.S. Department of Energy, and computational facilities were provided by the National Energy Resources Scientific Computer Center.


format: text

text:
# 1 Introduction

## 1 Introduction

####

Plasma and Astrophysical Plasma physicists support the possible existence of electromagnetic stationary waves with parallel $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$ fields and consequently having a null Poynting vector . K.R. Brownstein  points out that these waves may emerge as solutions of the vector equation $`\stackrel{}{}\times \stackrel{}{V}=k\stackrel{}{V}`$ ($`\stackrel{}{V}`$ is a vector field and $`k`$ is a positive constant). Brownstein considers this equation for the vector potential $`\stackrel{}{A}`$ as

$$\stackrel{}{}\times \stackrel{}{A}=k\stackrel{}{A},$$
(1.1)
with a particular solution

$$\stackrel{}{A}=a\left[i\mathrm{sin}kz+\widehat{j}\mathrm{cos}kz\right]\mathrm{cos}\omega t$$
(1.2)
and takes the associated electric and magnetic field as

$$\stackrel{}{E}=\frac{1}{c}\frac{\stackrel{}{A}}{t}=ka\left[i\mathrm{sin}kz+\widehat{j}\mathrm{cos}ka\right]\mathrm{sin}\omega t$$
(1.3)
and

$$\stackrel{}{B}=\stackrel{}{}\times \stackrel{}{A}=ka\left[\widehat{i}sinkz+j\mathrm{cos}kz\right]\mathrm{cos}\omega t.$$
(1.4)
The Brownstein fields (1.3) and (1.4) satisfy Maxwell equations and the usual vacuum free-wave equation. Moreover, fields, $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$ and the vector potential, $`\stackrel{}{A}`$, are parallel everywhere. The associated electromagnetic wave has a null Poynting vector and, so, does not propagate energy. The behavior of this wave is just like the phenomenology suggested by Plasma and Astrophysical Plasma Physics.

But it is pointed out  that this model is not complete, at least for Astrophysical Plasma. There it is argued that in an Astrophysical Plasma, gravitation must be taken into account in a way that the gravitational background can break the parellelism between the fields $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$. Consequently, the corresponding eletromagnetic wave does not have a null Poynting vector. It is not a stationary wave and propagates energy.

In this work, we consider the electromagnetic and gravitional coupling and we analyze, in a gravitational background, the corresponding Maxwell equations. We get the associated free-wave equation for fields $`\stackrel{}{E}`$ and $`\stackrel{}{B}`$ and discuss the corresponding electrostatic regime. That is, the situation in a gravitational background in which the field $`\stackrel{}{E}`$ does not induce the field $`\stackrel{}{B}`$ and vice-versa.

And we show that, in a particular situation, electric field $`\stackrel{}{E}`$ may depend on time but does not induce magnetic field $`\stackrel{}{B}`$ with or without sources. Finally we show particular electrostatic and magnetostatic solutions for Maxwell equations without sources in a gravitational background. These solutions are such that the electric field, $`\stackrel{}{E}_0`$ and the magnetic field, $`\stackrel{}{B}_0`$, both satisfy the corresponding free-wave equations. Moroever these fields behave as the phenomenology in Plasma Astrophysics has suggested, that is, $`\stackrel{}{E}_0`$ and $`\stackrel{}{B}_0`$ are parallel fields. Thus, the associated electromagnetic wave is stationary, has a null vector Poynting and does not propagate energy.

## 2 Electromagnetic and Gravitational Coupling

####

The action for the gravitaional and electromagnetic coupling is written as

$$𝒮=\sqrt{g}\left(\frac{1}{4}F_{\mu \nu }F^{\mu \nu }\right)d^4s,$$
(2.1)
requiring stationary action $`(\delta S=0)`$ the corresponding Maxwell inhomogeneous equations in a gravitational background are

$$𝒟_\mu F^{\mu \nu }=𝒥^\mu ,$$
(2.2)
where the covariant derivative above is given as:

$$𝒟_\mu F^{\mu \nu }=_\mu F^{\mu \nu }+\mathrm{\Gamma }_{\mu \lambda }^{\mu \nu }+\mathrm{\Gamma }_{\mu \lambda }^\nu \mathrm{\Gamma }^{\mu \lambda }=𝒥^\nu .$$
(2.3)
The corresponding homogeneous Maxwell equations are

$$𝒟_\mu F_{\nu \rho }+𝒟_\nu F_{\rho \mu }+𝒟_\rho F_{\mu \nu }=0.$$
(2.4)
The connection terms of this equation cancel each other in such a way that this equation is the usual homogeneous Maxwell equation

$$_\mu F_{\nu \rho }+_\nu F_{\rho \mu }+_\rho F_{\mu \nu }=0.$$
(2.5)
Finally, to analyse these equations we adopt the F.R.W. cosmological metric

$$dS^2=dt^2\left(a(t)\right)^2\left\{\left(1Ar^2\right)^1dr^2+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\phi ^2\right\},$$
(2.6)
where the term $`a(t)`$ is the scale factor and the constant $`A`$ may assume values $`A=1,0,1`$. Each value represents the associated curvature of F.R.W. spatial metric.

The second term $`\mathrm{\Gamma }_{\mu \lambda }^\mu `$ of the covariant derivative equation (2.3) can be written as

$$\mathrm{\Gamma }_{\mu \lambda }^\mu =\frac{h}{x^\lambda }$$
(2.7)
where

$$h=\mathrm{ln}\sqrt{\stackrel{~}{g}}$$
(2.8)
and $`\stackrel{~}{g}`$ is the metric determinant. In terms of the electric field, $`E^i=F^{0i}`$, and the magnetic field $`B^i=\epsilon ^{ijk}F_{jk}`$, the Maxwell eq. (2.4) and (2.3) are explicitly given by

$$\stackrel{}{}\stackrel{}{E}=\rho _{(\stackrel{}{x})}+\stackrel{}{}h\stackrel{}{E},$$
(2.9)
$$\stackrel{}{}\stackrel{}{B}=0,$$
(2.10)
$$\stackrel{}{}\times \stackrel{}{E}=\frac{\stackrel{}{B}}{t},$$
(2.11)
$$\stackrel{}{}\times \stackrel{}{B}=\stackrel{}{J}(\stackrel{}{x})+\frac{\stackrel{}{E}}{t}\frac{h}{t}\stackrel{}{E}+\stackrel{}{}h\times \stackrel{}{B}.$$
(2.12)
Taking the divergence of equation (2.12) and using eq. (2.9) we get the continuity equation in a gravitational background

$$\frac{\rho }{t}\frac{h}{t}\rho +\stackrel{}{}\stackrel{}{J}\stackrel{}{}h\stackrel{}{J}=0.$$
(2.13)
which is expressed in a covariant way as

$$𝒟_\mu J^\mu =0,$$
(2.14)
where the covariant derivative is $`𝒟_\mu =_\mu _\mu h`$.

Without electromagnetic sources $`(\rho =0;\stackrel{}{J}=0)`$ the free electromagnetic field equations in a gravitational background are:

$$\stackrel{}{}\stackrel{}{E}=\stackrel{}{}h\stackrel{}{E},$$
(2.15)
$$\stackrel{}{}\stackrel{}{B}=0,$$
(2.16)
$$\stackrel{}{}\times \stackrel{}{E}=\frac{\stackrel{}{B}}{t},$$
(2.17)
$$\stackrel{}{}\times \stackrel{}{B}=\frac{\stackrel{}{E}}{t}\frac{h}{t}\stackrel{}{E}+\stackrel{}{}h\times \stackrel{}{B}.$$
(2.18)
Taking the curl of equation (2.17), using equations (2.15) and (2.18) and, since $`\stackrel{}{}h`$ does not depend explicitly on time, we get the free-wave equation in a gravitational background for the electric field:

$$^2\stackrel{}{E}\frac{^2\stackrel{}{E}}{t^2}=\stackrel{}{}(\stackrel{}{}h\stackrel{}{E})\frac{\left({\displaystyle \frac{h}{t}}\stackrel{}{E}\right)}{t}\stackrel{}{}h\times (\stackrel{}{}\times \stackrel{}{E}).$$
(2.19)
Similarly, taking the curl of eq. (2.18), and using eq. (2.16) and, since $`{\displaystyle \frac{h}{t}}`$ does not depend explicitly on the spatial coordinates, we get the free-wave equation in a gravitational background for the magnetic field:

$$^2\stackrel{}{B}\frac{^2\stackrel{}{B}}{t^2}=\stackrel{}{}\times (\stackrel{}{}h\times \stackrel{}{B})\frac{h}{t}\frac{\stackrel{}{B}}{t}.$$
(2.20)
## 3 Analysis of Maxwell Equations in a Gravitational Background

####

Now we analyze the necessary conditions to obtain electrostatic and magnetostatic solutions in the gravitational background, that is, the conditions that the electric field $`\stackrel{}{E}`$ and the magnetic field $`\stackrel{}{B}`$ must satisfy in order that the electric field $`\stackrel{}{E}`$ does not induce the magnetic field $`\stackrel{}{B}`$ and vice-versa.

¿From eq. (2.11) it is clear that the gravitational field does not modify the Faraday’s law and so it is simple to conclude that the magnetic field $`\stackrel{}{B}`$ must be stationary in order not to induce an electric field $`\stackrel{}{E}`$. On the other hand, the Ampère Law (2.12) has been modified by gravitation. Thus, it is simple to see that even if the electric field is stationary $`\left({\displaystyle \frac{\stackrel{}{E}}{t}}=0\right)`$ it may induce a magnetic field $`\stackrel{}{B}`$ from the gravitation term $`{\displaystyle \frac{h}{t}}\stackrel{}{E}`$. For example, if we have a non-null stationary electric field $`\left({\displaystyle \frac{\stackrel{}{E}}{t}}=0;\stackrel{}{E}0\right)`$ and no current density $`(\stackrel{}{j}=0)`$, eq. (2.18) becomes $`\stackrel{}{}\times \stackrel{}{B}={\displaystyle \frac{h}{t}}\stackrel{}{E}+\stackrel{}{}h\times \stackrel{}{B}`$ and so it implies that we necessarily have a non-null magnetic field $`\stackrel{}{B}`$, since $`\stackrel{}{B}=0`$ is not solution for this equation. From eq. (2.12) it is clear that even if the electric field $`\stackrel{}{E}`$ depends on time it shall not induce a magnetic field $`\stackrel{}{B}`$ provided the condition below is satisfied:

$$\frac{\stackrel{}{E}}{t}\frac{h}{t}\stackrel{}{E}=0.$$
(3.1)
¿From condition (3.1) and from eq. (2.9) it is simple to verify that the electric field $`\stackrel{}{E}(\stackrel{}{r},t)`$ and the charge density $`\rho (t)`$ can be written as:

$$\stackrel{}{E}(\stackrel{}{r},t)=(a(t))^3\stackrel{}{E}_0(\stackrel{}{r}),$$
(3.2)
and

$$\rho (\stackrel{}{r},t)=(a(t)^3)\rho _0(\stackrel{}{r});$$
(3.3)
where $`\stackrel{}{E}_0(\stackrel{}{r})`$ is a stationary vector field and $`\rho _0(\stackrel{}{r})`$ is stationary “charge density”. Furthermore the current density $`\stackrel{}{J}`$ must satisfy the equation:

$$\stackrel{}{}\stackrel{}{J}\stackrel{}{}h\stackrel{}{J}=0$$
(3.4)
If the conditions (3.1)-(3.4) are verified and if the magnetic field $`\stackrel{}{B}`$ is static, the electrostatic and magnetostatic field-equations become

$$\stackrel{}{}\stackrel{}{B}=0,$$
(3.5)
$$\stackrel{}{}\times \stackrel{}{B}=\stackrel{}{j}+\stackrel{}{}h\times \stackrel{}{B},$$
(3.6)
$$\stackrel{}{}\stackrel{}{E}_0=\rho _0(\stackrel{}{r})+\stackrel{}{}h\stackrel{}{E}_0,$$
(3.7)
$$\stackrel{}{}\times \stackrel{}{E}_0=0.$$
(3.8)
Any solution of equation (3.4)-(3.6) for the charge current density $`\stackrel{}{J}(\stackrel{}{r},t)`$ and for the magnetic field, $`B(\stackrel{}{r},t)`$, can be combined with any solution of equations (3.2), (3.3), (3.7) and (3.8). So under these conditions the magnetic field, $`B(\stackrel{}{r},t)`$, does not induce the electric field, $`E(\stackrel{}{r},t)`$, and vice-versa. This way we call these equations magnetostatic and electrostatic, but we must remember the electric field is not static but depends on time according to eq. (3.2).

## 4 Solutions for Maxwell Equations without Sources in a Gravitational Scenery

####

We now consider the Maxwell equations (2.9)-(2.12) without sources $`(\rho =0;\stackrel{}{J}=0)`$ with the condition that the electric field $`\stackrel{}{E}(\stackrel{}{r},t)`$ satisfies eq. (3.1) and that the magnetic field $`\stackrel{}{B}`$ does not depend explicitly on time. The Maxwell equations for the electrostatic and magnetostatic field-equations without sources in a gravitational field (3.5)-(3.8) are given by:

$$\stackrel{}{}\stackrel{}{B}=0,$$
(4.1)
$$\stackrel{}{}\times \stackrel{}{B}=\stackrel{}{}h\times \stackrel{}{B},$$
(4.2)
$$\stackrel{}{}\stackrel{}{E}_0=\stackrel{}{}h\stackrel{}{E}_0,$$
(4.3)
$$\stackrel{}{}\times \stackrel{}{E}_0=0.$$
(4.4)
A particular solution in spherical coordinates is

$$E_0(\stackrel{}{r},t)=\frac{\alpha }{r\mathrm{sin}\theta }\stackrel{}{E}(\stackrel{}{r},t)=\alpha \frac{(a(t))^3}{r\mathrm{sin}\theta }\widehat{\phi }$$
(4.5)
and

$$\stackrel{}{B}(\stackrel{}{r},t)=\beta \frac{r}{\sqrt{1Ar^2}}\widehat{\phi }$$
(4.6)
where $`\alpha `$ and $`\beta `$ are constants and $`\widehat{\phi }`$ is the unit azimutal vector: $`\widehat{\phi }=(\widehat{i}\mathrm{sin}\widehat{\phi }+\stackrel{}{j}\mathrm{cos}\phi )`$.

A simple substitution show that these fields satisfy all the Maxwell equations without sources (2.15)-(2.18). It is interesting to point out that the electric and magnetic fields are parallel and satisfy the free-wave eq. (2.19) and eq. (2.20). The electromagnetic wave has parallel electric and magnetic fields, has a null Poynting vector, it is a stationary wave and it does not propagate energy as suggested by the Astrophysical Plasma phenomenology.

The term $`\stackrel{}{}h`$ and the unitary azimutal vector $`\widehat{\phi }`$ are perpendicular vectors, so $`(\widehat{\phi },\stackrel{}{}h,\widehat{\phi }\times \stackrel{}{}h)`$ form a complete set and any vector can be written as:

$$\stackrel{}{V}=V_\phi \widehat{\phi }+V_h\stackrel{}{}h+V_{\phi h}(\widehat{\phi }\times \stackrel{}{}h)$$
(4.7)
The vector components $`V_\phi `$, $`V_h`$ and $`V_{\phi h}`$ may depend on time and spatial coordinates. General solutions $`E(\stackrel{}{r},t)`$ and $`\stackrel{}{B}(\stackrel{}{r},t)`$ for Maxwell equations in a gravitational background using the ansatz (4.7) are being considered.

We conclude that in astrophysical plasma it is important to consider electromagnetic and gravitational coupling and this coupling modifies the free-wave equations for electric and magnetic fields.

### Acknowledgements:

####

The authors acknowledge J.A. Helaÿel-Neto for discussions and comments. Thanks are also due to G.O. Pires Pires for reading the manuscript. Carlos Pinheiro is partially supported by the Conselho Nacional de Pesquisa – CNPq/Brazil.


text:
# 1 Angular distributions of the two proton components from pp collisions included in the simulation, for low secondary proton energies: Parametrization of ref [] data (dashed line), and parametrization from ref [] (solid line). In the upper right plot, the sum of the two components is compared to the data from [].

ISN Grenoble report 00/48

Origin of the high energy proton component below the geomagnetic cutoff in near earth orbit

L. Derome, M. Buénerd, A. Barrau, A. Bouchet, A. Menchaca-Rocha 111Permanent address: Instituto de Fisica, IFUNAM, Mexico, and T. Thuillier

Institut des Sciences Nucléaires, IN2P3, 53 av. des Martyrs, 38026 Grenoble cedex, France

Abstract: The high flux proton component observed by AMS below the geomagnetic cutoff can be well accounted for by assuming these particles to be secondaries originating from the interaction of cosmic ray protons with the atmosphere. Simulation results are reported.

The existence of a high flux proton component below the earth geomagnetic cutoff (GC), was reported recently by the AMS collaboration . The observed spectrum has an intensity about equal to the cosmic primary proton flux above GC. It has a maximum at low momenta and extends from the momentum threshold of the spectrometer at $``$0.4 GeV/c, up to the GC around 10 GeV/c in the equatorial region, where the high momentum tail merges with the primary proton spectrum above GC.

It has been demonstrated long ago that such subGC component cannot be part of the primary proton flux . It has to be a secondary product of the primary cosmic ray (CR) flux on earth. This note reports on an attempt both to identify which part of the incident CR flux on earth generates this spectrum and by which mechanism, and to investigate its dynamical status with respect to the trapping process observed for other circum-terrestrial populations of particles.

A similar phenomenon was observed previously with the very high intensity low energy proton flux measured in the Van Allen belts, at kinetic energies typically below 100 MeV and for distances typically beyond 1.2 earth radius, which was interpreted as originating from secondary neutrons produced on atmospheric nuclei by CR protons . Many low energy neutrons are produced by nuclear evaporation, or spallation, or preequilibrium nuclear emission, which can constitute the source term driving the equilibrium of this population.

The situation is quite different for the population of higher momentum protons, beyond 0.5 GeV/c, observed here. It is too high in energy to be accounted for by the above quoted mechanism, neither could it be by neutrons from the direct charge exchange reaction $`pnnp`$, induced by high energy CR protons of nuclei, because of highly unfavored kinematics and small cross sections of the diffractive production mechanism, and very small probability of neutron decay below the altitude of AMS.

The two most significant sources of secondary protons over the observed energy range are the $`p+𝒜`$ and 4He$`+𝒜`$ collisions, $`𝒜`$ standing for atmospheric nuclei, the CR flux of higher nuclear masses being too small to contribute significantly.

Considering primary protons as a likely dominant source of the observed subGC protons, two extreme situations are met over the explored range of latitudes: In the equatorial region, primary and secondary spectra do not overlap, incident protons are above GC, then with E$``$10 GeV, and produce a spectrum of secondaries (below GC) with E$``$10 GeV, whereas in the polar region (no GSC) primary and secondary protons are totally mixed in a spectrum extending down to the 0.1 GeV kinetic energy range. This latter feature implies that a variety of reaction mechanisms, from nuclear in the low energy range to subnucleonic hadronic in the upper range, could contribute to generate the observed subGC population.

Experimentally, the $`p𝒜pX`$ reaction on light nuclei has been studied in several experiments over a broad range of incident energies , (see  for inclusive $`pppX`$ studies) providing a body of data rich enough for systematic trends to be identified and for a consistent physical picture of the dominant mechanisms to emerge. Two components were observed in the measured proton spectra, which could be qualitatively identified to two different regimes of collision. One is forward peaked, dominated with leading particle effects, and corresponds to peripheral collisions and low orders of multiple intranuclear proton-nucleon collisions. The closeness of the experimental forward proton multiplicites measured in $`p𝒜`$  and $`pp`$  collisions support this picture and shows in addition that protons measured in the first reaction are dominantly hadronic products of $`pp`$ collisions. This is also consistent with the known proton mean free path in nuclei . The other component corresponds to more central collisions and higher order type of multiple intranuclear collisions. In the following, these two components will be referred to as quasi elastic (QE) and deep inelastic (DI) respectively, following ref . In ref , the first component was observed for low energy protons only in the angular range allowed by quasi free scattering kinematics, whereas the other component was observed also at backward angles in the region of zero rapidity, kinematically not accessible to free proton-nucleon processes. The two components have been found to exhibit significantly different dependence on the nuclear target mass, supporting the above description .

In this context, the following features of the reactions relevant to the present study, have been observed: a) The invariant cross section for proton production versus transverse kinetic energy, can be described with a simple exponential over the transverse mass, with a slope weakly dependent on the rapidity . b) At backward angles, a limiting production mechanism has been observed for incident energies above 8.5 GeV/c , the cross-section showing little incident energy dependence, with the onset of the limiting behaviour occuring as low as 1-2 GeV kinetic energy per nucleon for protons produced within the range 0.4-1 GeV/c . It was shown in ref  that the target-like proton production can be described with a simple parametrization over the full angular range (see ref  and the quoted refs in this paper for the theoretical aspects).

The inclusive spectrum of protons at the altitude of AMS (390-400km) has been calculated by means of a computer simulation program built to this purpose. CR particles are generated with their natural abundance and momentum distributions  corrected for the solar modulation effect , at a distance of 5 terrestrial radii, or 2 Störmer length units  in order to ensure the generation point to be always in appropriate initial conditions for the earth magnetic field. They are propagated inside the earth magnetic field , using 4th order adaptative Runge-Kutta integration of the equation of motion. Particles are allowed to interact with atmospheric nuclei (mainly 14N and 16O, see ref  for the model of atmosphere) and produce secondary protons with cross sections and multiplicities as discussed below. Each secondary proton is then propagated and allowed to collide as in the previous step. A reaction cascade can thus develop through the atmosphere. The reaction products are counted when they cross the virtual sphere at the altitude of the AMS spectrometer, upward and downward. Particles undergo energy loss by ionisation before and after the interaction. Multiple scattering effects have not been included at this stage. Each event is propagated until the particle disappears by either colliding with a nucleus, or being stopped in the atmosphere, or escaping to outer space beyond twice the production altitude. A 2000 s time out protection in the program has never been called for the simulated sample. It must be noted that particles are counted each time they cross the sphere of detection altitude. The contributions of trapped particles are thus weighted statistically with their numbers of crossings, which increases in proportion their contribution to the final spectrum. Some empirical cuts over the initial particle kinematics allowed significant computer time saving by rejecting kinematically irrelevant events.

The secondary nucleon spectrum generated has to cover two orders of magnitudes in kinetic energy, between about 100 MeV and 10 GeV. The main component of proton production cross section was obtained by means of analytical relations fitted to the 14.6 GeV $`p+Be`$ data  using the slope parameters given in this reference for the transverse mass distribution, whereas the slowly varying rapidity distributions were fitted with a polynomial. The scaling properties have been checked with the FRITIOF/PYTHIA (Lünd) event generator . Since this generator is not expected to account for the very low energy and backward proton emission (target-like to negative rapidities), this latter component was incorporated using of the parametrization given in ref . The respective contributions to the total multiplicity-weighted proton production cross-section, were 352 mb for the QE component and 88 mb for the DI components. Cross sections on atmospheric nuclei were renormalized from the original data or parametrizations obtained on different nuclei, using ratios of geometrical cross sections. Figure 1 shows the two proton components as a dashed line (QE), and solid line (DI) from $`p+C`$ at 7.5 GeV. For 125 MeV protons, it is seen that the data  are pretty well reproduced by the sum of the two components, providing sound grounds to the present calculations. For each event, the proton multiplicity was generated using a Poisson distribution with the experimental mean value $`\mu `$=1.7 from ref . Neutron spectra and multiplicities were taken identical to proton’s.

A sample of 6 107 particles was generated over the terrestrial sphere (120 hours of CPU computing time), of which 2.7 106 reached the atmosphere, corresponding to a sampling time of the incident flux of the order of 0.5 ps. A cut requiring the accepted events to be inside the angular acceptance of AMS  was applied to the detected protons. Figure 2 shows some basic features of the simulated sample. The four histograms show respectively, the distribution of the altitude of the interaction points (upper left), the thickness of matter crossed by particles (upper right), the lifetime of particles between production and detection (lower left), and the rank of the event in the atmospheric collision cascade (lower right).

The distribution of the proton production altitude has a maximum around 25 km (mean value 39 km), and then decays approximately exponentially, following the atmospheric density profile. The thickness of matter crossed by the particles is distributed over a range extending up to about 600 g/cm2, i.e., about 10 collision lengths, as it could be expected from the known range of the incident particles. The lifetime of the particles extends up to around 600 s. The insert shows the distribution between 0 and 10 s. The lifetime vs momentum correlation obtained (not shown) is in qualitative agreement with the backtracing calculations reported in .

It is seen on the lower right histogram that more than 50% of the detected protons are not produced in the first collision, and originate from up to the 10th collision generation in the atmospheric cascade for a given incident proton. An expected correlation is observed between collision rank and production altitude.

Figure 4 displays an important feature of these results with the scatter plot of the detection latitude versus the number of crossings of the detection altitude. It shows that events with large crossing multiplicities are strongly correlated with low latitudes, thereby showing that a population of trapped particles confined in the equatorial region $`\pm `$0.3 rad is expected to be observed experimentally.

In addition, the following other features are observed:
\- A population of longlived particles, with flight times up to several 102s, is found in the polar region, more than 95% of them having only one or two crossings, i.e., having very long flight path between two crossings. The energy distribution of these particles is the same as that of particles at lower latitudes. Their trajectories then make large loops away from earth.
\- The angular distribution of the detected protons with respect to the incident proton direction, resulting from the folding of the (forward peaked) cross section of the proton production reaction with the (outward peaked) momentum acceptance of the magnetospheric transport system, is found to be approximately evenly distributed over the forward hemisphere and decaying rapidly at backward angles.

Figure 4 shows the experimental kinetic energy distributions of downward (left) and upward (right) protons measured for some bins of latitude, compared to the results of the simulation. No free parameter is used for normalization to the data: The calculated results are entirely determined by the physics input to the calculation. It is seen that the agreement between the data and the simulation result is remarkably good, at all latitudes and for both the inward and outward flux. In particular, the cutoff region is particularly well reproduced, which indicates that the processing of the particle dynamics and kinematics is good. The shaded histograms in the figure correspond to secondary particles in the simulation. The fraction of events originating from the DI component of the proton production cross section described previously, vary from about 10% in the equatorial region up to 25% in the polar region, with a momentum distribution peaking at low kinetic energy and distributed below 500 MeV.

It must also be noted that the broad enhancement observed experimentally in the 1 GeV region in the subGC spectrum for the equatorial latitude compared to the other latitudes, is quite well reproduced. It has been verified that this enhancement is mainly due to the population of trapped particles above 10 crossings as shown on figure 4. The energy distribution of this population is confined below about 3 GeV since particles circling vertically in the earth magnetic field around the mean AMS altitude (380 km) have a Larmor radius of approximately 110 km per GeV/c.

The calculated flux at very low momenta shows a more or less significant deficit compared to the data for most latitudes. Since more than 50% of the detected particles originate from secondary nucleon interaction with the atmosphere, the concerned incident energy range extends down to values where contributions of nuclear processes not taken into account by the event generator, could account for the observed deficit. A defect of the particle trapping processing at low momenta is also possible however.

The 4He cosmic flux should also contribute to the studied population. This contribution should roughly scale with the ratio of incoming flux of 4He and protons, since after the primary interaction, the atmospheric cascade should develop the same way. It is therefore expected to be about one order of magnitude smaller than the cosmic proton contribution , which would not change significantly the results. Detailed investigations are in progress on the above points.

The AMS measurements constitute a new input to the calculations of atmospheric neutrino production (see refs  for example). New calculations are clearly required to evaluate the effect of the observed subGC component of the proton flux on the neutrino flux. It must be observed that the two proton populations, above and below GC, do not have the same interaction probability and angular distribution, and then will not affect equally the atmospheric neutrino flux.

In summary it has been shown that the high intensity flux of protons observed below the geomagnetic cutoff by the AMS experiment can be well reproduced by assuming that this flux originates from the interaction of the primary proton Cosmic Ray flux with the atmosphere. The technical means developed in this work are being applied to the investigation of the other particle populations observed by AMS. They are also potentially useful for similar calculations concerning other experiments such as the atmospheric antiproton component in satellite and balloon experiments, as well as for many other related processes and particle populations. Detailed results will be reported in a more comprehensive paper in preparation.


text:
# Muon capture on nuclei with 𝑁>𝑍, random phase approximation, and in-medium renormalization of the axial-vector coupling constant

## Abstract

We use the random phase approximation to describe the muon capture rate on 44Ca,48Ca, 56Fe, 90Zr, and 208Pb. With 40Ca as a test case, we show that the Continuum Random Phase Approximation (CRPA) and the standard RPA give essentially equivalent descriptions of the muon capture process. Using the standard RPA with the free nucleon weak form factors we reproduce the experimental total capture rates on these nuclei quite well. Confirming our previous CRPA result for the $`N=Z`$ nuclei, we find that the calculated rates would be significantly lower than the data if the in-medium quenching of the axial-vector coupling constant were employed.

The capture of a negative muon from the atomic $`1s`$ orbit,

$$\mu ^{}+(Z,N)\nu _\mu +(Z1,N+1)^{}$$
(1)
is a semileptonic weak process which has been studied for a long time (see, e.g., the reviews by Walecka  or Mukhopadhyay  and the earlier references therein). The total capture rate has been measured for many nuclei ; in some cases the partial capture rates to specific states in the daughter nucleus have been determined as well.

The nuclear response in muon capture is governed by the momentum transfer which is of the order of the muon mass. The energy transferred to the nucleus is restricted from below by the mass difference of the initial and final nuclei, and from above by the muon mass. The phase space and the nuclear response favor lower nuclear excitation energies. There is an intimate relation between the inclusive muon capture rate and the cross section for the antineutrino-induced charged-current reactions; both are governed by the same nuclear matrix elements and proceed from the same initial to the same final states.

Since the experimental data are quite precise, and the theoretical techniques of evaluating the nuclear response in the relevant regime are well developed, it is worthwhile to see to what extent the capture rates are theoretically understood. This is not only interesting per se, but should be viewed as a more general test of our ability to describe semileptonic weak charged-current reactions.

A step in this direction was undertaken by us several years ago in Ref.  where the continuum random phase approximation (CRPA) was used to describe muon capture on the $`N=Z`$ nuclei 12C, 16O and 40Ca. We showed that the method allows us to reproduce the experimental total capture rates on these nuclei to better than 10% using the free nucleon weak form factors and two different residual interactions. In particular, it was not necessary to apply the in-medium quenching of the axial vector coupling constant. This is contrary to various well-known indications that the axial-vector coupling constant $`g_A`$ in nuclear medium is reduced from its free nucleon value of $`g_A=1.26`$ to the value of $`g_A1`$ when one analyses the data on beta decay between low-lying states of the $`(sd)`$ shell nuclei  and $`(pf)`$ shell nuclei . In addition, the “missing Gamow-Teller strength” problem, as revealed in the interpretation of the forward-angle $`(p,n)`$ and $`(n,p)`$ charge-exchange reactions , is also often quoted as evidence for quenching of $`g_A`$. Note that the Gamow-Teller (GT) strength is concentrated in the giant GT resonance at excitation energies not very far from the energies involved in the muon capture, although this latter process is usually dominated by the transitions to the negative parity spin-dipole states. That is so in particular in the double-magic nuclei 16O and 40Ca where the GT strength is strongly suppressed, while in 12C it is essentially exhausted by the transition to the analog of the $`T=1`$ state at 15.11 MeV in 12C, whose contribution to the muon capture rate has been subtracted in . It is thus of interest to inquire whether a similar quenching applies in muon capture over a broader range of nuclei. In this way, we can test whether the quenching is a general phenomenon, applicable to more than just GT transitions.

In this paper we extend the previous calculation  to heavier nuclei, in particular nuclei with the neutron excess, i.e., with a nonvanishing value of the initial isospin. In  the Continuum Random Phase Approximation (CRPA) was used, a method shown to be successful for the description of the nuclear response to weak and electromagnetic probes . The method combines the usual RPA treatment with the correct description of the continuum nucleon decay channel. For heavier $`N>Z`$ nuclei that method is computationally quite demanding. Moreover, it can describe only transitions to states in the daughter nucleus above the particle emission threshold. In order to evaluate the total capture rate one has to include also the transitions to bound states, which contribute relatively more, and are typically not experimentally separated in the heavier nuclei.

If one is interested in the total capture rate, the numerically simpler standard RPA is just as good. As an additional bonus, it avoids the distinction between the bound and unbound states. That the two methods, CRPA and standard RPA, are equivalent for the present purpose is demonstrated in Fig. 1 for the case of 40Ca. We show in the upper part the differential capture rate as a function of the excitation energy in the final nucleus. While the CRPA (for details see Refs. ) is characterized by a continuous curve, nonvanishing everywhere above the threshold, the standard RPA is characterized by the ‘picket fence’, since there is a finite number of discrete final states. The similarity of both methods is even better seen in the lower part of Fig. 1 which shows the integrated rate up to a given excitation energy $`E^{}`$. There appears to be a slight systematic shift of a few MeV (caused by the bound state contribution, presumably), but the final capture rates, and the typical excitation energies, are remarkably similar. Thus, we use the standard RPA for the evaluation of the muon capture in the selected $`N>Z`$ nuclei.

For the calculation in the present work we used the phenomenological Landau-Migdal force with parameters that has been shown to be applicable for a wide range of nuclei . All single-particle states below the Fermi level were included and two oscillator shells above it were taken into account. Using again 40Ca as a test case, we checked that adding or subtracting in the calculation few subshells above the Fermi level does not visibly change the muon capture rate. However, in 208Pb enlarging the single particle space leads to an increase of the capture rate by about 5%. The free nucleon form factors were used to describe the weak nuclear current. In particular, the unquenched axial vector coupling constant $`g_A(0)=1.26`$ was used.

Muon capture also depends on the induced pseudoscalar hadronic weak current. At the free nucleon level the corresponding coupling constant is determined by the Goldberger-Treiman relation

$$F_P(q^2)=\frac{2M_pg_A(0)}{m_\pi ^2q^2},$$
(2)
where $`m_\pi `$ is the pion mass and $`g_A(0)=1.26`$. (In muon capture one often uses a dimensionless quantity $`g_P=m_\mu F_P(q^2)`$ at the relevant momentum transfer $`q^20.9m_\mu ^2`$, such that $`g_P8.4`$ for free protons.) In nuclear medium $`F_P`$ can be again renormalized, and this renormalization does not necessarily obey the Goldberger-Treiman relation . We have shown in our previous work that the total muon capture rates are not sensitive enough to the various choices of $`F_P`$ renormalization. Consequently, throughout this work we use the Goldberger-Treiman relation.

The calculated total capture rates are collected in Table I and compared with the data . For comparison we also show in Tab. I the earlier results  for the $`N=Z`$ nuclei, evaluated with the same residual Landau-Migdal force. Among the nuclei in Table I, 56Fe is the only one with a substantial contribution of the Gamow-Teller transitions, which are known to be quenched. Hence the $`1^+`$ part of the capture rate was quenched by the empirical factor of 2.54 obtained by comparing the total experimental and calculated GT strength. This reduced $`1^+`$ rate was then combined with the unquenched rate from the other multipoles in column 3 of Table I. Clearly, a better agreement with experiment for all considered nuclei is achieved when the full value of $`g_A`$ is used. The quenched value leads to an obvious underestimate of the muon capture rate.

In Fig. 2 the fractional contributions of different multipoles are shown for 16O, 48Ca, 90Zr, and 208Pb. For the closed shell nucleus 16O the negative parity $`1^{}`$ and $`2^{}`$, and to lesser extent $`0^{}`$, multipoles dominate, as expected. In the intermediate mass nuclei, 48Ca and 90Zr, the dominance of these multipolarities is less pronounced. Finally, in 208Pb the positive parity $`1^+`$ and $`2^+`$ multipoles give the largest contribution, since the negative parity proton hole-neutron particle states are blocked. Note, however, that these $`1^+`$ and $`2^+`$ states correspond to the 2$`\mathrm{}\omega `$ excitations involving different major shells, not the usual 0$`\mathrm{}\omega `$ Gamow-Teller transitions.

The dominance of dipole transitions for the $`\mu ^{}`$ capture process in the $`p`$ and $`s,d`$ shell nuclei is a well known phenomenon. It can be exploited in the generalized Goldhaber-Teller model , where all dipole and spin-dipole strength is concentrated in a single collective state. The Goldhaber-Teller model has been applied to the muon capture already a long time ago . Here we have repeated the calculation with one modification. Instead of using as the energy of the single collective state the energy $`E_D`$ of the giant dipole resonance in the initial nucleus, we placed the strength into a state at the excitation energy $`E^{}=E_DE_{IAS}`$ in the final nucleus, where $`E_{IAS}`$ is the energy of the isobar analog state with $`T=1`$ in the initial $`N=Z`$ nucleus. Such an assignment puts the strength close to the centroid of the excitation spectrum obtained in the RPA. The method is rather crude, since more detailed calculations clearly show that the spin-dipole strength is spread over a sizable energy interval. Nevertheless, the results displayed in Table II, obtained again with the full value of $`g_A(0)=1.26`$, support our conclusion that there is no quenching in these nuclei of the operators of the weak current that change parity.

First forbidden beta decays, in particular those with $`\mathrm{\Delta }I^\pi =0^{}`$, have been often analyzed as a source of information on the enhancement caused by the meson exchange currents. In the context of the work reported here it is worthwhile to quote the work of Warburton and Towner  who analyzed 18 first forbidden beta decays in the lead region. They used a truncated shell model and found that the $`\mathrm{\Delta }I^\pi =1^{}`$ transitions, which are unaffected by the meson exchange currents, do not require any quenching. In fact, their fit to an overall quenching factor for the first forbidden matrix elements results in $`sq_1=0.98\pm 0.05`$, compatible with unity, i.e., with no quenching. That analysis, totally independent of our evaluation of the muon capture, again supports the conclusion about the apparent absence of an appreciable quenching of the parity-changing (i.e., first forbidden) weak current operators.

The unique second forbidden beta decays are governed by a single operator $`r^2[Y_2\sigma ]^{\lambda =3}`$. For only a handful of them the partial decay rates are known. When analyzed , these transitions do not allow one to draw any definitive conclusion about the possible quenching of the corresponding strength. However, as pointed out above, $`\mu ^{}`$ capture rate in 90Zr, and 208Pb, is strongly affected by the 2$`\mathrm{}\omega `$ transitions, whose operators are related to the second forbidden beta decays. Since the agreement between the experimental and calculated rates in these two nuclei, and in particular in 208Pb, is not as good as in the other cases, we cannot make a strong statement regarding the quenching of the positive parity “second forbidden” multipoles based on the muon capture calculations.

In conclusion, the present analysis shows that the CRPA and SRPA methods are capable of describing the total $`\mu ^{}`$ capture rates quite well in a large range of nuclei. The dependence of the muon capture rate on the isospin, the so-called Primakoff rule , is also reasonably well reproduced. There is no indication of the necessity to apply any quenching to the operators responsible for the $`\mu ^{}`$ process. Thus our findings indicate that any in-medium quenching of the axial current matrix elements appears to be restricted to the 0$`\mathrm{}\omega `$ spin changing operators, i.e. to the Gamow-Teller operator.

###### Acknowledgements.

This work was supported in part by the Swiss National Science Foundation, the Danish Research Council, and by the U.S. Department of Energy, Contract #DE-F603-88ER-40397.


text:
# Equipartition and Virial theorems in a nonextensive optimal Lagrange multipliers scenario

## I Introduction

Tsallis’ thermostatistics (TT)  is by now recognized as a new paradigm for statistical mechanical considerations. It revolves around the concept of Tsallis’ information measure $`S_q`$, a generalization of Shannon’s one, that depends upon a real index $`q`$ and becomes Shannon’s measure  for the particular value $`q=1`$.

The Equipartition theorem, first formulated by Boltzmann in 1871, and the Virial one, due to Clausius (1870)  are two pillars of classical physics. They were discussed within the TT framework in , in terms of the Curado-Tsallis nonextensive normalizing treatment , and also in , by recourse to escort probabilities . In the present work we are going to revisit these classical subjects in connection with a recently advanced scheme for dealing with Tsallis’ thermostatistics , that seems to yield illuminating insights into classical themes.

Tsallis’ thermostatistics involves extremization of Tsallis’ entropy by recourse to the celebrated technique of Lagrange. A key TT ingredient is the particular way in which expectation values are computed  and, in such a respect, several proposals have been considered during the last ten years . No matter which recipe one chooses (for $`q1`$), classical phenomena are always reproduced by TT in the limit $`q1`$ . Recently, a new algorithm (to be henceforth referred to as the OLM-one) has been advanced to such an extremizing end that diagonalizes the Hessian associated to the Lagrange procedure  and yields the optimal Lagrange multipliers associated to the input expectation values. A diagonal Hessian enormously facilitates ascertaining just what kind of extreme the Lagrange method is leading to .

A key point in TT considerations is the following one: the entropy constant $`k`$ is usually identified with Boltzmann’s $`k_B`$. However, the only certified fact one can be sure of is “$`kk_B`$ for $`q1`$, which entails that there is room to choose $`k=k(q)`$ in any suitable way. It is seen  that if one chooses

$$k=k_B\overline{Z}_q^{q1},$$
(1)
where $`\overline{Z}_q`$ stands for the partition function, in conjunction with the OLM formalism , the classical harmonic oscillator determines a specific heat $`C_q=k_B`$, so that the classical Gibbs’ result arises without the need to invoke the limit $`q1`$. Additionally, this treatment is able i) to reproduce thermodynamics’ zero-th law , a feat that had eluded previous TT practitioners, and ii) the classical $`q`$-independent value for the mean energy of the ideal gas . It is then not unreasonable to conjecture that many other classical results might be obtained by TT without invoking the $`q1`$ limit.

In the present effort we re-discuss the equipartition theorem and the virial one from an OLM viewpoint. We also revisit the ideal gas problem as an application. More precisely, we perform an OLM treatment of the subject following the canonical ensemble strictures.

## II Brief review of OLM formalism in a classical scenario

The most general classical treatment requires consideration of the probability density $`p(𝐱)`$ that maximizes Tsallis’ entropy

$$\frac{S_q}{k}=\frac{1𝑑𝐱p^q(𝐱)}{q1},$$
(2)
by recourse to the Lagrange’s technique, subject to the foreknowledge of $`M`$ generalized expectation values

$$O_j_q=\frac{𝑑𝐱p^q(𝐱)O_j(𝐱)}{𝑑𝐱p^q(𝐱)},$$
(3)
where $`O_j(𝐱)`$ ($`j=1,\mathrm{},M`$) denote the $`M`$ relevant dynamical quantities (the observation level ), $`q\mathrm{}`$ is Tsallis’ nonextensivity index, $`𝐱`$ is a phase space element ($`N`$ particles in a $`D`$-dimensional space), and $`k`$ is the entropy constant, akin to the famous Boltzmann one $`k_B`$, employed in the orthodox statistics.

Tsallis’ normalized probability distribution , is obtained by following the well known MaxEnt route . Instead of effecting the variational treatment of , involving $`M+1`$ Lagrange multipliers $`\lambda _j`$ (associated to constraints given by the normalization condition together with the $`M`$ equations (3)), the OLM technique follows an alternative path  with Lagrange multipliers $`\lambda _j^{}`$: one maximizes Tsallis’ generalized entropy $`S_q`$ (2)  subject to the modified constraints (“centered” generalized expectation values)

$`{\displaystyle 𝑑𝐱p(𝐱)}1`$ $`=`$ $`0`$ (4)
$`{\displaystyle 𝑑𝐱p(𝐱)^q\left(O_j(𝐱)O_j_q\right)}`$ $`=`$ $`0.`$ (5)
The resulting probability distribution reads

$$p(𝐱)=\frac{f(𝐱)^{\frac{1}{1q}}}{\overline{Z}_q},$$
(6)
where

$$\overline{Z}_q=𝑑𝐱f(𝐱)^{\frac{1}{1q}}.$$
(7)
and

$$f(𝐱)=1(1q)\underset{j}{\overset{M}{}}\lambda _j^{}\left(O_j(𝐱)O_j_q\right)$$
(8)
is the so-called configurational characteristic.

Although the Tsallis-Mendes-Plastino (TMP) procedure originally devised in  overcomes most of the problems posed by the old, unnormalized way of evaluating Tsallis’ generalized mean values , it yields probability distributions that are self-referential, which entails some numerical difficulties. The complementary OLM treatment of  surmounts these hardships. Inspection of (6) shows that the self-reference problem has vanished.

One shows in  that the relation

$$𝑑𝐱p^q(𝐱)=\overline{Z}_q^{1q},$$
(9)
valid under TMP, still holds. Eq. (9) allows one to connect the Lagrange multipliers $`\lambda _j`$ of the TMP procedure  with the corresponding OLM $`\lambda _j^{}`$ via

$$\lambda _j^{}=\frac{\lambda _j}{\overline{Z}_q^{1q}},$$
(10)
and to write the entropy as

$$S_q=k\mathrm{ln}_q\overline{Z}_q,$$
(11)
where $`\mathrm{ln}_qx=(1x_q^{1q})/(q1)`$ has been used.

If now, following  we define

$$\mathrm{ln}_qZ_q^{}=\mathrm{ln}_q\overline{Z}_q\underset{j}{}\lambda _j^{}O_j_q,$$
(12)
and additionally use

$$k^{}=k\overline{Z}_q^{1q},$$
(13)
then

$`{\displaystyle \frac{S_q}{O_j_q}}`$ $`=`$ $`k^{}\lambda _j^{}`$ (14)
$`{\displaystyle \frac{}{\lambda _j^{}}}\left(\mathrm{ln}_qZ_q^{}\right)`$ $`=`$ $`O_j_q.`$ (15)
Equations (14) and (15) constitute the basic Information Theory relations to build up Statistical Mechanics à la Jaynes . Notice that, when $`k`$ is given by (1), Equation (13) leads to $`k^{}=k_B`$ .

Remembering that

$$\frac{S_q}{O_j_q}=k\lambda _j,$$
(16)
one is straightforwardly led to

$$k^{}\lambda _j^{}=k\lambda _j,$$
(17)
which entails that the intensive variables are identical in both alternative pictures, TMP and OLM .

As a special instance of Eq. (17), for the Canonical Ensemble it adopts the appearance

$$k^{}\beta ^{}=k\beta =\frac{1}{T}.$$
(18)
Looking at Eq. (18) one gathers that the temperature $`T`$ is the same for both sets of Lagrange multipliers.

Before tackling the equipartition theorem, some preliminary results are needed, that we discuss next.

## III Normalization considerations

Replacing (6) into (9) we obtain for the partition function the expression

$$\overline{Z}_q=𝑑𝐱f(𝐱)^{\frac{q}{1q}},$$
(19)
and, comparing it to (7), we deduce the (at first sight surprising) relationship

$$𝒟=𝑑𝐱f(𝐱)^{\frac{q}{1q}}𝑑𝐱f(𝐱)^{\frac{1}{1q}}=0,$$
(20)
valid for all $`q`$. That this is indeed so can also be gathered by first recasting (20) in the form

$$𝒟=𝑑𝐱\left[f(𝐱)^{\frac{q}{1q}}\left(1f(𝐱)\right)\right],$$
(21)
and then using Eq. (8) to obtain

$$𝒟=(1q)\underset{j}{}\lambda _j^{}𝑑𝐱f^{\frac{q}{1q}}\left(\widehat{O}_j\widehat{O}_j_q\right),$$
(22)
that, remembering (6), has to vanish identically on account of (5).

We introduce now the factor

$$F(q)\frac{𝑑𝐱f(𝐱)^{\frac{1}{1q}}}{𝑑𝐱f(𝐱)^{\frac{q}{1q}}}=1,$$
(23)
in order to be in a position to use it later on.

It is to be stressed that $`F(q)=1`$ only within the so-called “normalized” framework of . If one uses, instead, the Curado-Tsallis normalization , this is not so ($`F(q)1`$). The OLM formalism does employ the normalized treatment. It must be pointed out that the $`F(q)=1`$-result has also been obtained within a purely quantal, Green function scheme by Lenzi, Mendes, and Rajagopal .

## IV Generalized equipartition and virial theorems

In classical statistical physics, a Hamiltonian dynamical system is described by an appropriate phase space probability distribution $`p(r_i,p_i).`$ Just one assumption will be made on the probability density (PD): that it depend on the phase space variables $`(r_i,p_i)`$ only through the Hamiltonian $`H(r_i,p_i)`$.

Tsallis’ normalized probability distribution  is obtained by maximizing Tsallis’ generalized entropy $`S_q`$ given by Eq. (2), subject to the constraints

$`{\displaystyle 𝑑\mathrm{\Omega }p(r_i,p_i)}`$ $`=`$ $`1`$ (24)
$`{\displaystyle 𝑑\mathrm{\Omega }p^q(r_i,p_i)\left(H(r_i,p_i)H_q\right)}`$ $`=`$ $`0,`$ (25)
where the generalized expectation values

$$H_q=\frac{𝑑\mathrm{\Omega }p^q(r_i,p_i)H(r_i,p_i)}{𝑑\mathrm{\Omega }p^q(r_i,p_i)}$$
(26)
are assumed to be a priori known. They constitute the macroscopic information at hand concerning the system. $`d\mathrm{\Omega }`$ stands for the corresponding phase space volume element

$$d\mathrm{\Omega }=(1/(N!h^{DN}))\underset{i=1}{\overset{DN}{}}dr_idp_i,(i=1\mathrm{}DN),$$
(27)
where $`h`$ is the linear dimension (i.e. the size) of the elementary cell in phase space, and we assume $`_{i=1}^{DN}dr_i=V^N`$ with $`V`$ the system’s volume.

The resulting probability distribution reads

$$p(r_i,p_i)=\overline{Z}_q^1\left[1(1q)\beta ^{}\left(H(r_i,p_i)H_q\right)\right]^{\frac{1}{1q}},$$
(28)
where

$$\overline{Z}_q=𝑑\mathrm{\Omega }\left[1(1q)\beta ^{}\left(H(r_i,p_i)H_q\right)\right]^{\frac{1}{1q}}.$$
(29)
Let $`A(q_{i,}p_i)`$ denote a generic dynamical quantity. Its generalized mean value is given by

$$A_q=\frac{𝑑\mathrm{\Omega }A(r_{i,}p_i)\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{q}{1q}}}{𝑑\mathrm{\Omega }\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{q}{1q}}}.$$
(30)
Let us now assume that our Hamiltonian $`H`$ can be split into two pieces according to

$$H=g_0+g$$
(31)
where $`g`$ is a homogeneous function of degree $`\gamma `$ of $`L`$ canonical variables of which, say $`\nu `$, are generalized coordinates while the remaining ones are $`\mu =L\nu `$ generalized momenta

$$g=g(r__1,\mathrm{},r_\nu ,p_1,\mathrm{},p_\mu ),$$
(32)
and $`g_0`$ does not depend upon these variables. According to Euler’s theorem  we have

$$\gamma g=\underset{i=1}{\overset{\nu }{}}r_i(g/r_i)+\underset{i=1}{\overset{\mu }{}}p_i(g/p_i)$$
(33)
so that the generalized mean value of $`g`$ reads

$$g_q=\frac{1}{\gamma }\left[\underset{i=1}{\overset{\nu }{}}r_i(g/r_i)_q+\underset{i=1}{\overset{\mu }{}}p_i(g/p_i)_q\right].$$
(34)
We shall now discuss in some detail one generic term of this equation. Consider

$$r_k(g/r_k)_q=\frac{𝑑\mathrm{\Omega }r_k(g/r_k)\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{q}{1q}}}{𝑑\mathrm{\Omega }\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{q}{1q}}},$$
(35)
a multi-dimensional $`(2DN)`$ integral. Let us evaluate the integral over $`r_k`$ ranging between $`r_a`$ and $`r_b`$. These values are given by the well-known Tsallis’ cut-off condition : the probability distribution vanishes in those regions of phase space that would make (28) a negative quantity. For the numerator of (35) we have

$$J_{r_a}^{r_b}𝑑r_kr_k(g/r_k)\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{q}{1q}},$$
(36)
so that

$$r_k(g/r_k)_q=\frac{\mathrm{}J𝑑r_1\mathrm{}𝑑r_{k1}𝑑r_{k+1}\mathrm{}𝑑p_{DN}}{𝑑\mathrm{\Omega }\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{q}{1q}}}.$$
(37)
In order to obtain $`J`$ we proceed to an integration by parts. To do so we first notice that, on account of Equations (31) and (32)

$$\frac{}{r_k}\left\{\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{1}{1q}}\right\}=\beta ^{^{}}\frac{g}{r_k}\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{q}{1q}},$$
(38)
which, since the integrated part will vanish because of the above mentioned cut-off condition, leads to

$$J=\frac{1}{\beta ^{^{}}}𝑑r_k\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{1}{1q}}.$$
(39)
Insertion of (39) into (37) yields

$$r_k(g/r_k)_q=\frac{1}{\beta ^{^{}}}\frac{𝑑\mathrm{\Omega }\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{1}{1q}}}{𝑑\mathrm{\Omega }\left[1\beta ^{^{}}(1q)\left(HH_q\right)\right]^{\frac{q}{1q}}},$$
(40)
but, on account of (23) one has

$$r_k(g/r_k)_q=\frac{1}{\beta ^{^{}}}.$$
(41)
It is apparent that each term in the sums appearing in (34) will yield a contribution of the type (41), so that

$$g_q=\frac{L}{\gamma \beta ^{^{}}}=\frac{L}{\gamma }k_BT$$
(42)
which is a generalized version of the equipartition theorem . We have thus arrived to a new an interesting result. Contrary to current belief, we see that the classical result is attained, independently of the $`q`$-value.

According to the canonical equations of motion, $`H/r_i=\dot{p}_i.`$ Hence (41) leads to the statement

$$\underset{i=1}{\overset{DN}{}}r_i\dot{p}_i_q=DNk_BT$$
(43)
which is Clausius’ virial theorem . Again, no dependence upon $`q`$ is to be detected.

As a simple application, consider now the classical ideal gas. We necessarily reproduce the classical results, by virtue of the above considerations: we deal with $`N`$ particles confined within a $`D`$dimensional box of volume $`V`$. Thermodynamical equilibrium at temperature $`T`$ is assumed. The Hamiltonian is

$`H={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{𝐩_i^2}{2m}}+{\displaystyle \underset{i=1}{\overset{N}{}}}U_{wall}(𝐫_i),`$

where $`U_{wall}`$ is a contribution due to the constraints (walls of the container). The constraint forces come into existence when the gas particles collide with the walls. It thus follows that

$$\underset{i=1}{\overset{N}{}}𝐫_i\frac{H}{𝐫_i}=\underset{i=1}{\overset{N}{}}𝐫_i𝐅_{wall}^{(i)}$$
(44)
where $`𝐅_{wall}^{(i)}`$ stands for the force on the i-particle due to the walls of the box.

Following now a well-known argument , it is possible then to express the last member of the above equation in terms of the volume $`V`$ and pressure $`P`$

$`{\displaystyle \underset{i=1}{\overset{N}{}}}𝐫_i𝐅_{wall}^{(i)}`$ $`=`$ $`P{\displaystyle _{wall}}𝐫(d𝐬)=P{\displaystyle _{wall}}𝐫𝑑𝐬`$ (45)
$`=`$ $`P{\displaystyle (𝐫)d^3r}=DPV,`$ (46)
where $`d𝐬`$ denotes a (vector) surface element of the box in Eq. (44). We get

$$\underset{i=1}{\overset{N}{}}𝐫_i\frac{H}{𝐫_i}=DPV$$
(47)
Now, by taking the canonical average and using the virial theorem (unmodified by the nonextensive scenario) (43), we obtain

$$PV=Nk_BT,$$
(48)
i.e., the equation of state for the perfect gas. No $`q`$-dependence is detected.

For the ideal gas , the total energy $`E`$ is a homogeneous quadratic function $`(\gamma =2)`$ of $`DN`$ momenta, which allows one to write, according to (42)

$$H_q=\frac{1}{2}DNk_BT.$$
(49)
## V Conclusions

Classical thermostatistics has been the subject of the present effort.

We have tackled some key issues, namely,

* Virial theorem,
* Equipartition theorem,
* Equation of state of the ideal gas,

and shown that they are, contrary to present belief, reproduced by Tsallis’ thermostatistics independently of the value adopted by the index $`q`$. The present work lends further credence to the hypothesis that most important classical statistical results might be reproduced by Tsallis’ thermostatistics for all $`q`$-values.

###### Acknowledgements.

The financial support of the National Research Council (CONICET) of Argentina is gratefully acknowledged. F. Pennini acknowledges financial support from UNLP, Argentina.