text:
# Tracer dispersion in two-dimensional rough fractures

## I Introduction

Tracer spreading in flows between parallel walls has received considerable attention since the celebrated work of Taylor. This problem is of fundamental importance in a variety of fields, and in particular transport processes through artificial or natural porous media. In general, the dispersion of a tracer carried along by a flowing fluid in a medium of disordered structure, such as hydrocarbon or water reservoirs, involves a combination of convection-dispersion through the microscopic pore space of the rock itself and through macroscopic channels such as fractures. In the simplest case of Poiseuille flow in a Hele-Shaw cell, which in several applications is used to model fractures , the vanishing velocity near the solid gives rise to a large dispersion, quadratic in the Péclet number . However, one frequently encounters systems in which the channel aperture varies along the flow direction, or in which the channel walls have some rugosity. This is the case for fractured rocks, in which fractured boundaries can be usually described as correlated, self-affine fractals. The roughness exponent, is usually found to be close to $`0.8`$, insensitive to the material and fracturing process . Most of the studies dealing with varying channel aperture model the fluctuations as slowly varying periodic functions . Only a few works deal with more realistic models of fractures, such as random rugosities perpendicular to the flow  or self-affine fractures .

In the present paper we present numerical simulations of tracer diffusion and dispersion in self-affine fractures. The simulations are two-dimensional, (two-dimensional fracture surfaces), and the fluid flow and tracer spreading are computed using the lattice-Boltzmann (LB) method. We will first discuss, in section II, the implementation of the LB method to simulate diffusive transport. We shall propose a new boundary condition which improves the accuracy and validity of the method, particularly when describing narrow fractures. In section III we will compare the proposed boundary condition with previous ones, in two different cases. Diffusion (section III A) and Taylor dispersion (section III B) in a two-dimensional straight channel.

After reviewing some basic facts about self-affine surfaces (section IV), we investigate, in section V, how the complex geometry of fracture surfaces affects the diffusive transport (zero mean flow velocity). We also study the dependence of diffusivity on the size of the gap between complementary fracture surfaces. We will show that it is possible to derive an analytic expression, which accounts for the geometric effects on the diffusivity, in the limiting case where the roughness associated with the fracture surface is small compared to the mean aperture.

In section VI we study dispersive transport in two different situations. First, when the two surfaces are simply displaced normally to the mean fracture plane. In this case, we will show that a description in terms of a parameter $`\mathrm{\Lambda }`$, related to dynamically connected pore space, accounts for a large part of the enhancement due to the presence of low-velocity zones in the fracture channel. Secondly, we show that when there is a lateral shift as well as the normal displacement, an increase in dispersion is obtained. Finally, we present a novel result, showing that dispersion is diminished in the presence of convection at low Péclet numbers. This effect is found when the surfaces are shifted and the convective transport is weak enough.

## II Lattice-Boltzmann model with the BGK collision operator

Our goal is to study various aspects of transport in fractures which are sensitive to the fracture roughness. Since we consider convection and dispersion in a highly irregular geometry, the lattice-Boltzmann 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 boundary condition rules should be imposed. We use a version of the LB model first proposed by Qian et. al , with a cubic lattice in 3 dimensions and 19 velocities (D3Q19 in the terminology of ). The collision operator is approximated by a single relaxation parameter $`\lambda `$, the Bhatnagar-Gross-Krook model, and the local equilibrium distribution given in  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. As usual, we define the lattice spacing as the unit of length, and the time step in the simulations as the unit time. In what follows, quantities will be measured in terms of the lattice-spacing and simulation time-step. Note that, since we are concerned with incompressible flows, we do not need to introduce a dimension of mass.

The model used for miscible fluids is a straightforward extension of the model restricted to describe simple fluids. Following the convention of Rothman and Zaleski , we will distinguish between particle types by assuming that the particles are colored. The pseudo-equilibrium color distribution given in  is used, and a single but new relaxation parameter $`\lambda _D`$ (BGK approximation) is used in the equation that describes the advection and diffusion of color. In addition to mass and momentum, color is also a conserved quantity.

The basic variables of the model are the distribution functions $`N_i`$, corresponding to the mean occupation number of particles in the direction $`i`$ at a given node, and $`\mathrm{\Delta }_i`$, describing the relative amount of color . Since the evolution of the total population is independent of the color of the particles, the hydrodynamic equations for mass and momentum may be obtained as in the simple fluids case. In other words, the time evolution of the $`N_i`$’s is independent of the $`\mathrm{\Delta }_i`$’s. These evolution equations are then coupled to the evolution of color, through the local fluid velocity. Thus, the diffusive behavior of the fluids is superimposed on the underlying Navier-Stokes dynamics.

The fact that that the information on mass density and flow is decoupled from the information on color, opens the possibility of imposing different boundary conditions (BC) on the populations describing fluid flow ($`N_i`$) and those describing advective-diffusive transport ($`\mathrm{\Delta }_i`$). Specifically, for the non-slip BC on the solid surface, we shall use the simplest implementation of particle exclusion - the bounce-back rule (BB), where the particles incident on the boundary are propagated back into the directions from which they came. On the other hand, for color concentration a different macroscopic BC is desired, i. e. zero color-gradient normal to the solid surface.

### A Diffusion boundary conditions

As mentioned, we will study tracer dispersion in narrow gaps between self-affine solid surfaces. Thus, in this situation, the BC imposed at the solid surfaces becomes a crucial aspect of the simulation method.

In previous works concerning lattice BGK models for miscible fluids, BC imposed on color concentration at solid boundary were not distinguished from those used to simulate fluid flow. In these situations, the relevance of the BC at solid surfaces varies depending on the particular system under consideration. In some cases, when simulating bulk processes, as in , there is no need to treat the BC separately, and periodic boundary conditions on all distribution functions may be used. In other cases, fully three dimensional problems in Helle-Shaw cells are described in terms of two dimensional LB model by the modification of the forcing to account for the viscous drag of the top and bottom plates of the cell . Therefore, molecular diffusion in the transverse direction is not described and similar BC can be applied to all fields. In the others, where fluid flows in complex geometries  or in narrow channels bounded by solid surfaces , the above-mentioned distinction in BC would be desirable.

Here we will use a different set of BC for color concentration, to ensure zero color-gradient normal to the solid surface. To this end, we implement a mirror-reflection condition or bounce-forward rule (BF), where upon collision only the normal component of the particle velocity is changed. In Figure 1 we show a schematic view of the BF collision rule. The last case was chosen arbitrarily, due to its very simple implementation, among the different possibilities that satisfy the mirror-reflection-type condition. For instance, it would be also consistent to split the incoming concentration into the two neighboring sites.

Physically, BB has the (approximate) effect of making both components of the fluid velocity vanish at the solid, consistent with the non-slip condition. A passive tracer may, however, diffuse along a solid boundary, so the relevant condition is that only the normal component of the flux vanishes at a solid. the bounce-forward rule is the obvious LB realization of the physical boundary condition.

To answer the question whether this BC improves the LB method for miscible fluids we shall present the results of numerical simulations using both BB and BF in two different situations: diffusion and Taylor hydrodynamic dispersion in a two-dimensional straight channel.

## III Comparison between BB and BF rules: numerical simulations

We will exhibit a comparison between BB and BF, in two applications, diffusion and Taylor dispersion in a two-dimensional channel. These situations are particularly important in terms of addressing the validity and accuracy of the above-mentioned boundary conditions, particularly if one is interested in further use of them in complex geometries.

### A Diffusion in a two-dimensional straight channel

Diffusion was studied varying the channel width $`H`$ from $`4`$ to $`40`$ in grid sites, in a channel of length $`L=512`$, with periodic BC used at both ends. In the LB model the diffusion coefficient in bulk is given by,

$$D_m=\frac{1}{3}\left(\frac{1}{\lambda _D}\frac{1}{2}\right)$$
(1)
varying the relaxation parameter $`\lambda _D`$ we have also studied diffusive transport for a set of bulk diffusivities ranging from $`D_m=8.8\times 10^4`$ to $`D_m=1.5`$.

In all the simulations for both boundary rules, the obtained depth-averaged color-concentration displays a Gaussian distribution, and the mean square displacement grows linearly in time. However, the vertical dependence of the concentration, as well as the diffusivity values, strongly depend on the BC used. While concentration remains vertically-homogeneous using BF, in Fig. 2 we show that, using BB yields the undesired effect of a transient vertical variation of color. Two different situations can be distinguished in Figure 2, namely, close to the mean position of tracer particles $`x`$, or far from it. (Note that $`x`$ is constant in time). At $`x`$ (Fig. 2a), the initial concentration gradient in drives tracer particles. Due to the bounce back rule, particles close to the surface remain there longer, and the relative tracer concentration builds up. Shortly after, particles close to the surface diffuse toward the center of the channel, and eventually the concentration gradient in vertical direction is smoothed out. On the other hand, far from the mean position of tracer particles the situation is inverted. Bounce back at solid surfaces, makes the tracer particles mostly arrive from the center of the channel, and a vertical concentration gradient develops. Thus, after several time steps, particles diffuse from the center toward the solid surfaces. At large times, the vertical gradient vanishes an concentration becomes homogeneous across the channel.

Using BB, the diffusion coefficient strongly depends on the vertical size of the channel. Diffusive transport close to the solid surfaces is reduced due to the BB rule, yielding a smaller diffusion coefficient compared to bulk diffusivity. This effect becomes negligible when the gap is large enough or, when the diffusivity is sufficiently small. In Figure 3 we show the dependence of the diffusion coefficient (relative to $`D_m`$), on the size of the gap and bulk diffusivity. We also show that simulations using BF give $`D=D_m`$, for any size and $`D_m`$.

As mentioned, the undesired effects due to BB boundary conditions become less significant at small diffusivities. When using $`D_m=8.8\times 10^4`$ the deviation of the numerical value from the expected one is within $`1\%`$. However, the use of small diffusion coefficients allows large concentration gradients, which may result in numerical instability. Discrepancies between theory and simulations have been reported, in LB model of miscible fluids, for values of $`D_m`$ below $`10^4`$. Also, for values of $`D_m`$ below $`10^4`$, oscillations between negative and positive values of concentration have been observed . Therefore, the parameter space where the model is fairly independent of the particular choice of the boundary condition rule, is very close to being numerically unstable.

### B Taylor dispersion

Finally, we compute the asymptotic hydrodynamic dispersion when a mean flow is set within the channel. In this case, dispersion has two different contributions, one due to molecular diffusion and the other due to spatial variations of the fluid velocity, namely Taylor dispersion . Asymptotic analysis of Taylor dispersion in a two-dimensional channel of constant width $`H`$ and mean flow velocity $`U`$ gives the exact formula for the longitudinal dispersion coefficient,

$$D_{}=D_m+\frac{1}{210}\frac{U^2H^2}{D_m},$$
(2)
Introducing the Péclet number, which accounts for the relative magnitude of the convective and diffusive transport, $`\mathrm{Pe}=\mathrm{UH}/\mathrm{D}_\mathrm{m}`$, the previous equation can be rewritten as,

$$D_{}=D_m\left(1+\frac{\mathrm{Pe}^2}{210}\right).$$
(3)
Concerning the numerical simulations, whether the BF rule is implemented or not, the Taylor regime for the longitudinal dispersion is expected to hold. However, the effective diffusivity is slowed down when BB collision at solid surfaces is used. Thus, as the diffusivity is smaller, the characteristic homogenization time in the transverse direction becomes larger. In other words, we may expect Eq. (2) to hold, but where $`D_m`$ is replaced by the effective diffusion coefficient in the channel $`D`$, as measured in diffusion simulations using BB rule. (section III A). In Figure 4 we compare the numerical results with Eq. (2). There is a very good agreement between the numerical results and the Taylor theory in both cases, using BB and BF as boundary rules. Let us remark that solid and dashed lines corresponds exactly to Eq. (2) with no adjustable parameter, and the only difference between them is the molecular diffusivity value. These results confirms the previous one, showing that the BB rule diminishes the diffusive transport and consequently enlarges the longitudinal hydrodynamic dispersion. A departure from Eq. (2) at large Péclet numbers can also be observed, where the largest discrepancy is $`0.3\%`$. Similar results, showing slightly smaller numerical values than theoretical ones have been previously reported .

As a conclusion to this section we may say that the effect of using the same BC on fluid and color distributions slows down the diffusive transport. This effect becomes appreciable when the system is narrow enough or the diffusion coefficient is large.

## IV Self-affine narrow fractures

We briefly review here the mathematical characterization of self-affinity. A more detailed discussion can be found in our previous work . We consider a rock surface without overhangs, whose height is given 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 (in contrast to self-similar surfaces which stretches the coordinates with equal ratios). Experiment indicates that 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.,

$$z(x,y)=\lambda ^\zeta z(\lambda x,\lambda y)$$
(4)
where $`\zeta `$ is the roughness or Hurst exponent . In all cases the roughness exponent is chosen as the experimentally-observed value $`\zeta =0.8`$.

We shall emphasize the limiting situation of narrow fractures, in which the two surfaces are very close to each other. Consider the situation in which a rock of lateral size $`L`$ is fractured and the two surfaces are simply displaced by a distance $`HL`$, perpendicular to the mean fracture plane, with no relative shift. The fluctuations (the difference between the maximum and minimum values of $`z`$) scales as $`RL^\zeta `$. If $`H`$ is the width of the fracture, then the limiting situation corresponds to $`RH`$.

In narrow fractures, the correlation between the two sides is an important feature for the dispersion process, 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$$
(5)
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`$.

The notion of fractal dimension also differs from the case of self-similar surfaces. It has been shown , that one needs in general several distinct notions. We are interested in the local and global box dimensions of vertical plane sections. Using the self-affine scaling law for the correlation function,

$$\sigma _z^2(\mathrm{\Delta }x)=[z(x,y)z(x+\mathrm{\Delta }x,y)]^2=\varphi (\mathrm{})\left(\frac{\mathrm{\Delta }x}{\mathrm{}}\right)^{2\zeta }$$
(6)
where $`\mathrm{}`$ is a microscopic length, say a grain size, such that

$$\sigma _z^2(\mathrm{})=\varphi (\mathrm{})\mathrm{}^2,$$
(7)
we can estimate,

$$\mathrm{\Delta }Z(\mathrm{\Delta }x)z(x,y)z(x+\mathrm{\Delta }x,y)\mathrm{}\left(\frac{\mathrm{\Delta }x}{\mathrm{}}\right)^\zeta .$$
(8)
For length scales $`\mathrm{\Delta }x\mathrm{}`$ it can be seen from the previous relation that $`\mathrm{\Delta }Z\mathrm{\Delta }x`$ (given that $`\zeta `$ is smaller than one). This so called local property yields a local box counting dimension $`D_{BL}=2\zeta `$. On the other hand, for length scales $`\mathrm{\Delta }x\mathrm{}`$, it is clear that $`\mathrm{\Delta }Z\mathrm{\Delta }x`$. Thus, the global box counting dimension is $`D_{BG}=1`$. As we consider $`\mathrm{}`$ as the lower cut-off of the self-affine behavior in real fractured system, we will focus on the global regime where the box dimension is one. This fact largely determines the convective-diffusive behavior, as we will discuss later.

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 mean flow, when present, 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.

We use statistically self-affine surfaces 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 .

## V Diffusion in narrow fractures

In this section we study diffusion in two-dimensional self-affine fractures. The approach will be based on the tortuosity concept . It is a common characteristic in fluid transport in porous media, and in fractures in particular, that the actual path followed by the fluid is very tortuous. In our fractured two-dimensional systems, a purely geometrical definition of tortuosity can be done, by considering the ratio of the shortest continuous paths between any two points within the fracture to the length of the system projected on to the main plane . It is clear that, for narrow fractures, this ratio will approach the ratio between the length of the surface profile in the $`xz`$ plane to the distance $`\mathrm{\Delta }x`$ between the two points. Then, if $`l_e`$ is the true length between the two points separated a distance $`l=\mathrm{\Delta }x`$, we can write for the tortuosity,

$$T=\frac{l_e}{l}$$
(9)
The tortuosity factor $`T`$ accounts for the reduction in diffusivity, due to irregular geometries, in a very simple way. Along the actual average-path of tracer particles, diffusion behaves as in free bulk, where the diffusion coefficient is the free bulk molecular diffusion coefficient $`D_m`$. On the other hand, diffusive spreading along the $`x`$ axis is reduced, due to the difference between $`l_e`$ and $`l`$. In terms of $`T`$, this difference can be accounted for by writing,

$$D=\frac{l^2}{2t}=\left(\frac{l}{l_e}\right)^2\frac{l_e^2}{2t}=\frac{D_m}{T^2}$$
(10)
(some authors prefer to define tortuosity as $`T=(l_e/l)^2`$, or as the inverse of this definition ).

For the previous equation to be meaningful one should have a constant value of $`T`$, that is, independent of $`l`$. Whether or not $`T`$ is constant clearly depends, for narrow fractures, on the dimension of the vertical plane sections of the surface. In general, the effective length $`l_e`$ depends on $`l`$ through $`l_el^{D_B}`$, where $`D_B`$ is the dimension of the surface profile. Therefore, while at local length scales $`D_B=D_{BL}=2\zeta `$ and $`T`$ depends on $`l`$, at global length scales, considered in this work, $`D_B=D_{BG}=1`$ and the tortuosity factor is constant . For the numerically-generated self-affine surfaces used in this work, we tested the dependence of the effective length $`l_e`$ on $`l`$. Using the method of ‘dividers’  to measure $`l_e`$, we obtain, in all cases, a linear relation between the dividers opening $`ϵ`$ and $`l_e(ϵ)`$, i.e. $`D_B=1`$.

The fact that $`T`$ is constant has some important consequences on diffusive transport. First, the relation between mean square displacement and time should be linear (after perhaps a transient time). Second, the distribution of tracer concentration should be asymptotically Gaussian, and independent of the initial distribution of tracer. Finally, the correction due to tortuosity is strictly geometrical and therefore it should be independent on the actual value of the free bulk molecular diffusivity $`D_m`$.

In Fig. 5 we show the linear dependence of the mean square displacement on time, for different values of the diffusion coefficient. Many other simulations, with different values of $`D_m`$ and $`H`$ were performed, and in all cases a linear relation was obtained. In Fig. 6 we show how the same Gaussian distribution is approached at long times, starting from three different initial distribution of tracer particles. Finally, in Fig. 7 we show that the tortuosity factor is a strictly geometrical property, being clearly independent on the diffusivity of the tracer particles.

### A Tortuosity dependence on the fracture width

In the previous section, we discussed the linear relation between the effective length $`L_e`$ and $`L`$, and its consequences on diffusive transport. We showed that the effect of surface roughness on diffusion can be accounted for by a purely geometrical property of the system, the tortuosity $`T`$. However, it remains to analyze the dependence of $`T`$ on geometrical parameters that describe the system. Of particular importance, is the dependence of $`T`$ on the width of the fracture $`H`$.

The theoretical analysis presented here will closely follow the kind of approximation used in , where the fracture is divided into a sequence of quasi-linear segments at varying orientation angles. First, we estimate a typical size $`\xi _{}`$ in the direction of the mean flow over which the fluctuations in the vertical direction are small compared to the effective aperture of the channel. A segment of length $`\xi _{}`$ is roughly straight when $`\mathrm{\Delta }Z(\xi _{})\mathrm{}(\xi _{}/\mathrm{})^\zeta `$ is a small fraction of $`H`$, which yields

$$\xi _{}\mathrm{}\left(\frac{H}{\mathrm{}}\right)^{1/\zeta }$$
(11)
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`$. Thus, we write for the total length of the channel

$$L_e=\underset{i=1}{\overset{N}{}}\xi _{}^i=\xi _{}\underset{i=1}{\overset{N}{}}\mathrm{cos}^1(\theta _i)$$
(12)
Finally, taking into account 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 tortuosity,

$$T=\frac{L_e}{L}=\mathrm{cos}^1(\theta )$$
(13)
Since $`\zeta <1`$, the channels have small vertical fluctuations, and we can give a simple estimate of the cosine as

$$\frac{1}{\mathrm{cos}\theta }=\frac{\sqrt{\xi _{}^2+\mathrm{\Delta }Z^2(\xi _{})}}{\xi _{}}1+\frac{1}{2}\left(\frac{\mathrm{\Delta }Z(\xi _{})}{\xi _{}}\right)^2$$
(14)
Replacing this approximation in Eq. (13), we obtain,

$$T\left[1+\frac{1}{2}\left(\frac{\sigma _z(\xi _{})}{\xi _{}}\right)^2\right]$$
(15)
where $`\sigma _z^2(\xi _{})=\mathrm{\Delta }Z^2(\xi _{})`$ (see Eqs. (6) and (8)).

A more precise evaluation can be obtained based on a Gaussian distribution of heights , as supported by experimental measurements  and, in fact, the actual distribution given by our numerical procedure for generating self-affine surfaces. In this case, the angular average is given by

$$\mathrm{cos}^1\theta =p(\mathrm{\Delta }Z)\sqrt{1+\left(\frac{\mathrm{\Delta }Z}{\xi _{}}\right)^2}=x^{1/2}U(\frac{1}{2};2;x)$$
(16)
where $`x=(\xi _{}^2/2\sigma _z^2(\xi _{}))`$, and $`U(a;b;x)`$ is the confluent hypergeometric function of the second kind. Therefore, we can get a more convincing evaluation of the tortuosity of a narrow two-dimensional self-affine fracture, using the leading terms in the asymptotic representation of $`U`$ for large $`x`$ and Eq. (11) for the value of $`\xi _{}`$,

$$T=\left[1+\left(\frac{\varphi (\mathrm{})}{\mathrm{}^2}\right)\left(\frac{H}{\mathrm{}}\right)^{\frac{2\zeta 2}{\zeta }}\right]$$
(17)
To compare the previous relation with the numerical results, let us first recast it in terms of the diffusion coefficient,

$$\frac{D_mD}{D_m}2C_1\left(\frac{H}{\mathrm{}}\right)^{\frac{2\zeta 2}{\zeta }}$$
(18)
where we have added an adjustable parameter $`C_1`$, which is expected to be of order one.

In Fig. 8 we present the decrease in the diffusivity due to the tortuosity of the channel as a function of the distance between the opposite surfaces $`H`$. We find a good agreement for the predicted exponent $`(2\zeta 2)/\zeta =1/2`$ (the roughness exponent is $`\zeta =0.8`$). The only adjustable parameter is the coefficient $`C_1`$, which is found to be $`C_10.9`$ in good agreement with the expected value.

In Fig. 9 we present the numerical results obtained for smaller fluctuations in the surface height. In this case it is clear that the exponent differs from the predicted value. We believe that, this discrepancy comes from the discretization of the surface height which introduces a larger error the smaller the fluctuations are. However, we also show in Fig. 9, the theoretical correction to diffusive transport computing the tortuosity factor directly through the numerically-generated surfaces instead of using the asymptotic analysis. In this case a good agreement is recovered and the only adjustable parameter $`C_12`$ is again of order one (and in agreement with previous results ). Similar results were obtained using several intermediate values for the surface fluctuations in height, where the reduction in diffusivity is always in agreement with the correction due to the tortuosity of the channels. In addition, we found that the decrease in diffusivity is, in all cases, well described by a power-law (where the power-law exponent increases slightly with the amplitude of the surface fluctuations), for which we have no explanation.

## VI Dispersion in narrow fractures

In this section we address the case of Taylor hydrodynamic dispersion in narrow fractures. First, we will analyze the case where the two sides of the fracture are displaced normally to the mean fracture plane, and then we will turn to the case where there is a lateral shift as well.

When the two complementary surfaces are simply displaced vertically by a distance $`HL`$, the vertical aperture of the system is constant everywhere. Nevertheless, the flow field differs from one in a straight channel, due to variations in the local width of the channel normal to the mean flow direction. In Fig. 10 we show a set of streamlines inside a fracture, where the effect of the varying effective aperture of the fracture is evident (aperture normal to the mean flow). Moreover, in  we show that the complex geometry of the fracture gives rise to low-velocity zones (close to depressions and corners), reducing the permeability of the system. In order to describe the dispersion process we need to obtain a measure of the fraction of the system that is subject to convection. To this end we will make use of the $`\mathrm{\Lambda }`$ parameter , which is directly related to transport, and measures the dynamically connected part of the pore space in porous media . Following Ref. , we write $`\mathrm{\Lambda }`$ in terms of the permeability and tortuosity of the system,

$$k=\frac{\mathrm{\Lambda }^2}{12T^2}$$
(19)
Note that defining a characteristic length, as the ratio between pore volume $`V`$ and surface area $`S`$ would yield $`V/S=H/T\mathrm{\Lambda }`$, which does not depend on the effects exerted on the fluid flow by the complex geometry of the fracture.

Now, assuming that we have a winding channel of length $`L_e`$, effective aperture for convective transport $`\mathrm{\Lambda }`$ and actual aperture $`A=H/T`$ (given by volume conservation), we may apply the same reasoning as in Section III B to get the dispersion coefficient. After replacing $`H`$ by the width of the effective channel $`H/T`$, $`U`$ by the mean velocity in the convective part of the channel $`UH/\mathrm{\Lambda }`$, and taking into account the tortuosity, we obtain

$$D_{}=\frac{D_m}{T^2}\left(1+\frac{Pe^2}{210}\left(\frac{H}{\mathrm{\Lambda }T}\right)^2\right)$$
(20)
Thus, to see how well the $`\mathrm{\Lambda }`$ parameter can be used to estimate the dispersion, we will compare the values obtained both from dispersion measurements (Eq. (20)) and from the permeability (Eq. (19)).

In Figure 11 we show the numerical results obtained for the dispersion coefficient when varying the injection rate. The linear behavior shows that, as expected, Taylor dispersion is governing tracer spreading. The fracture gap is $`H=16`$, the length of the system is $`L=512`$, and the relaxation parameter used is $`\lambda _D=1.9`$. The molecular diffusivity and the lambda parameters obtained from the the best fit are: $`D_m=(8.86\pm 0.01)\times 10^3`$ which differs by only $`1\%`$ from the theoretical value; and $`\mathrm{\Lambda }=13.4\pm 0.1`$. On the other hand, from the flow rate computed in numerical simulations, and estimating the permeability by Eq. (19), we get $`\mathrm{\Lambda }=15.3\pm 0.2`$.

We now turn to shifted surfaces. In the presence of a small lateral shift $`d`$ between complimentary surfaces, most of the previous discussion remains valid. The local aperture now varies with position and $`H`$ is the average aperture of the channel. The difference in height between surfaces at any point $`x`$ along the channel is of order $`d^\zeta `$, while the vertical excursion of the fracture between points separated a distance $`\mathrm{\Delta }x`$ is $`(\mathrm{\Delta }x)^\zeta `$. Therefore, at large length scales, $`\mathrm{\Delta }xd`$, the fracture may be considered as a winding channel of length $`L_e`$ and effective aperture $`\mathrm{\Lambda }`$ and where the ratio between the length of the channel $`L_e`$ and the system length $`L`$ is the tortuosity factor.

In Figure 12 we show the dispersion coefficient, as a function of the Péclet number. The two fracture surfaces are vertically displaced by $`H=16`$ and laterally shifted by $`d=8`$. As expected, we observed a linear dependence, in agreement with Taylor-like spreading. From the best fit of the numerical results we obtain $`D_m=(8.5\pm 0.1)\times 10^3`$ which differs less that $`3\%`$ from the theoretical one ($`\lambda _D=1.9`$); and $`\mathrm{\Lambda }=14.6\pm 0.2`$. Whereas, by means of flow rate data computed in numerical simulations and Eq. (19) for the permeability, we obtain $`\mathrm{\Lambda }=11.7\pm 0.5`$.

Even though estimated values of $`\mathrm{\Lambda }`$ are in fairly good agreement (within a $`20\%`$ discrepancy), it is also clear that the $`\mathrm{\Lambda }`$-parameter fails to completely predict the enhancement of dispersion due to the complex geometry of the fractures. Nevertheless, we believe that the presence of low-velocity zones, is the only possible feature that accounts for the enhancement in the spreading of tracer particles.

Let us note that the uncertainty in the computed dispersion coefficient is considerably larger than in the case where there is no lateral shift. It has been shown, in three-dimensional fractures and under lubrication approximation for the velocity field, that a lateral shift yields geometric dispersion for tracer particles advected along the flow . This effect is due to different mean velocity along different streamlines. In our case, the two-dimensional nature of the fractures prevents the presence of geometric dispersion, given that the height-averaged velocity is constant throughout the system. However, an analogous effect appears when averaging over different fractures, giving rise to a large dispersion in the computed mean velocity and dispersion coefficient. In fact, from the previous discussion an uncertainty proportional to the mean velocity is expected when averaging the dispersion coefficient (corresponding to a geometric dispersion term in fully three-dimensional fractures).

In order to validate the dispersion results obtained by means of the LB method, we shall now compare them with results computed via a Monte-Carlo (MC) approach to dispersion process. In the MC method, one follows the displacement of a large number of particles, or random walkers, moving in a two-dimensional fracture. The motion of each particle is a combination of the effects of molecular diffusion and convection (we assume, as in the LB simulations that tracer particles moves independently of each other). In time $`\mathrm{\Delta }t`$, a particle is displaced according to

$$\mathrm{\Delta }𝐱=𝐮(𝐱)\mathrm{\Delta }t+\widehat{𝐧}(4D_m\mathrm{\Delta }t)^{1/2}$$
(21)
where $`𝐮(𝐱)`$ is the velocity field obtained by the LB method, $`\widehat{𝐧}`$ is a unit vector with random orientation, and the amplitude of the random steps has been chosen so that the variance of any coordinate is $`2D_m\mathrm{\Delta }t`$. Boundary condition at solid surfaces are implemented as in , where those random steps that would take the particle outside the channel are suppressed. The sequence of steps is repeated while recording the distance $`\mathrm{\Delta }x`$ from the initial position of the particle. The process is repeated for a large number of particles and average values are computed.

We found an excellent agreement between the two methods. In Fig. 13 we compare the mean square displacement as a function of time obtained using LB and MC methods. The agreement is evident, and it can be seen that after releasing $`10^5`$ particles in the MC simulations, the noise is negligible. Both simulations corresponds to a fracture of mean aperture $`H=16`$, length $`L=512`$, and a lateral shift between surfaces $`d=8`$. It is also interesting to note, in Fig. 13, a change in the slope at time $`t20000`$, which approximately corresponds to the characteristic time for transverse diffusion across the aperture $`\tau _DH^2/2D_m15000`$. Therefore, this marked change in slope is showing the transition towards Taylor-like spreading at times larger than $`\tau _D`$.

### A Dispersion at small Péclet numbers

In a straight channel, the presence of convection increases the dispersion of tracer particles. The same effect is usually found in porous media, including fracture systems. However, in the case where surfaces are laterally shifted, and at very small Péclet number, we observed the opposite behavior. That is, the presence of convective transport inside the fractures reduces the dispersion of tracer particles. This effect can be observed in Fig. 12, at very small Pe the dispersion coefficient grows as the flow rate decreases. In order to validate this observation we simulate the dispersion process with two other methods. The first one is the MC method presented before. The second method is a variation of the MC where, instead of using the velocity field computed by means of LB simulations, we use lubrication approximation for velocities. In the lubrication approximation the velocity $`u(x,z)`$ is given by,

$$u(x,z)=\frac{6Q}{H^3(x)}z(H(x)z)$$
(22)
where $`H(x)`$ is the local aperture of the fracture.

In Figure 14 we show the results obtained using the three different methods in a range of Pe from $`0`$ to $`4`$. The agreement between methods is excellent and all show an initial decrease in dispersion due to the presence of weak convective transport. We believe that, at very low velocities, the main effect of convective transport is to carry out tracer from low-velocity zones, and thus, reducing dispersion. This explanation is supported by the fact that lubrication approximation yields very similar results. (Note that this approximation assumes local Poiseuille profile without accounting for possible stagnant zones). Unfortunately, the small Péclet numbers at which this surprising effect arises, makes it very difficult to be experimentally measured.

## VII Conclusions

In first place, we have presented a new set of boundary conditions to describe diffusive transport within the lattice-Boltzmann method. We tested this new boundary rule in two situations, simple diffusion and Taylor dispersion in two-dimensional straight channels. We showed that the proposed bounce-forward rule improves the accuracy of the method and does not possess the undesirable effects of the bounce-back rule, i.e., dependence of the diffusion coefficient on the aspect ratio of a straight channel and transient concentration gradients near solid walls. Even though the accuracy might be recovered in the BB case using small diffusivities, we showed that this option leads the simulations towards the numerical instability border.

We then turn to study diffusive transport in two-dimensional self-affine fractures. First, we showed that the slow-down in diffusive transport can be accounted for by the purely-geometric tortuosity factor. Our numerical simulations have verified all the implications of this result, that is, linear spreading on time, Gaussian distribution of tracer and independence of geometrical effects on the actual value of the diffusivity. Secondly, using analytic arguments in the limit of small aperture fluctuations, we have obtained an expression for the tortuosity in terms of the fracture gap and the Hurst exponent characterizing the fracture surface. Numerical simulations verify the validity of the theoretical approach even when the discrete nature of the surface seems to affect the asymptotic scaling law.

Finally we studied tracer dispersion in fractures in the presence or not of a lateral shift between complementary surfaces. In both cases we showed that tracer spreading can be described as analogous to Taylor-like dispersion in a straight channel. In the case without lateral shift we showed that, the enhancement in dispersion can be mostly understood assuming an effectively reduced aperture for fluid transport, due to the rugosity of the surface. The $`\mathrm{\Lambda }`$-parameter measuring this effective aperture, and computed from permeability measurements, was shown to be similar to that estimated from dispersion measurements. We also presented a novel result, showing a decrease in tracer dispersion when convective transport is set in the fracture. This last observation was obtained in the frame of three different numerical approaches.

## Acknowledgments

We thank J. P. Hulin, F. Plouraboué for many fruitful discussions and M. Tanksley, N. Rakotomalala and D. H. Rothman for useful discussions on the lattice-Boltzmann method. G. D. thanks M. Tirumkudulu and I. Baryshev for their helpful comments. G. D. was partially supported by CONICET Argentina and The University of Buenos Aires. 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:
# Nature and Nurture: a model for soft gamma-ray repeaters

## 1. Introduction

Soft Gamma-ray Repeaters (hereafter SGRs) and Anomalous X-ray Pulsars (hereafter AXPs) are two groups of enigmatic sources. They share the following properties: 1. They all have long rotation periods (clustered within 5-12 seconds) and large spin-down rates (see, e.g. Mereghetti & Stella 1995; Kouveliotou et al. 1998, 1999); 2. Most of them are associated with supernova remnants, indicating that they are young objects (for reviews, see Hurley 1999; Mereghetti 1999); 3. No optical, infrared or radio counterparts have been identified (e.g. Eikenberry & Dror 2000; Lorimer & Xilouris 2000); 4. They all have soft persistent pulsed X-ray emission with luminosities of $`L_\mathrm{x}10^{35}10^{36}\mathrm{ergs}\mathrm{s}^1`$, well in excess of the spin down energy of these sources (e.g. Thompson 2000 for a review). The main difference between both types of the objects is that SGRs show occasional soft gamma-ray bursts while AXPs do not. It is also found that SGRs usually have larger proper motion velocities than AXPs according to their relative positions with respect to the cores of their supernova remnants (Hurley 1999). The main characteristics of the SGR bursts include: 1. Most of the bursts have super-Eddington luminosities with $`L_\mathrm{b}10^{38}10^{42}\mathrm{erg}\mathrm{s}^1`$; 2. The fluence distribution of the bursts is a power-law, and there is no correlation between the burst intensity and the time intervals between the bursts (Gögüs et al. 1999; 2000); 3. Two giant flares have been detected from SGR 0526-66 (the March 5, 1979 event) and SGR 1900+14 (the August 27, 1998 event), which share some common properties (see Thompson 2000 for a review); 4. Most bursts have soft spectra with characteristic energy around 20-30 keV.

The popular model for SGRs and AXPs is the magnetar model, which can account for almost all the phenomena listed above (Duncan & Thompson 1992; Thompson & Duncan 1995, 1996; Thompson 2000). However, the differences between SGRs and AXPs are not straightforwardly interpreted since these objects are not intrinsically different objects within the magnetar picture. It also remains unclear how some other issues, e.g., the non-systematic discrepancy between the characteristic ages derived assuming dipolar spindown and the ages of the associated supernova remnants, no clear positive dependence between $`L_\mathrm{x}`$ and the polar surface field strength $`B_\mathrm{p}`$, etc., can be properly addressed. On the other hand, a fossil-disk-accretion model for AXPs recently emerges from the independent studies by Chatterjee et al. (Chatterjee, Hernquist & Narayan 2000; Chatterjee & Hernquist 2000) and Alpar (1999, 2000). The neutron stars in such a scenario have normal magnetic fields as the Crab pulsar. The model can interpret the AXP phenomenology well, but the bursts from the SGRs are difficult to interpret. On the observational ground, Marsden et al. (2000) observed that the SGRs and the AXPs are located in a much denser environment than the normal pulsars. They hence argue that the peculiar behaviors of the SGRs and AXPs may be due to their “nurture” from the environment rather than due to their special “nature” (i.e. magnetars) as compared with the normal pulsars. However, no plausible idea was proposed to connect the “nurture” to the phenomenology of these sources, especially the bursting behavior of the SGRs.

In this Letter, we attempt to propose a model to understand the bursting behavior of the SGRs without introducing the magnetar idea. We propose that the central objects of the SGRs are “bare” strange stars with normal magnetic fields ($`10^{12}10^{13}`$ G). We assume that these strange stars are born directly from supernova explosions from some massive progenitors, and they have experienced a spindown history as that having been proposed for the AXPs within the fossil disk model (Chatterjee et al. 2000; Alpar 2000). According to this model, some fallback materials from the supernova ejecta will form a fossil disk around the strange star. The SGRs/AXPs are just such strange/neutron stars that have experienced the “propeller” phase ($`r_\mathrm{c}r_\mathrm{m}<r_\mathrm{l}`$), and are now in the “tracking” phase ($`r_\mathrm{c}r_\mathrm{m}<r_\mathrm{l}`$) when infall of the materials onto the surface is possible and the star is X-ray bright. Here $`r_\mathrm{l}`$, $`r_\mathrm{m}`$ and $`r_\mathrm{c}`$ are the light cylinder, the magnetospheric radius, and the corotating radius, respectively. In our picture, AXPs may be still neutron stars. We will attribute the SGR bursts to their occasionally collisions with some comet-like objects in the dense environment of the SGRs. We will show how various SGR properties as reviewed above could be accounted for within this picture. Our model differs from some other strange star SGR models (e.g. Alcock et al. 1986b; Cheng & Dai 1998; Dar & de Rujula 2000).

## 2. The model

Strange stars (Haensel, Zdunik & Schaeffer 1986; Alcock, Farhi & Olinto 1986a) are hypothetical objects based upon the assumption that strange quark matter is more stable than nuclear matter (Witten 1984; Farhi & Jaffe 1984). Though the existence of such stars are still subject to debate, some evidence in favor of strange stars has recently been collected (e.g. Li et al. 1999a, 1999b; Titarchuk & Osherovich 2000). Strange stars can be either bare or have normal matter crusts (Alcock et al. 1986a). They can be formed directly during or shortly after some supernova explosions when the central density of the proto-neutron stars is high enough to induce phase conversion (e.g. Dai, Peng & Lu 1995; Xu, Zhang & Qiao 2000). If a strange star is born directly from a supernova explosion, it is likely that the star might be almost bare (Xu et al. 2000). Some radio pulsars may be such strange stars with exposed bare quark surfaces (Xu, Qiao & Zhang 1999).

There are three main motivations for us to choose (bare) strange stars rather than neutron stars to interpret the SGRs. 1. A prominent feature of the SGR bursts is their super-Eddington luminosities. This feature has been regarded as a strong support to the magnetar model, since superstrong magnetic fields may considerably suppress the Thompson cross section and consequently raise the Eddington limit to several orders of magnitude higher (Paczynski 1992; Thompson & Duncan 1995). However, the luminosities of the most luminous events, e.g., the initial spike of the March 5 event with $`L10^{44}\mathrm{ergs}\mathrm{s}^1`$, are still above the enhanced Eddington limit. An important merit of bare strange stars is that they are not subject to Eddington limit at all since the bulk of the star (including the surface) is bound via strong interaction rather than gravity (Alcock et al. 1986a). This presents a clean interpretation to the super-Eddington luminosities of the SGRs, as long as the impacts are not in the polar cap region where the accretion flow from the fossil accretion disk is channeled. 2. As criticized by Thompson & Duncan (1995), the impacting model for neutron stars suffers the baryon contamination problem. The impact may load too much baryonic matter to cause adiabatic dilution of photons in an expanding fireball to energies well below the hard X-ray and $`\gamma `$-ray band. A bare strange star can naturally evade such a criticism, since the infall matter will be essentially converted into strange quark matter within a very short period of time ($`10^7`$ s, Dai et al. 1995) when they penetrate into the star. A new-born bare strange star may have a very thin normal matter atmosphere (Xu et al. 2000), which is far less than the amount required to pollute the fireball. 3. Observationally, SGRs tend to have larger proper motion velocities ($`1000\mathrm{k}\mathrm{m}\mathrm{s}^1`$) than normal pulsars and AXPs. Though we do not attempt to propose a detailed “kick” theory in the present Letter, we note that the formation of a strange star rather than a neutron star may potentially pose some possibilities to interpret the large proper motion velocities of SGRs. Present kick theories invoke either hydrodynamically-driven or neutrino-driven mechanisms (Lai 2000). For the former, the kick arises from presupernova g-mode perturbations amplified during the core collapse, leading to asymmetric explosion (Lai & Goldreich 2000). We note that the formation of a strange star is a two-step process, i.e., the formation of a proto-neutron star and phase conversion. Neutrino emission in the second step could be significantly asymmetric since the phase conversion may be off-centered due to the initial density perturbation (Dong Lai, 2000, personal communication). An off-centered transition condition may be also realized in the presence of an electron-neutrino-degenerate gas in a proto-neutron star (Benvenuto & Lugones 1999). Thus the phase transition process may give an additional kick to achieve a higher velocity. More detailed investigations are desirable to verify these proposals.

We now describe the model in more detail. We assume that the progenitor of a strange star is surrounded by a huge spherical comet cloud which is similar to the Oort Cloud in the solar system. They may be formed during the formation of the massive star, and have almost finished gravitational relaxation. Since the progenitor of a strange star should have a mass larger than $`10M_{}`$, we expect that the radius of the Oort Cloud in the progenitor system may be one order of magnitude larger than the solar value ($`2\times 10^{13}`$ km Weissman 1990)., i.e., $`r_\mathrm{o}2\times 10^{14}`$ km. Supernova explosion blast waves will not destroy these comet clouds (Tremaine & Zytkow 1986). The luminous UV/optical emission from the progenitor is also unlikely to evaporate the comets. Although the radiation flux received by the Oort Cloud comets of the massive star should be about a factor of 30 higher than that received by the Solar Oort Cloud comets, the existence of copious “Kuiper Belt” comets in the solar system (which is 4 orders of magnitude closer to the sun than the Oort Cloud) hints that comets can withstand shining with much higher luminosities. The influence of nearby stars may be also not prominent due to the same reason, even if SGRs are associated with luminous star clusters (e.g. Vrba et al. 2000). Using the typical proper motion velocity of the SGRs, $`V_{_{\mathrm{SGR}}}10^3\mathrm{km}\mathrm{s}^1`$, and the typical supernova remnant age, $`t_{_{\mathrm{SGR}}}10^4`$ yr, the distance that a SGR has traveled since its birth is $`r3\times 10^{14}`$ km, remarkably consistent with the distance of the Oort Cloud $`r_\mathrm{o}`$. Thus the age clustering of the SGRs near $`10^4`$ yr is simply due to that this is the age when a lot of impacts are available. The lack of bursts from the AXPs may be due to their much smaller proper motion velocities, and probably also their different nature, i.e., neutron stars. Although SGR 1806-20 has a smaller projected proper motion velocity ($`V_{}100\mathrm{k}\mathrm{m}\mathrm{s}^1`$), we assume that it has a similar velocity as other SGRs, with a large velocity component along the direction of the line-of-sight. The capturing rate could be estimated as $`\dot{N}\pi (2GM_{}/V_{_{\mathrm{SGR}}}^2)^2V_{_{\mathrm{SGR}}}n_\mathrm{c}`$, where $`n_\mathrm{c}`$ is the number density of the comets within the Oort Cloud. To produce a bursting rate of 1 yr-1, $`n_\mathrm{c}`$ is required to be $`(10^{22}10^{23})\mathrm{km}^3`$. This is about 4 orders of magnitude higher than the inferred comet number density in the solar Oort Cloud \[$`(10^{26}10^{27})\mathrm{km}^3`$, Weissman 1990\], but about 3-4 orders of magnitude lower than the inferred number density in the Kuiper Belt of solar system \[$`(10^{18}10^{20})\mathrm{km}^3`$, Weissman 1990\]. Keeping in mind that the mass density of the Oort Cloud and the number density of the comets may be enhanced due to accretion from the dense environment in the supernova remnants (Marsden et al. 2000) and that the number density quoted for the solar system might be a lower limit (Weissman 1990), the required $`n_\mathrm{c}`$ may be not unreasonable. Some SGRs have more frequent bursting rate. This may be due to that the strange star has captured a denser small-scale comet cloud.

When the strange star passes through its Oort Cloud, it may capture some small-scale clouds and make them circulate around it within its rest frame111The fossil disk around the star may be a good perturber of the comets, which enhances the chances of captures., and the comets within the cloud will be occasionally accreted onto the strange star surface. The different bursting luminosities (or more precisely the different energies for different bursts) correspond to different masses of the impacting objects. During each impact, the energy released is a sum of the gravitation energy and the phase conversion energy. The former has an efficiency of $`\eta _{\mathrm{grav}}=GM/(Rc^2)`$, which is $`20.6\%`$ for typical strange star parameters, and the latter has an efficiency of $`\eta _{\mathrm{conv}}=\mathrm{\Delta }ϵ/(930\mathrm{M}\mathrm{e}\mathrm{V})`$, where $`\mathrm{\Delta }ϵ`$ is the energy per baryon released during the phase conversion. The value of $`\mathrm{\Delta }ϵ`$ is rather uncertain which depends on unknown QCD parameters (e.g. MIT bag constant, strange quark mass and the coupling constant for strong interaction). Some recent calculations (e.g. Bombaci & Datta 2000) show that $`\mathrm{\Delta }ϵ100`$ MeV may be reasonable, and we will adopt this value for indicative purpose. The deviation of this value from the exact value is not important since this only reflects slightly different required comet masses. We thus get $`\eta _{\mathrm{conv}}11\%`$. Assuming that about one half of the energy will be brought away by neutrinos, the total $`\gamma `$-ray emission efficiency is $`\eta _\gamma (\eta _{\mathrm{grav}}+\eta _{\mathrm{conv}})/216\%`$. Thus the repeating bursts with $`L_b10^{38}10^{42}\mathrm{erg}\mathrm{s}^1`$ and typical bursting time $`0.1`$ s correspond to the comet masses within the range of $`7\times (10^{16}10^{20})`$ g. These are reasonable values for comet masses. The so-called giant flare requires an object (an asteroid or a comet) with a mass of several $`10^{24}`$ g. In view that the giant flares are rather rare, it is reasonable to suppose that such large objects may exist in some dense clouds. Notice that all the luminosities quoted above are derived under the assumption of isotropic emission. For impacting events discussed here, during which the emission is anisotropic, the required comet masses may be lowered by a factor of 10-100. There is no mass distribution data available for the solar comets, but we expect that the distribution should be a power law (cf. Pineault & Poisson 1989). This is because the stars, which also belong to a gravitationally self-organized system but in a larger scale, have a well-known Salpeter’s power law mass distribution222Observationally the giant flares belong to the high end of the power-law fluence distributions.. The bursting intervals depend on the spatial distribution of the comets within their orbits, thus there should be no correlations between the luminosity of a burst and the waiting time before or after this burst. All these are in excellent agreement with the statistics of the SGR bursts (Gügüs et al. 1999, 2000). Adopting the typical comet mass as the lowest value of the power-law distribution, the comet number density inferred above gives a total comet mass of about $`0.1M_{}`$, not unreasonable due to the same reasons discussed before.

When a comet falls into the strange star magnetosphere, it will endure tidal distortion and compression so that they are elongated dense solid objects when they reach the strange star surface (Colgate & Petschek 1981). Because they are globally neutral solid bodies, these comets will not be channeled to the polar cap regions where the asymptotic accretion flow from the fossil-disk takes place. This ensures the super-Eddington luminosity emission from a bare strange star. The large Coulomb barrier above the bare quark surface (Alcock et al. 1986a) will not prevent the object from penetrating into the quark core. The rising rate of the energy released from the falling object is similar to the rising rate of the density from vacuum to solid iron (Howard, Wilson & Barton 1981; Katz, Toole & Unrul 1994), so that the rising time of the bursts could be of sub-millisecond to millisecond order, consistent with the observations of the giant flares (Hurley et al. 1999). The duration of the hard spike observed in the giant flares corresponds to the continue infall time of the object, which is of the order of 0.1-1 s (e.g. Katz et al. 1994). The August 27 giant flare from SGR 1900+14 has slightly smaller total energy but both longer rising time and longer duration of the initial spike than the March 5 event of SGR 0526-66. This may be understood by assuming that the falling object of the March 5 event is an asteroid while that of the August 27 event is a comet, both with a similar mass. During an impact, both gravitational energy and phase transition energy will be released in a sufficiently short period of time. Since there is no baryon contamination for a bare strange star, the energy will be mainly released as photons and neutrinos. Soon an optically thick pair fireball will form via the processes such as $`\gamma \gamma `$ (Thompson & Duncan 1995) and $`\gamma E`$ (Usov 1998) processes near a bare quark surface. The magnetic field will confine this pair plasma, and the soft fading tail of the giant flares can be due to contraction of this pair bubble (Thompson & Duncan 1995; Katz 1996). For the accretion case discussed here, the energy deposited into the pair bubble is continually supplied, which is different from the abrupt-release case in the magnetar model. Thus the required magnetic field for confinement is less demanding, i.e. $`B>(2L_\mathrm{b}/R^2c)^{1/2}=8\times 10^{10}\mathrm{G}L_{44}^{1/2}R_6^2`$ (Katz 1996). The trapped pair plasma has a characteristic temperature of $`T23\mathrm{k}\mathrm{e}\mathrm{V}`$, and the emergent spectrum is roughly a blackbody with absorption, which is almost independent on the size of the impacting object (Katz 1996). All these match the SGR phenomenology well.

Sometimes the accreting matter is not solid, but is an ionized plasma. In such cases, the effect of the large Coulomb barrier should be carefully investigated. The kinetic energy of a proton when it is accreted onto the strange star surface is $`E_\mathrm{k}Gm_\mathrm{p}M/2R100`$ MeV. However, when materials are accreted as fluid, it is possible that the kinetic energy will be radiated away via heat before hitting the surface. In the accretion column, the scale of the shock wave zone is dependent on the accretion rate. It is found that when the accretion luminosity is less than $`4\times 10^{36}\mathrm{erg}\mathrm{s}^1`$, the deceleration of the accreting fluid can be neglected (Basko & Sunyaev 1976). This is true for SGRs and AXPs since the quiescent X-ray luminosities of these objects are only $`10^{35}10^{36}\mathrm{erg}\mathrm{s}^1`$. The Coulomb barrier of a bare strange star is $`E_\mathrm{C}=(3/4)V_\mathrm{q}15`$ MeV, where $`V_\mathrm{q}^3/3\pi ^220`$ MeV is defined as the quark charge density inside the quark matter (Alcock et al. 1986a). Thus the accreting fluid, including that from the fossil-disk, can also penetrate into the strange quark core. This ensures the bare strange star picture conjectured in this paper. The accreted matter at the polar cap will undergo phase transition and release some extra energy. It is unclear whether the slightly harder spectra of the quiescent emission of the SGRs with respect to the AXPs is caused by phase transition (we suppose AXPs to be neutron stars). The enigmatic precursor of the August 29 event of SGR 1900+14 (Ibrahim et al. 2000) may be due to infall of an extended ionized cloud which is followed by a solid object.

Depending on the impacting angles during the captures, the small-scale comet clouds may have various orbital periods and eccentricities, so that the precipitation onto the star surface is expected to be periodic, especially when the comets are clustered into a clump in the orbit rather than being spread over the orbit. In fact, SGR 0526-66 has been reported to have a 164-day period in bursts (Rothschild & Lingenfelter 1984). Its present quiescence may be because the previous comet cloud has been depleted due to many cycles of precipitations. If it becomes active again, a different period is expected since it may have captured a different cloud. SGR 1900+14, on the other hand, has experienced three active periods during 1979 (Mazets et al. 1981), Jun. - Aug. 1992 (Kouveliotou et al. 1993), and May 1998 - Jan. 1999 (Gögüs et al. 1999). The active periods are short, and the interval between the first two is roughly twice of that between the last two (about 6 year). This makes us to suspect a 6-year period for SGR 1900+14 activity. According to this picture, there should be some bursts in 1986. But this is within the “detection gap” of the SGR bursts when there is no gamma-ray mission in space before BATSE was launched. Thus we expect that SGR 1900+14 should become active again during 2004-2005. This will give a definite test to our model. SGR 1806-20 activity does not have a clear periodicity. However, a plausible 733-day period is found from its timing residual (Woods et al. 2000). This might be due to that comets are almost spread over the whole orbit and the spin-down of the strange star is perturbed by this comet orbit.

In our model, a fossil-disk is assumed to interpret the spin-down behavior and the quiescent emission. It is expected that emission (especially during the bursts) should have some interactions with the disk with certain optimal geometric configurations. The chance to see such interactions should be small due to the small size of the disk. The 6.4 keV emission line from the August 19 burst of SGR 1900+14 (Strohmayer & Ibrahim 2000) may be due to the disk’s re-processing the bursting emission.

## 3. Discussions

In this Letter, we propose that the peculiar behaviors of the SGRs are due to both their Nature (bare strange stars) and their Nurture (the Oort Cloud in the dense environment). Instead of invoking the magnetar hypothesis, we adopt the strange star hypothesis to interpret some interesting features of the SGRs. It is worth pointing out that the periodic activity does not depend on the nature of the central star. Although some authors argue that the bursting phenomenology (e.g. super-Eddington luminosity) can be also interpreted by colliding comets with a neutron star (e.g. Katz 1996), we think that a bare strange star is a cleaner interpretation due to the reasons discussed above. An important criterion to differentiate our model from the magnetar model is the activity period. If SGR 1900+14 will be turned on again in 2004, the magnetar model is then not favored, since it may be hard to find a mechanism to trigger the magnetic field decay instabilities periodically. If it turns out that some problems (e.g. super-Eddington luminosity, baryonic contamination, and large proper motion velocity) are not solvable within the neutron star impacting model, the bursts from SGR 1900+14 in time then present a support to the strange star hypothesis and will bring profound implications for fundamental physics.

According to this picture, there might be some other bare strange stars which may also have super-Eddington bursts when they collide with comet-like objects. However, they must have passed through the Oort Cloud and/or in a much less dense environment, so that the chance to detect repeating bursts is rare. Single bursting events are possible and they may account for a small portion of the short, soft bursts in BASTE data. The association of SGR 1806-20 with a radio plerion may not be compatible with the present picture, but recent results indicate that the non-thermal radio core of the supernova remnant G10.0-0.3 may be associated with another luminous blue variable rather than with the SGR (Hurley 1999).

BZ acknowledges stimulated discussions with Matthew Baring, Zigao Dai, Alice Harding, Dong Lai, David Marsden and Vladimir Usov, and informative email contacts with Pinaki Chatterjee, Chryssa Kouveliotou, Tan Lu and Chris Thompson. RXX thanks supports from NNSF of China (19803001). GJQ acknowledges supports from NNSF of China, the Climbing Project of China, and the Research Fund for the Doctoral Program Higher Education.


text:
# Zigzag equilibrium structure in monatomic wires

\[

## Abstract

We have applied first-principles density-functional calculations to the study of the energetics, and the elastic and electronic properties of monatomic wires of Au, Cu, K, and Ca in linear and a planar-zigzag geometries. For Cu and Au wires, the zigzag distortion is favorable even when the linear wire is stretched, but this is not observed for K and Ca wires. In all the cases, the equilibrium structure is an equilateral zigzag (bond angle of 60o). Only in the case of Au, the zigzag geometry can also be stabilized for an intermediate bond angle of 131o. The relationship between the bond and wire lengths is qualitatively different for the metallic (Au, Cu and, K) and semiconducting (Ca) wires.

\]

Very recently, using a combination of transmission electron microscopy (TEM) and scanning tunneling microscopy, Ohnishi et al.  were able to visualize the formation of monatomic gold wires. In particular, one of their images shows a wire of four gold atoms forming a bridge between two gold tips. The atoms in the wire were spaced at intervals of 3.5-4.0 Å, a distance much larger than that of nearest-neighbors in bulk gold (2.88 Å). In spite of these large interatomic distances, the wire was observed for more than two minutes, showing clearly its stability. Observation of gold monatomic chains with a length of four or more atoms was independently claimed by other experimental group . In this work, Yanson et al. do not provide direct structural information, but the very long and flexible contacts formed in the last plateau of the conductance curves (with a conductance very close to one quantum, 2e2/h) were associated to the presence of monatomic chains. These new truly one dimensional structures open a great number of opportunities to test our understanding of the interactions in solids. However, the very large interatomic distances found in the experiment represent by themselves a puzzle. In ref.  we already addressed this issue with the aid of first-principles density-functional calculations. Our estimate of the maximum interatomic distance that the gold wires can hold before breaking was $``$2.9 Å, in clear disagreement with the commented experimental observations. Several recently published references also report similar conclusions . The cause for such a discrepancy between theory and experiment was not clear. However, in ref.  we proposed a simple explanation for the puzzle posed by the TEM experiments. First, the gold wires were found to exhibit an equilibrium zigzag structure. Second, the calculated rotational energy barrier was low enough for the small zigzag wires, supported between gold tips, to be rotating at room temperature. In these conditions, only half of the atoms would be clearly visualized, the other half would produce a fuzzy image, and the apparent interatomic distance observed would correspond to that between second-neighbors.

Zigzag structures for monatomic wires have already been observed and proposed by several authors. Whitman et al.  observed the formation of long, isolated zigzag chains when Cs was deposited over GaAs(110) and InSb(110) surfaces. March and Rubio  also proposed zigzag monatomic chains as the building blocks for highly expanded liquid Rb and Cs. A zigzag distortion was also predicted by Chaplik  for a chain of electrons held on the surface of liquid helium by the field of a charged filament.

In this paper we use first-principles density-functional calculations to study the stability and electronic properties of linear and zigzag wires, for Cu, K, and Ca, comparing the results with those of Au . Our standard density functional theory (DFT)  calculations have been performed with the SIESTA program . Core electrons are replaced by scalar-relativistic norm-conserving pseudopotentials , plus a pseudocore to account for non-linear core corrections in the exchange-correlation energy . We used a doble-$`\zeta `$ polarized basis set of numerical atomic orbitals, obtained from the solution of the atomic pseudopotentials with an excitation energy of 50 meV . Real and reciprocal-space integrations were performed with a 200 Ry-cutoff grid and with 20 inequivalent k-points, which guarantee a convergence better than 1 meV/atom. We present results obtained with the local density approximation (LDA)  for Au and Cu, and with the generalized gradient approximation (GGA)  for Ca and K. However, we stress that all the results are very similar using LDA and GGA. We have performed test calculations for the bulk phases of the studied materials, obtaining results in good agreement with experiment and with previous DFT calculations.*** The calculated equilibrium lattice parameter, bulk modulus and, cohesive energy are: 4.11 Å, 194 GPa and, 4.40 eV/atom for Au, 3.56 Å, 210 GPa and, 4.88 eV/atom for Cu, 5.51 Å, 20.3 GPa and, 2.20 eV/atom for Ca, and 5.29 Å, 4.3 GPa and, 0.94 eV/atom for K. In order to minimize the interactions between the wires and their periodic images, we have allowed a distance of 20 Å

between wires in different cells.

Fig 1 shows the energies of monatomic wires with a planar zigzag and linear structure, as a function of the wire length, for different elements. To facilitate the comparison, lengths are given in units of the linear-chain equilibrium length, and energies are relative to its minimum. These are: 2.56 Å and 2.09 eV/atom for Au, 2.25 Å and 2.16 eV/atom for Cu, 4.10 Å and 0.42 eV/atom for K, and 4.00 Å and 0.98 eV/atom for Ca. The evolution of the bond length in the zigzag structure can be seen in Fig. 2, where we have again used reduced units. Hereafter we will use $`x`$ to refer to the reduced length, and $`d`$ for the reduced bond distance of the wires (see the inset of Fig. 2).

First, we will briefly comment on the data for Au, displayed in Fig 1a. They present two peculiar behaviors: i) The planar zigzag is more stable than the linear geometry, even before the latter reaches its optimum bond length, and ii) the zigzag structure presents two energy minima as a function of the wire length. While a zigzag structure is naturally expected under compression, the first point mentioned above implies that the zigzag distortion, which increases the bond length, is favorable even when the linear wire is stretched. This is quite surprising in a metallic system like gold where, in principle, we do not expect any important bond directionality effects. Fig. 1b shows that this result also holds for Cu wires. As clearly shown in the insets of Fig. 1a and b, in both cases the linear configuration is, at its equilibrium length ($`x`$=1), already $``$0.1 eV/atom less stable than the zigzag structure. In spite of this energy gain, Fig. 2 shows that for both materials, at $`x`$=1, the nearest-neighbors distance increases in $``$5 % when going from the linear to the zigzag conformation. The behavior of the energy shown in panels (c) and (d) of Fig. 1, corresponding to Ca and K, is quite different. In those curves we do not find any sign (at least within the precision of the calculation) of the zigzag distortion until the linear chain starts to be compressed. This is what we could call a ‘trivial’ zigzag distortion: the energy minimum of the linear wire approximately becomes an inflection point for the energy of the zigzag configuration. This is in fact the kind of behavior predicted with, for example, a simple pairwise potential.

All the curves in Fig. 1 present a clear minimum near $`x0.5`$. This minimum corresponds to a zigzag formed by equilateral triangles, where each atom has four nearest-neighbors at the same distance. In Fig. 2 the equilateral configuration corresponds to the points where each curve intersects the dashed line, i.e. where $`d=2x`$.

The second peculiarity of the Au wires mentioned above, i.e. the presence of two minima, is not observed for the other species studied. In ref.  we explained this behavior with a simple jellium model, enphasizing the relationship between the energy vs length behavior, and the crossing of the Fermi level by different bands. For linear wires of Cu and Au with $`x<1.1`$, a doubly degenerated $`\pi `$ band (with mainly $`d`$ character) rises through the Fermi level (Fig. 3a). In the case of Au linear wires, this was already reported in ref. , and later confirmed in ref. . The high density of states at the Fermi level, associated to this band, destabilizes the linear chain, and is the origin of the ‘anomalous’ zigzag behavior of Au and Cu (as opposed to the ‘normal’ behavior of Ca and K), which begins when the crossing occurs.

For $`x>0.7`$, Au, Cu, and K zigzag-wires have a single band crossing of the Fermi level (Fig. 3b), at the Brillouin zone edge. This situation might, in principle, favor a Peirls distortion. However, for the wire lenghts studied here, we have not found any measurable dimerization. As the wire is contracted below $`x0.7`$ (the exact positions are indicated by arrows in Fig. 1) a new conduction band crosses the Fermi level. This is shown in Fig. 3c for the Cu wire. For the three elements (Au, Cu and, K), this new band crossing precedes a change in the behavior of the energy, and may be interpreted as the begining of a chemical bond between second neighbors.

The case of the Ca zigzag-wire is different. The wire is semiconducting for most of the length interval studied here, and the bands cross the Fermi level at a much lower $`x0.45`$ (below the equilibrim length) making it metallic. At $`x0.75`$ its band structure suffers another qualitative change: the gap, that was direct at $`\mathrm{\Gamma }`$ for larger lengths, begins to move towards the edge of the two-atom Brillouin zone, still being approximately direct.

While the energy curves of K and Ca zigzag-wires were relatively similar, the behavior of the bond distance, displayed in Fig. 2, is completely different for metallic and semiconducting wires. In the case of the metallic wires (Au, Cu and, K) the bond distance increases as the zigzag wire is contracted. This can be understood as a consequence of the metallic bonding: as the second-neighbour distance decreases, the first-neighbor distance must increase to maintain the optimum value of the effective atomic coordination. In fact, the calculated behavior is qualitatively reproduced by semiempirical metallic potentials like the glue model , or the embedded atom method , or even a jellium model where the volume per atom remains fixed. However, for the semiconducting wire (Ca), the nearest-neighbor distance behaves just in the opposite way, decreasing as the wire is compressed. It is worth noting here that the predictions of a simple pair potential, like the Lennard-Jones , lie exactly in between these two regimes, i.e. the nearest-neighbors distance would remain almost constant once the wire assumes the zigzag configuration. In Fig. 2, it is also clearly visible the abrupt change in the behavior of the bond length of K wires induced by the band crossing at $`x0.7`$.

In conclusion, we have applied first-principles density-functional calculations to study the energetics, and the elastic and electronic properties of monatomic wires of different elements, in a linear and a planar zigzag geometry. Surprisingly, for Cu and Au wires, a zigzag distortion is favorable even when the linear wire is stretched. This phenomenon is not observed for K and Ca wires, and might be related with the presence of $`d`$ bands (with a high density of states associated) at the Fermi level for the linear conformation of Au and Cu wires. For all the studied materials, the equilibrium structure is an equilateral zigzag (bond angle of 60o). Only in the case of Au, the zigzag geometry can also be stabilized for an intermidiate bond angle of 131o. The behavior of the nearest-neighbor bond length is qualitatively different in the metallic (Au, Cu and, K) and semiconducting wires (Ca).

D.S.P. is grateful to Richard M. Martin for advice and support. This work was supported by Grants No. DOE 8371494, DEFG 02/96/ER 45439. J.M.S. acknowleges a sabbatical fellowship from Spain’s SEUID.


text:
# On the Significance of the Quantity 𝐴² MPI-PhT/00-40

## 1 Introduction

Physical quantities should be gauge invariant. At first glance this might seem to imply that only expressions involving the fields ($`E`$ and $`B`$ in electromagnetism) and not the potentials ($`A`$) should appear in physically meaningful quantities, and in fact this is usually true. However, this logic can be misleading. A well known case in point is the loop integral $`\mathrm{𝐀𝐝𝐱}`$. Although only $`A`$, and not the fields, appears explicitly in this construction, it is so devised that it leads to a gauge invariant and indeed an interesting object.

We would like to point out the interest of another quantity constructed from $`A`$ itself: the volume integral of $`A^2(x)`$. One may come upon this thought when considering the role of condensates in quantum field theory. Vacuum condensates have been a useful way to understand and characterize the dynamics of QCD and other field theories. The most famous example is perhaps the quark condensate:

$$0|\overline{q}q|00,$$
(1)
where $`q`$ stands for light $`u`$\- or $`d`$-quarks. In the realistic case of negligibly small quark mass, a nonvanishing value of the quark condensate signals spontaneous breaking of chiral symmetry.

In the framework of the QCD sum rules  one also used the concept of the gluon condensate

$$0|\alpha _s(G_{\mu \nu }^a)^2|00.$$
(2)
Here the non-vanishing value of the condensate signifies not the breaking of a symmetry but rather the presence of nonperturbative fields in the vacuum.

This gluon condensate would appear to be the simplest quantity characterizing nonperturbative vacuum fields. It has dimension $`\mathrm{d}=4`$, leading one to assume that the leading non-perturbative corrections in the QCD sum rules at large external momentum $`Q`$ are of order $`0|\alpha _s(G_{\mu \nu }^a)^2|0/Q^4`$.

Now there is of course an even simpler candidate for a condensate, namely just the square of the vector potential: $`A^2`$. This is of dimension $`\mathrm{d}=2`$. However, such expressions seem not to be allowed since they appear gauge noninvariant . That is, one tends to think that physically meaningful quantities must involve only the fields and not the potentials, and that an expression like $`A^2`$, involving only potentials, could not be meaningful. However this is not necessarily true, as we would now like to illustrate on the simple example of magnetostatics.

## 2 Magnetostatics

Consider a situation with some magnetic field $`𝐁`$ present in space. There is a considerable amount of freedom in the choice of $`𝐀`$. However since there is a nonzero magnetic field $`𝐁=\times 𝐀`$, we know some nonzero $`𝐀`$ must be present; $`𝐀`$ cannot be zero everywhere. Now consider the volume integral of $`A^2(x)`$. It is a positive quantity and cannot be zero. It must then have some minimum value. Therefore of all the possible $`𝐀`$ configurations which yield the given $`𝐁`$ the one (or the ones) with the smallest integral of $`𝐀^2`$ has in a sense an invariant significance. We would then like to examine the possible significance of the resulting quantity, the volume integral of $`A^2(x)`$ at its minimum value. We will call this $`A_{min}^2`$ .

The connection between the “minimum $`A^2`$” requirement and a more familiar gauge condition may be seen as follows. Suppose for a given field congifuration that $`𝐀^2d^3x`$ is at its minimum value; then under a gauge transformation it is stationary. Considering $`𝐀𝐀+\varphi `$ for infinitesimal $`\varphi `$ we have $`𝐀\varphi d^3x=0`$ and integrating by parts

$$\varphi 𝐀d^3x+surfaceterms=0.$$
(3)
Since $`\varphi `$ is arbitrary we conclude that, up to the surface terms and the question of local minima in $`A^2`$, the “minimum $`A^2`$” condition is equivalent to the familiar gauge condition

$$𝐀=0.$$
(4)
Not surprisingly the “minimum $`A^2`$” requirement is connected with that gauge condition which is invariant, i.e. makes no reference to any particular direction.

Furthermore it appears that $`A_{min}^2`$ is sensitive to, or measures in some way, topological features of the system under consideration. This is suggested by the comparison of two situations, both with no magnetic field. Let one be simple empty space with $`𝐁=\mathrm{𝟎}`$, while the other has a non-trivial topology with the presence of a tube or string containing magnetic flux, like a “cosmic string” or a vortex in superconductivity. In the first case we have simply no $`A`$ and so $`A_{min}^2`$ $`=0`$. In the second case, due to the flux $`\mathrm{\Phi }`$ in the tube or string

$$\mathrm{𝐀𝐝𝐱}=𝐇𝑑𝐬\mathrm{\Phi },$$
(5)
and $`A`$ cannot be zero in the surrounding space even though the magnetic field is absent. This example suggests that $`A_{min}^2`$ can signal the presence of non-trivial topological structures.

The logical situation concerning $`A_{min}^2`$ resembles somewhat that of the question of the energy of a particle in relativity. The energy of a particle is of course a frame dependent quantity. However the minimum energy, which is the energy in the rest frame, has an invariant meaning, namely the mass. In going to the rest frame of the particle we do make a certain choice of frame, but nevertheless the mass is an undeniably meaningful quantity .

Of course the mass also has an explicitly invariant expression, $`m^2=E^2P^2`$, and the loop integral $`𝐀𝑑x`$ can, via Stokes theorem, be expressed in term of the fields. Analogously, is there an a expression for $`A_{min}^2`$ directly in terms of the fields?

Indeed there is the vector relation

$$𝐀^2(x)d^3x=\frac{1}{4\pi }d^3xd^3x^{}\frac{[\times 𝐀(x)][\times 𝐀(x^{})]}{|𝐱𝐱^{}|}+\frac{1}{4\pi }d^3xd^3x^{}\frac{[𝐀(x)][𝐀(x^{})]}{|𝐱𝐱^{}|}$$

$$+surfaceterms.$$
(6)
Each of the two terms is positive, hence (up to the surface term question) we can minimize the integral of $`𝐀^2`$ by choosing $`𝐀=0`$. With this choice the integral of $`𝐀^2`$ is minimal in accord with our above remarks, and is expressed only in terms of the magnetic field $`\times 𝐀`$:

$$A_{min}^2=\frac{1}{4\pi }d^3xd^3x^{}\frac{𝐁(𝐱)𝐁(𝐱^{})}{|𝐱𝐱^{}|}+surfaceterms.$$
(7)
Thus we can trade, so to speak, apparent locality for explicit gauge invariance.

It will be seen that the arguments of this section carry over to four (or more ) dimensions directly, as long as the metric is euclidean. For example in four dimensions Eq becomes

$$A^2(x)d^4x=\frac{1}{2\pi ^2}d^4xd^4x^{}\frac{[F_{\mu \nu }(x)][F_{\mu \nu }(x^{})]}{(xx^{})^2}+\frac{1}{2\pi ^2}d^4xd^4x^{}\frac{[_\mu A_\mu (x)][_\nu A_\nu (x^{})]}{(xx^{})^2}$$

$$+surfaceterms$$
(8)
## 3  Quantum Field Theory

Returning now to quantum field theory and vacuum condensates, we would like to examine the suggestion that $`A_{min}^2`$ , now the expectation value of an operator, is sensitive to, or measures the presence of topological structures in some way.

A simple model we can investigate in this regard is “photodynamics”, i.e. the theory with the Lagrangian density

$$L=\frac{1}{4e^2}(F_{\mu \nu })^2.$$
(9)
This model can be studied in two realizations, compact and noncompact. While the noncompact realization is just the theory of free photons, it is known that the compact realization has nontrivial properties, including a phase transition near $`e^21`$ with a condensation of magnetic monopoles  (for review see, e.g. ). Since the monopoles are the sources of non-zero magnetic flux, we would expect $`A_{min}^2`$ to be sensitive to the phase transition.

We can test these ideas in a numerical simulation by considering the difference of $`A_{min}^2`$ calculated in the two realizations. We take the noncompact theory, given by the action (we work in four euclidean dimensions)

$$S_{non}(F)=\frac{1}{4e^2}d^4x(F_{\mu \nu })^2,$$
(10)
and the compact theory where

$$S_{com}(F)=\frac{1}{2e^2}d^4x[1cos(F_{\mu \nu })^2],$$
(11)
and we would like to examine the difference

$$\zeta (e^2)=𝒟AA^2e^{S_{com}}𝒟AA^2e^{S_{non}}.$$
(12)
where $`A^2`$ is at its minimum for each gauge equivalent configuration.

We do this in a lattice formulation, using a $`12^4`$ lattice with periodic boundary conditions. $`A^2`$ is then measured in units of the lattice spacing. $`𝒟A`$ is normalized so that $`𝒟Ae^S=1`$, with $`A(x)`$ running essentially from $`\mathrm{}`$ to $`+\mathrm{}`$ in the noncompact case and from $`\pi `$ to $`+\pi `$ in the compact case. The “minimum $`A^2`$” condition is enforced by an iterative procedure: given a certain $`A`$ configuration on the lattice links, a gauge function $`\alpha (x)`$ (giving a new potential, $`A\alpha `$) is repeatedly adjusted so as to reduce the volume integral of $`A^2`$ . Each pass works outward from an arbitrary lattice point and the procedure stops when the reduction is less than a certain amount.

Fig. 1 shows the results of the numerical simulation. The sharp jump in $`\zeta (e^2)`$ at the phase transition supports the idea that $`A_{min}^2`$ is a measure of the presence of the monoples and their associated strings, present for $`e^21`$. The fact that $`\zeta `$ jumps to zero is related to the particularly simple aspect of this model, that the small $`e^2`$ sectors of the compact and noncompact theory have the same behavior.

With these numerical calculations we have studied the ground state. When one inserts an external monopole, it is also possible to show the differing response of $`A_{min}^2`$ in the two theories by analytic arguments .

Many open and interesting questions remain, particularly concerning the nonabelian case and the role of a d=2 condensate in QCD. We hope to deal with some of them in future work .

## Acknowledgments

We are thankful to V.A. Novikov and M.I. Polikarpov for useful discussions.