text:
# On the Rapid Estimation of Permeability for Porous Media Using Brownian Motion Paths

## I Introduction

Much theoretical effort has been expended in attempts to either estimate or bound the fluid permeability of a porous medium, given its average statistical properties. Techniques used for this include the volume-to-surface-area ratio , the $`\mathrm{\Lambda }`$ parameter , percolation ideas , and concepts from the theory of Brownian motion such as mean survival distance or mean survival time (inverse reaction rate) .

One basic class of theoretical models allows one to study either packed beds or consolidated porous media such as sandstone. These models consist of ensembles of equal-sized impermeable spherical inclusions immersed in a completely permeable medium. When used to model a packed bed, these spheres are non-overlapping; when used to model a consolidated porous medium, they are randomly located and freely overlapping. A second basic class of models for packed beds and/or porous media consist of periodic arrays of impermeable spherical inclusions immersed in a freely permeable medium. Depending on the ratio between the diameter of the spherical inclusions and the separation between neighboring inclusions, and the periodic lattice chosen, these models may consist either of overlapping or non-overlapping spheres.

The use of Monte Carlo methods to solve elliptic or parabolic partial differential equations (PDEs) is grounded in the mathematical subject called probabilistic potential theory. It is well known that the probability density function of a diffusing particle obeys the diffusion equation. The time-integrated, or steady-state probability density function obeys the Laplace equation. These facts have been greatly generalized. One can solve a large class of elliptic PDE’s by simulating the corresponding class of (in general) biased diffusion processes . The diffusion Monte Carlo method is defined by the mapping of boundary value problems involving elliptic PDE’s into the corresponding diffusion problems, and solving them by using Monte Carlo methods. This method is powerful because it can be implemented computationally in a massively parallel manner.

Previous work involving one of us (J.G.) develops a class of diffusion Monte Carlo algorithms, the Green’s function based first-passage (GFFP) methods that are very efficient for studying the properties of random or highly irregular two-phase media. This class of problems involves:

$``$ Interfacial surfaces that may pass very close to each other, thus creating channels that are narrow, but nevertheless important for understanding the bulk properties of the material.

$``$ Spatial intermittency: large regions with no interfacial surface; smaller regions with a high density of interfacial surface. Intermittency may be due e.g. to polydispersity in pore diameter or void size, or to a dispersion, in solution, of very large irregular inclusions, such as macromolecules.

$``$ Interfacial surfaces that may include sharp corners and edges. The difficulties these features cause in solving the Laplace equation are clear (In electrostatic terminology, they cause singularities in charge density.).

In particular, the class of GFFP algorithms,:

$``$ Uses a purely continuum description of both sample and diffusing particle trajectory, with no discretization either in space or in time and hence no discretization error. This is essential for calculations of permeability in low-porosity materials.

$``$ Uses a library of Laplace Green’s functions that allow a diffusing particle both to make large ‘jumps’ across regions containing no interfacial surface, and to be absorbed from a distance, at an absorbing boundary, in a single jump.

$``$ Uses ‘charge-based’ rather than ‘voltage-based’ methods of calculation, i.e., retain only the location of absorption events, or equivalently surface charges. This allows one to perform trajectory simulation using importance sampling (because each trajectory yields one surface charge); thus ensuring that even singularities associated with corners and edges will be efficiently sampled.

This method is an adaptation and refinement of the "walk on spheres" (WOS) method of diffusion Monte Carlo . This method has been used since at least 1956 to treat diffusion problems . It was first applied in random-media problems by Zheng and Chiew ; see also . Other recent studies have also generalized WOS to diffusion in two-phase media modeled as arrays of small homogeneous squares  and cubes ;these studies use the simplest nontrivial Green’s functions to determine first-passage probability. The Einstein relation may be used to obtain effective electrical conductivities  from the bulk diffusion coefficient of these models.

The present method is more efficient than WOS for diffusing particles very close to absorbing boundaries; as we show here, this makes it ideal for the study of low-porosity materials. It has been applied to solve the scalar Laplace equation, which governs both electrostatic and diffusion problems.

Electrostatic applications involve the calculation of electrostatic capacitance for very general bodies . This is an important set of problems because probabilistic potential theory implies fundamental relations between the capacitance and the lowest eigenvalue of the Laplacian in the Dirichlet problem . This in turn leads to efficient bounds or approximations for a wide variety of physical quantities in many-body systems, for example, the quantum mechanical scattering length in nuclear matter, or the classical rate of heat release in composite media.

Diffusive applications of this class of methods involve solutions of the Smoluchowski equation, either in the vicinity of a single irregular object, (e.g., a ligand diffusing near a reactive site on a macromolecule) , or in a system containing many absorbing objects, (e.g., systems with distributed trapping) .

Because so many physical phenomena are governed by vector Laplace equations (e.g. low-velocity fluid dynamics) , and tensor Laplace equations (e.g. solid mechanics) , it is natural to seek methods for extending this previous work to deal with these two classes of problems . An approach to the first class is provided by the method of angle-averaging . This method provides an approximate relation between the translational hydrodynamic friction $`f`$ of an arbitrary body, and the electrostatic capacitance $`C`$ of a perfect conductor of identical size and shape. For arbitrary nonskew objects with a no-slip boundary condition, using angular averaging of the Oseen tensor, Hubbard and Douglas  generalized the exact Stokes-Einstein result for the translational hydrodynamic friction of a sphere,

$$f_{sphere}=6\pi \eta R,$$
(1)
where $`\eta `$ is the solvent viscosity and $`R`$ is the radius of the sphere, to the translational hydrodynamic friction of an arbitrary impermeable body to give:

$$f=6\pi \eta C$$
(2)
where $`C`$ is the corresponding capacitance. This formula is exact in the "free tumbling" approximation, i.e., in the case that there is no coupling between the translational and rotational motion, and that there is sufficient noise, or disorder, in the flow to cause the object to occupy all possible orientations equally. But extensive tests show the method is accurate well beyond this set of cases.

A central insight of the work described here is this: the angle-averaging method may be applied to approximate the translational friction coefficient of an "object", even if the object is not a connected set, e.g., if it is a collection of inclusions that constitute the impermeable or matrix phase of a sample of porous media. A group of well-studied models for porous media fit this description. The inclusions may either be randomly placed, or set in a periodic lattice. The former models have properties similar to many commercially important porous media; the latter model a set of media formed by grain consolidation. The inclusions may be freely overlapping or impenetrable, i.e., unable to overlap. Finally the (spherical) inclusions may be uniform (monodisperse) in radius, or they may be polydisperse.

In this paper, we present two classes of accurate, computationally efficient methods of calculating permeabilities for these models and models like them. These models combine rapid efficient methods of simulating Brownian motion  with a pair of methods for deriving the permeability from the statistics of Brownian particles diffusing near a sample of the porous medium. The trajectories of these particles initiate outside the sample and terminate on contacting the porous matrix. The spherical collection of small inclusion spheres will occasionally be referred to as the porous sample in subsequent discussions. The first method, the "penetration depth" method, associates the fluid-dynamic penetration depth with a specific property of the Brownian paths, then uses the standard relation between penetration depth and permeability to calculate the latter. The second method, the "unit capacitance" method, involves using Brownian paths to calculate an effective capacitance for the sample, then relates the capacitance, via angle-averaging theorems, to the translational hydrodynamic friction of the sample. Finally, a result of Felderhof  is used to relate the latter quantity to the permeability of the sample.

This paper is organized as follows: Section II sketches the method used here to obtain the electrostatic capacitance using Brownian paths. Section III develops the penetration depth method of calculating permeability, in which we obtain the penetration depth from Brownian paths, and obtain the permeability from the penetration depth. Section IV describes the unit capacitance method of calculating permeabilities, and applies it to the random models of packed beds and consolidated porous media. Section V applies both methods to determining the permeabilities of the random models of packed beds. Section VI applies both methods to determining the permeabilities of model packed beds composed of periodic arrays of impermeable spheres. Section VII applies both methods to determining the permeabilities of polydisperse randomly packed bed. Section VIII contains our conclusions and suggestions for further study.

## II Calculation of the Capacitance of a Sample of General Shape Using Brownian Paths

In this section, we sketch the algorithm for calculating the electrostatic capacitance of an irregularly shaped conducting body by using Brownian trajectories that initiate on a spherical launch surface surrounding the conducting body, and terminate on contact with that body. We focus this derivation in several ways in order to apply it to samples of porous media. If a packed bed model is either randomly packed or densely packed, any smooth sample boundary will intersect some of the spherical inclusions that constitute the sample; we provide two methods for dealing with this problem in section IV, the "effective radius method" and the "sharp boundary method". Difficulties will occur if the sample is chosen either too large or too small; we give an analysis that yields the acceptable qualitive range of sample sizes in section IV. For periodic models of porous media, a non-round, e.g. cubic, sample (and thus a non-round launch surface) might be more appropriate. The derivation provided here is readily generalized to allow a cubic launch surface; in forthcoming work we explore this possibility numerically.

The electrostatic potential $`\varphi (r)`$, at position $`r`$, in the presence of a conducting object, and the probability density $`p(r)`$ in the associated Smoluchowski diffusion problem, i.e., the probability associated with finding a diffusing particle at position $`r`$, in the presence of an absorbing object of size and shape identical to the conductor, are related by: $`\varphi (r)=1p(r)`$. This equation is derived by noting that its LHS and RHS obey the Laplace equation with identical boundary conditions.

Thus the capacitance $`C`$ of a conducting object is given by:

$$C=(4\pi )^1_\mathrm{\Omega }𝑑\sigma \varphi (r)=(4\pi )^1_\mathrm{\Omega }𝑑\sigma p(r).$$
(3)
where $`\mathrm{\Omega }`$ is the surface of the object.

The diffusion-controlled reaction rate $`\kappa `$ is defined as

$$\kappa =D𝑑\sigma p(r)$$
(4)
where the integral is evaluated on any closed surface containing the inner boundary $`\mathrm{\Omega }`$ (because of Green’s theorem any such surface may be used.).

The relation between the diffusion-controlled reaction rate, in the Smoluchowski problem, of particles absorbed upon contact with an object, and the electrical capacitance of a conductor with identical size and shape to that object is thus given by:

$$\kappa =4\pi DC$$
(5)
We define $`\beta _{\mathrm{}}(\stackrel{}{r};\mathrm{\Omega })`$ to be the probability that the diffusing particle started at $`\stackrel{}{r}`$ will be absorbed on $`\mathrm{\Omega }`$ and $`\gamma _{\mathrm{}}(\stackrel{}{r};\mathrm{\Omega })`$ the probability that the diffusing particle started at $`\stackrel{}{r}`$ will go to infinity, i.e., that it will not be absorbed in finite time. Thus:

$$\gamma _{\mathrm{}}(\stackrel{}{r};\mathrm{\Omega })=1\beta _{\mathrm{}}(\stackrel{}{r};\mathrm{\Omega }).$$
(6)
$`\gamma _{\mathrm{}}(\stackrel{}{r};\mathrm{\Omega })`$ can be shown to satisfy the Laplace equation . If the obvious boundary conditions, $`\gamma _{\mathrm{}}(\mathrm{};\mathrm{\Omega })=1`$ and $`\gamma _{\mathrm{}}(\stackrel{}{r};\mathrm{\Omega })=0`$ for $`\stackrel{}{r}`$ on $`\mathrm{\Omega }`$ are considered, it is clear that

$$p(r)=\gamma _{\mathrm{}}(\stackrel{}{r};\mathrm{\Omega }).$$
(7)
Using Eqs. 6 and 7, Eq. 4 becomes:

$$\kappa =D_r𝑑\sigma \beta _{\mathrm{}}(\stackrel{}{r};\mathrm{\Omega })$$
(8)
$$=Dr^2\mathrm{sin}\theta d\theta d\varphi \frac{}{r}\beta _{\mathrm{}}(r;\mathrm{\Omega }),$$
(9)
$$=r^2D\frac{}{r}\mathrm{sin}\theta d\theta d\varphi \beta _{\mathrm{}}(r;\mathrm{\Omega }),$$
(10)
Defining

$$\beta (r;\mathrm{\Omega })=\frac{1}{4\pi }\mathrm{sin}\theta d\theta d\varphi \beta _{\mathrm{}}(r;\mathrm{\Omega }),$$
(11)
we obtain from Eq. 10

$$4\pi r^2D\frac{}{r}\beta (r;\mathrm{\Omega })=\kappa .$$
(12)
For any spherical surface $`r=b`$ which contains the boundary $`\mathrm{\Omega }`$, integrating Eq. 12 from $`b`$ to $`r`$ with respect to $`r`$ and using the boundary condition $`\beta (\mathrm{};\mathrm{\Omega })=0`$, gives an expression for the reaction rate $`\kappa `$:

$$\kappa =4\pi Db\beta (b;\mathrm{\Omega })$$
(13)
Comparing Eqs. 5 and 13 yields an important formula for the electrostatic capacitance

$$C=b\beta (b;\mathrm{\Omega })$$
(14)
$$b\beta $$
(15)
where $`\beta `$ is the fraction of diffusing particles started at a random, i.e., angle-averaged, position on the launch sphere that are absorbed on the target. The unit capacitance method, as described here, uses this relation to determine the capacitance of a packed bed by performing Monte Carlo simulations for the quantity $`\beta `$. Note that no assumption is made in the above that the "object" being studied is a connected set; it may be taken to be the set of inclusions that constitute the matrix phase of a sample of porous media.

We calculate the quantity $`\beta `$ by performing simulations as follows: each random diffusion path is constructed as a series of first-passage propagation jumps, each from the present position of the particle to a new position on a first-passage surface drawn around the present position. The new position is sampled from a first-passage position distribution function $`f(x)`$ which gives the probability density associated with finding a diffusing particle leaving the present position and first contacting the first-passage surface at point $`x`$.

This simplest first-passage surface is the first-passage sphere, i.e., a large sphere centered around the present position of the diffusing particle, that does not intersect any of the inclusions. The first-passage position density on this surface is uniform, i.e., its distribution function is trivial. Using only this first-passage surface yields a trivial case of the GFFP method, namely the WOS method. In this method, whenever a diffusing particle gets close enough (within a fixed distance $`\delta _H`$) to one of the inclusions, it is taken to be absorbed.

We show here that for the problems under study it is more efficient to use more complex first-passage surfaces. We let the first-passage sphere intersect the nearest inclusion, and grow as large as possible, provided that it:

$``$ not intersect the next-nearest inclusion

$``$ not intersect any corners or edges of the nearest inclusion

The resulting first-passage surface is the portion of the first-passage sphere not contained in any inclusion. Its surface consists of part of the first-passage sphere, and part of the surface of the inclusion. The probability distribution function $`f(x,\stackrel{}{P})`$, for a diffusing particle starting at the center of the first-passage sphere and making first-passage at the point $`x`$ is in general quite nontrivial. (Here $`\stackrel{}{P}`$ is the set of geometric parameters that characterize the first-passage surface. Also, we assume for simplicity, that symmetry allows specification of the first-passage position $`x`$ by a single position parameter, the polar angle $`\theta `$.) Thus, $`(\stackrel{}{P},\theta )`$ is the full set of parameters that the first-passage distribution function depends on.

For a wide class of inclusions, the resulting first-passage surfaces are ’locally simple’ closed geometrical surfaces. These surfaces are those for which we can tabulate, as a function of several parameters, the corresponding first-passage distribution function. To obtain this function, we:

$``$ Calculate the gradient, normal to the absorbing surface, of the appropriate Laplace Green’s function for a grid of parameter and field values.

$``$ Calculate the indefinite integral of this gradient as a function of $`\mathrm{cos}\theta `$. Normalize the integral to be unity:

$$I(\stackrel{}{P},\mathrm{cos}\theta )=1$$
(16)
for $`\mathrm{cos}\theta =1`$.

$``$ Invert this relation for each set of parameter values $`\stackrel{}{P}`$ to obtain $`\mathrm{cos}\theta `$ as a function of $`\stackrel{}{P}`$ and $`I`$. The first-passage position, i.e., $`\mathrm{cos}\theta `$, can then be importance-sampled by choosing a random number $`\alpha [0,1]`$, setting $`I=\alpha `$, and interpolating this relation.

If the chosen first-passage position is on the portion of inclusion surface included in the first-passage surface, the particle is absorbed. Otherwise, it will reach a point on the rest of the first-passage surface. From there it makes its next jump. The diffusing particle continues until it is absorbed in the absorbing target or goes to infinity.

As shown in Fig. 1, the GFFP method is more efficient than the WOS method whenever a small distance $`\delta _H`$ must be used; this is true because the average number of moves required for absorption in the WOS method increases rapidly as $`\delta _H`$ decreases . Thus, it is more efficient for studying porous media either at low porosities or high polydispersity.

The set of tabulated first-passage distribution functions is built up by a process of ‘bootstrap diffusion Monte Carlo’, in which simple distribution functions are used to accelerate the simulations used to calculate more complex distribution functions, and so on. In practice, the set of distribution functions is limited mainly by our inability to store in cache memory a tabulated distribution function depending on more than three parameters. In the present study the inclusions are asymmetric lenses (in the sharp-boundary method) and spheres. These are readily treated by the GFFP method .

## III The Penetration Depth Method of Calculating Permeability from the Absorption Point Statistics of Brownian Trajectories

In this section, we describe the ‘penetration depth’ (PD) method for calculating permeability of an absorbing sample. First we relate the fluid-dynamic concept of "penetration depth" to a property of the paths of Brownian particles that diffuse from outside the sample of porous media and terminate on contact with it. Then we use the standard relation between the permeability of a sample and its penetration depth to calculate the latter.

A simple model calculation allows one to relate the permeability of a porous medium to the effective penetration depth of fluid into the medium. Consider a homogeneous porous medium in a half-space with constant permeability $`k`$. Fluid permeates this medium and has a constant velocity $`V_0`$ within it. If $`z`$ is the distance between the macroscopic fluid-medium interface and a point of the porous medium, the Debijf-Brinkman equation

$$P=\eta ^2𝐕\frac{\eta }{k}𝐕.$$
(17)
where $`𝐕`$ and $`P`$ are macroscopic velocity and pressure respectively,

reduces to the following,

$$\eta \frac{d^2V}{dz^2}=\frac{\eta }{k}V,$$
(18)
with the solution

$$V(z)=V_0\mathrm{exp}(z/\sqrt{k})$$
(19)
This result shows that $`\sqrt{k}`$ measures the distance that the flow effectively penetrates into the porous medium. We will call $`\sqrt{k}`$ the ‘penetration depth’.

We can use this result to measure the permeability of a sample by identifying the penetration depth with the difference $`l`$, between the average radial position at which the diffusing particles are absorbed and the actual sample radius, thus yielding the approximate relation:

$$k=l^2$$
(20)
The penetration depth defined here for diffusing particles is different from the mean linear survival distance: the latter is the average distance from the random starting point in the void phase to the absorption point; the former is the average of the radial component of this distance. It is possible that the mean linear survival distance can give an upper bound on permeability analogous to the Torquato-Kim upper bound  based on the mean survival time.

## IV The Unit Capacitance Method of Calculating Permeability of Random Systems

In this section, we describe an algorithm, the unit capacitance (UC) method, for calculating the permeability of a sample of random media. We calculate the capacitance of a sample (as described in Section II), use Eq. 2, from angle-averaging to give the translational hydrodynamic friction in terms of the capacitance, and finally use a relation, first published by Felderhof, between the translational hydrodynamic friction and the permeability to determine the latter.

Assuming that a given spherical sample is much larger than either the average distance between spherical inclusions or the correlation length associated with the statistics of the packed bed, it can be modeled as a homogeneous porous sphere with the appropriate porosity. One can solve the linear Stokes equation  for the translational frictional coefficient, $`f`$, of such a sphere by making the assumption that the excess pressure is linear in the fluid velocity, and by using the symmetry of the problem. The result is:

$$f=6\pi \eta RG_0(\sigma )\{1+\frac{3}{2\sigma ^2}G_0(\sigma )\}^1$$
(21)
where $`\eta `$ is the fluid viscosity and the function $`G_0(\sigma )`$ is given by:

$$G_0(\sigma )=1\frac{1}{\sigma }\mathrm{tanh}\sigma .$$
(22)
Here $`\sigma `$ is the dimensionless quantity defined by

$$\sigma =\frac{R}{\sqrt{k}},$$
(23)
where $`R`$ is the porous sphere radius and $`k`$ is the permeability.

Eliminating the translational frictional coefficient between Eqs. 2 and 21, one finds a relation between the capacitance and the permeability

$$\frac{C}{R}=G_0(\sigma )\{1+\frac{3}{2\sigma ^2}G_0(\sigma )\}^1.$$
(24)
Obtaining $`C/R`$ (unit capacitance) from simulation allows us to use Eq. 24 to calculate $`\sigma `$, and thus obtain the desired permeability estimate from Eq. 23.

An important finding of the present work is that, even though neither the graph of C/R vs. R, nor that of $`\sigma `$ vs. C/R contains a flat region, the composition of the two, giving our permeability estimate as a function of R, does contain a flat "plateau" region. For example, Fig. 2 shows that for a broad range of porosities there is a plateau range of sample radii over which the permeability estimate does not vary. The lower and upper limits of this region vary systematically.

In monodisperse homogeneous porous media models, the permeability is dependent only on the porosity (defined as the ratio of void space volume to total volume) and the inclusion sphere radius. By dimensional analysis, it can be expressed in the following functional form:

$$k=f(\varphi _1)a^2$$
(25)
where $`\varphi _1`$ is porosity and $`a`$ is inclusion sphere radius. The permeability estimates produced by the UC method will, in addition, depend on the parameter $`R/a`$. However, once this parameter is set larger than both the correlation length of the medium (very small in random media) and the average distance between inclusions, this dependency should become quite weak.

In Eq. 24, for a given porosity and inclusion sphere radius (i.e. the permeability is constant), if we increase the sample radius, the unit capacitance goes to unity as $`\sigma `$ goes to infinity. For each porosity, there is a range of sample radii which are much bigger than the average distance between spherical inclusions, but for which the ratio $`C/R`$ will be far enough from unity to permit us to interpolate the corresponding $`\sigma `$ value. We choose a radius that is in this range.

When C/R is close to unity, that is, for very low porosities, the $`\sigma `$ interpolation becomes unstable: a small change in unit capacitance causes a large change in $`\sigma `$. This makes the permeability estimate in the unit capacitance method difficult; thus the penetration depth method will prove to be more reliable at very low porosities.

A fundamental problem of this study is that of defining the sample boundary. To see this, note the following paradox: we study the bulk property (permeability) of a medium using diffusing particles that start outside a finite sample and are absorbed on its surface. These particles will in general be absorbed in a surface boundary layer that grows thinner as the porosity decreases. Thus, we face the problem of freeing our simulation from surface artifacts so that it can determine this bulk property, when the method is based on events that occur near the surface! In this paper, we show that this is indeed quite possible.

We generate samples of the random media models that we study using two different methods: the "effective-radius" method, and the "sharp-boundary" method. We describe each method in turn.

In the "effective-radius" method, if randomly overlapping, i.e., when Poisson statistics are used, we place inclusion sphere centers at random in the sample sphere of radius $`R`$. We continue to place centers until the density of sphere centers approaches $`\rho `$, where $`\rho `$ and the porosity $`\varphi _1`$ being modeled are related by:

$$\varphi _1=\mathrm{exp}[(4\pi /3)R^3\rho ]$$
(26)
If random non-overlapping statistics are used, we place centers sequentially at randomly chosen positions, requiring only that each center placed be far enough from the centers already placed that the corresponding inclusions not overlap. The resulting distribution is not identical to true random non-overlapping statistics, i.e., to what liquid-state theorists call the "hard-sphere" distribution, but studies of these two distributions show that differences in bulk properties emerge only at very low porosities. Even though an inclusion extends in part outside the sample sphere of radius $`R`$, it still contributes both to the density of centers $`\rho `$ and to the simulation.

Because of the lack of the contribution of the inclusion spheres in the region $`r[R,R+2a]`$, the region $`r[Ra,R+a]`$ will have local porosity less than the desired bulk porosity. This is an important concern: especially at lower porosities, the diffusing particles will seldom sample any part of the sample except for the boundary layer. Because of this effect, an effective sample radius, which is expected to be less than $`R`$, is used. For the non-overlapping case, we use for an effective radius the radius of a sphere, which, if it contained the same number of surviving inclusion centers as our sample, would have as its average porosity the same as the bulk porosity of our sample, i.e., the local porosity far away from the sample boundary. (In previous research , we used this choice of effective radius for both the overlapping and non-overlapping case.) For the overlapping case, the same procedure is used; the set of effective sample radii is then rescaled so that the effective radius at the critical porosity ($`\varphi _1^c`$=0.03) makes the calculated unit capacitance equal to 1. (Here by the term ‘critical porosity’, we mean the porosity below which fluid ceases to flow through the sample, i.e., the percolation threshold.)

In the "sharp-boundary" method, instead of using an ad hoc effective radius, the porous media sample is constructed as follows: we place the centers of inclusion spheres into a large sphere of radius ($`R+a`$) according to the chosen statistics, for a given porosity. We then define the actual sample by drawing a sphere of radius $`R`$ and allowing it to freely intersect inclusions already placed. The sample is then defined to be all of the void phase, all inclusions, and all fragments of inclusions, that are contained in this sample sphere of radius $`R`$. With this sharp boundary, the porosity in the actual sample boundary is maintained uniform up to the boundary.

The sampling of periodic grain consolidation models is also done according to the sharp boundary method: we take a spherical sample of radius $`R`$ after making periodic grain models in a cube which is slightly bigger than $`2R\times 2R\times 2R`$. For these models, the result will depend upon the position of the sample center; thus we average the results over ten samples, each with a randomly chosen center.

## V Studies of Randomly Packed Beds

In this section, we apply the two permeability estimation methods to the calculation of permeability of model packed beds, composed of non-overlapping, randomly placed, impenetrable spherical inclusions, and model porous media, composed of randomly placed, freely overlapping, impenetrable spherical inclusions. We compare our results with the available numerical solutions of the Stokes equation, and also with a number of theoretical bounds and estimates from the literature.

Our simulations show that, for the overlapping sphere model studied here, the two permeability estimation methods, the penetration depth method and the unit capacitance method, agree well within random error, estimated, in this case, as the relative standard deviation of the results. For the nonoverlapping sphere model, the two methods deviate as porosity decreases. For both models, we obtain 10 packing bed samples and average the ten permeability estimates.

As already discussed, the graph of permeability shows a very substantial plateau as a function of sample radius, i.e., a region of sample radius over which the permeability estimate shows almost no variation. Also, with the PD method the same property of sampling-size-independent estimation is observed (see Fig. 3).

Our estimation results for the randomly non-overlapping model are compared with some other data sets in Fig. 4. The other data sets used, numbered a) to d), are:

a) The Stokes’ law  for a dilute bed of spheres is

$$k=\frac{2}{9(1\varphi _1)}a^2$$
(27)
where $`a`$ is the inclusion sphere radius.

b) The Happel-Brenner approximation  is,

$$k=(\frac{2}{9\gamma ^3})(\frac{3\frac{9}{2}\gamma +\frac{9}{2}\gamma ^53\gamma ^6}{3+2\gamma ^5})a^2$$
(28)
where $`a`$ is the inclusion radius, and $`\gamma ^3(1\varphi _1)`$.

c) The upper bound of Torquato-Kim  is:

$$k=\frac{2}{9\tau }\frac{\varphi _1}{(1\varphi _1)}a^2$$
(29)
where $`D`$ is the diffusion constant and $`\tau `$ is the average survival time. The data for $`\tau `$ is taken from Ref. .

d) The empirical Kozeny-Carmen relation  in the general case is:

$$k=\frac{1}{180}\frac{\varphi _1^3}{(1\varphi _1)^2}a^2$$
(30)
The Stokes’ law, which is the simplest approximation for the dilute bed, gives estimates that are greater than all other data sets. The Happel-Brenner approximation and the Torquato-Kim upper bound are above our data. Our simulation data with the effective sample radius is above the Torquato-Kim upper bound. Our simulation data with the sharp sample boundary lies between the Happel-Brenner approximation and the empirical Kozeny-Carmen relation. Because our simulation data with the effective radius method violates the Torquato-Kim upper bound, the sharp boundary method is clearly to be preferred for this case.

Our estimation results for the overlapping case $`R=15.0`$ in Fig. 5 are compared with some other data. The estimation is not applied at very low porosities ($`\varphi _1`$ \< 0.1) and very high porosities ($`\varphi _1`$ \> 0.9).

The other data sets used, numbered e) to g), are as follows:

e) The Torquato-Kim upper bound on permeability  for an overlapping sphere bed is,

$$k\frac{D\tau }{F}$$
(31)
where $`\tau `$ is the average survival time and $`F`$ is the formation factor. Data for $`\tau `$ and $`F`$ are obtained from Refs.  and .

f) The Doi upper bound  on permeability for this system is

$$k=\frac{2}{3}a^2e^\gamma _0^1𝑑x(\frac{1}{3\gamma }+x(1x)^2)\mathrm{exp}(\gamma (\frac{3}{2}x\frac{1}{2}x^3)),$$
(32)
where $`\gamma `$ is given in terms of the porosity by $`\varphi _1`$ = exp(-$`\gamma `$).

g) The lattice-Boltzmann data , which agrees well with the Kozeny-Carmen equation,

$$k=\varphi _{1}^{}{}_{}{}^{3}/6s^2.$$
(33)
The lattice-Boltzmann method is a variant of the lattice-gas automata method. The specific surface area $`s`$ is given analytically by $`s=\mathit{4}\pi a^\mathit{2}\rho \varphi _\mathit{1}`$ where $`\rho `$ is the number of sphere centers per unit volume. At low porosity, finite size effects occur .

Simulation data with the sharp boundary sampling method are below the two upper bounds, the Doi upper bound and the Torquato-Kim upper bound. Note that the lattice Boltzmann result at high porosities and the actual data for monosized sphere beds at low porosities agree well with our simulations. However, the PD method overestimates permeabilities at very low porosities because this method does not take into account the phenomenon of percolation. Our data from the effective radius method is quite good for the overlapping case.

At extremely low porosities, estimation instability appears; when the porosity is very close to unity, simulation may give a sample unit capacitance slightly above 1.0 due to the Monte Carlo fluctuations. The unit capacitance data which exceed unit value are discarded. The simulations described here are large enough that this effect is quite small except at some low porosities of big porous samples like $`R=50.0`$.

## VI Studies of Periodic Grain Consolidation Models of Porous Media

In this section, with the sharp boundary sampling method, we apply the two methods of determining permeability that we have developed to the permeability calculation of model packed beds, composed of non-overlapping, periodic, impermeable spherical inclusions, and model porous media, composed of periodic overlapping, impermeable spherical inclusions. The difference between the two types of models consists only in the ratio of the sphere diameter to the separation between the spherical inclusions. We compare our results with the solution of the Stokes equation, using fluid dynamics codes, by Larson and Higdon , for simple cubic (SC), body-centered cubic (BCC), and face-centered cubic (FCC) lattices.

In all cases, 10 packed bed samples are used, and the results are averaged. Figures 6, 7 and 8 compare results from the two methods with the Larson-Higdon  calculations for sample radius $`R=15.0`$. At low porosities, the standard deviation for the unit capacitance estimation is higher because of the estimation instability. With the sampling radius $`R=15.0`$, two estimation results, UC and PD, show some deviations. The PD estimate is better than UC estimate. However, with $`R=50.0`$, the deviation between the two methods becomes smaller. At low porosities, due to the estimation instability of UC, no UC data are shown. It seems that for grain consolidation models UC estimation needs larger sampling due to the cubic lattice structure.

Based on these results, we increase the sampling size to $`R=50.0`$. As already discussed, we cannot use the UC method here. The results from the PD method, Figs. 9, 10, and 11 are in excellent agreement with the Larson-Higdon results.

The same property of sampling-size-independent estimation as the random media is observed in Fig. 12. It should be noted that, for high porosities, the necessary sample radius can be as large as $`R=50.0`$.

## VII Studies of Polydisperse Randomly Porous Media

In this section, we apply the two new methods with the sharp boundary sampling method to the calculation of the permeability of model porous media, composed of polydisperse overlapping, randomly placed, impermeable spherical inclusions. The inclusion sphere radii are chosen at random from the values {1.5, 3.5, 5.5, 7.5}. We compare our results with the available numerical solutions of the Stokes equation .

Figure 13 shows that the two sets of results agree well. At very low porosities, our methods give permeabilities greater than those given by the Stokes equation because we have not yet incorporated the phenomenon of percolation. This will become important when we perform more extensive study of polydisperse models; the percolation threshold seem to be reached at lower porosities in such systems.

## VIII Conclusions and Suggestions for Further Study

Our methods give very good results for all models of porous media tested. The penetration length method is better at low porosities; the unit capacitance method shows high standard deviations at low porosities due to the steepness of the unit capacitance curve. It is very unfortunate that so little high quality simulation data exists even for the simple models studied here. We note that these models have been standards for theoretical study for decades. Our method will predict permeabilities for a large class of homogeneous and isotropic porous media, in the medium and high porosity regimes $`(0.1<\varphi _1<0.9)`$.

The unit capacitance estimation discrepancy with sampling radius $`R=15.0`$ may be due to the importance, for periodic models, of using a launch surface with the geometry of the sample, e.g., a cubic launch surface for SC samples. We have laid the basis for this step here; it will require further study.

An important point of this method is that it is very fast compared with other methods. Our computations were in parallel Fortran, using MPI, on a PC cluster (32 node, 233 Mhz machine). Using 10 nodes, one porous sample per node, one million diffusing particles for each porous sample and porosity, and 10 porosities, it takes about 10 hours per set of calculations. Using e.g. boundary element/finite element methods to solve the Stokes equation in a sample of porous media can require the same amount of time to do this set of calculations for a single value of porosity (We note that the latter methods are, at present, more efficient than these methods developed here for obtaining the detailed flow field in a sample; this may be important in particular applications.).

ACKNOWLEDGMENTS

We give special thanks to Nicos Martys, Eduardo Glandt and Dan Rothman for their useful information and comments. We thank Joe Hubbard and Jack Douglas for useful discussions. We are especially grateful to Nicos Martys for sharing with us the raw simulation data used in reference .

Also, we wish to acknowledge NCSA, especially the Condensed Matter Physics section of the Applications Group, for access to their computational resources and collaborative support. In addition, we thank NCSA’s educational programs for training support to use NCSA computational resources.


format: text

text:
# Dust in active nuclei

## 1 Introduction

The properties of the circumnuclear gas are a key issue to understand the physics of Active Galactic Nuclei (AGNs). In particular, gas obscuration has important consequences on the classification of AGNs (Antonucci 1993), on their infrared emission (Granato et al. 1997, Pier & Krolik 1993) and also on the X-ray background (Setti & Woltjer 1989). AGNs are, optically-wise, divided in two main classes: type 1 AGNs, showing broad permitted emission lines, and type 2 AGNs, which only show narrow emission lines. The Unified Model assumes that AGNs of both classes host the same kind of nuclear engine and ascribes their differences solely to orientation effects with respect to an obscuring gaseous medium, possibly arranged in a torus-like geometry. For those lines of sights intercepting the obscuring torus, both the Broad Line Region (BLR, $`<1`$ pc in size) and the nuclear engine are obscured and only the much more extended Narrow Line Region (NLR) can be observed. This model has gained success from a large number of observational tests (see Antonucci 1993 for a review). In particular, X-ray observations have supported the unified scenario by discovering large columns of absorbing gas in type 2 AGNs (eg. Awaki et al. 1991, Turner et al. 1997, Maiolino et al. 1998, Risaliti et al. 1999). Also, spectroscopic observations in the infrared, where dust absorption is greatly reduced, detected broad permitted lines in several AGNs which are classified as type 2 in the optical.

However, the properties of the absorbing medium result to be more complex with a more quantitative analysis. Maccacaro et al. (1982) first noted, in a few AGNs, that the dust reddening affecting the BLR is significantly lower than expected from the NH measured in the X-rays, assuming a Galactic standard extinction curve and dust-to-gas ratio. Indications for a low value of $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ in the obscuring torus were also found by Granato et al. (1997) by modelling the IR emission of AGNs; they ascribed the low $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratio to the sublimation of dust at the inner face of the torus.

The goal of this paper is to observationally verify the low $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ and low $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ phenomena with a higher confidence (sects. 2 and 3). We also investigate dust spectral signatures which directly probe the properties of dust grains in the circumnuclear region of AGNs. Most of the interpretation of these observational effects is addressed in a companion paper (Maiolino et al. 2000b, paper II). Finally, note that in this paper we will often distinguish the $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ and the $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratios since, as discussed in paper II, in the circumnuclear region of AGNs dust reddening and obscuration are not necessary tied by the Galactic standard relation ($`\mathrm{A}_\mathrm{V}/\mathrm{E}_{\mathrm{B}\mathrm{V}}=3.1`$).

Throughout this paper we will assume a cosmology with $`\mathrm{H}_0=65`$ and $`\mathrm{q}_0=0.5`$.

## 2 Evidence for a low EB-V/NH

### 2.1 Sample selection and reddening determination

We have defined a sample of AGNs whose X-ray spectrum shows evidence for cold absorption (hence a measure of the gaseous column density NH along the line of sight) and whose optical and/or IR spectrum show at least two broad lines that are not completely absorbed by the dust associated to the X-ray absorber. This sample includes intermediate type (1.8–1.9) Seyferts, type 2 Seyferts with broad lines in the near-IR, and a few type 1–1.5 Sy characterized by cold absorption in the X-rays. We do not consider AGNs whose X-ray spectrum shows evidence for a warm absorber, since very likely the latter is not associated to the obscuring torus (Netzer 2000) and also because the gaseous column inferred for a warm absorber is strongly model-dependent (this is discussed further in Sect. 5). We do not consider cold absorbers with partial covering, since also in this case the absorber is not associated to the obscuring torus, but to matter very close to the X-ray source, possibly the BLR clouds (Maiolino 2000). In the cases of dual absorbers we assume the NH of the lower column density absorber, both to be conservative and because the higher density absorber has generally a low partial covering (Maiolino 2000). Tab.1 lists the sources in our sample.

The gaseous NH along the line of sight is derived directly by the photoelectric cutoff in the X-ray spectrum, provided that the signal-to-noise is high enough. Tab. 1 lists the NH derived for the sources in our sample along with the intrinsic 2–10 keV luminosity (i.e. absorption corrected). As mentioned above we excluded objects whose X-ray spectrum shows evidence for warm absorption as based on the presence of the absorption edges of OVII and OVIII at 0.74 and 0.87 keV respectively. For most of the objects in Tab. 1 the X-ray absorption is “cold”, in the sense that the chi-squared of the spectral fit is significantly better with the latter model than with a warm absorption model (this is discussed in some of the references reported in Tab. 1). For some objects the signal-to-noise of the X-ray spectrum is not high enough to discard warm absorption; this is certainly the case for SAX0045-25, SAX1218+29, SAX1519+65, AXJ0341-44, and Mkn231. Yet, the spectral shape in most of the latter cases is such that a fit with a warm absorber would require an unrealistically high colum of warm gas ($`10^{24}\mathrm{cm}^2`$, possibly with the exception of Mkn231) thus favoring the cold absorption model. Finally, even in those cases for which cold absorption provides a better fit with respect to a warm absorber this generally does not exclude that both components are present. Also, ionized gas along the line of sight might not be sampled by the OVII and OVIII absorption edges if the ionization stage is too low or too high (Kraemer et al. 1999, 2000, Brandt et al. 1996, Reynolds & Fabian 1995). However, the presence of ionized gas along the line of sight which is not detected in the X-ray spectra would make the total column of gas (neutral+ionized) higher than inferred assuming a single cold, component and, as we shall see, the problem of the reduced $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ would be even worst.

Ratios between broad components of the hydrogen lines compared to the intrinsic values give the amount of dust reddening affecting the BLR. However, radiative transport and collisional excitation effects in the extreme conditions of the BLR clouds ($`\mathrm{n}10^9\mathrm{cm}^3`$) can affect the standard hydrogen line ratios expected in the case B recombination. For instance, BLR models expect the H$`\alpha `$/H$`\beta `$ Balmer decrement to range from the “standard” ratio of 3.1 up to a factor of 3 higher (Rees et al. 1989, Netzer et al. 1985, Mushotzky & Ferland 1984), this is mostly due to the large optical thickness of H$`\beta `$ whose de-excitation transition has a high probability of being split into Pa$`\alpha `$ and H$`\alpha `$. Nonetheless, the H$`\alpha `$/H$`\beta `$ ratio observed in Sy1s and QSOs is often consistent with the standard case B value. Should the intrinsic ratio be higher, the observed H$`\alpha `$/H$`\beta `$ compared to the case B value provides at least an upper limit to the reddening. Ratios between infrared broad hydrogen lines (Pa$`\beta `$, Pa$`\alpha `$ and Br$`\gamma `$) provide a more reliable measure of the reddening, since they are much less affected by the radiative transport effects discussed above. In some cases the broad lines were not measured simultaneously and therefore variability might have affected the real ratios.

As a consequence of what discussed above we adopted the following criteria to select the broad line pairs to be used for the reddening determination. When more than two broad lines were available, to avoid problems related to the variability we chose those lines which were observed (nearly) simultaneously. The cases for which the broad lines were measured simultaneously are marked with a “\*” on column 6 of Tab. 1 (these cases are the majority). When more than two simultaneous broad lines were available we used the ratios involving only near-IR lines rather than H$`\alpha `$/H$`\beta `$. Finally, Pa$`\beta `$/H$`\alpha `$, Pa$`\alpha `$/H$`\alpha `$ and (in one case) H$`\beta `$/H$`\delta `$ were used only in a few cases. The line ratio adopted to estimate the reddening for each object is given on column 6 of Tab. 1.

To determine the reddening EB-V from the broad line ratios an extinction curve must be assumed. We have assumed the “standard” Galactic extinction curve (Savage & Mathis 1979). To derive the extinction we also assumed a foreground screen, which is a reasonable assumption given that the BLR is within the sublimation radius (and therefore mostly dust-free) and that the torus is probably more extended than the BLR (paper II). Column 2 of Tab. 1 lists the dust reddening toward the BLR derived for each object.

We also include in our study three AGNs (namely SAXJ0045-25, SAXJ1218+29, SAXJ1519+65) whose reddening was not inferred from the ratios of broad hydrogen lines but by a detailed fit to the continuum and of the EW of the H lines. More specifically, the optical spectrum (including shape, stellar features and broad H lines), along with near-IR to U-band spectroscopic measurements, were fitted with stellar population synthesis templates of various ages combined with an AGN template reddened by various degrees of extinction. The three free parameters required to fit the data were the relative contribution of the stellar and of the AGN component, the age of the stellar population and the reddening of the AGNs. For a more detailed description of the method see Maiolino et al. (2000a) and Vignali et al. (2000). Although the large number of observational constraints allow a relatively good determination of the three free parameters, the determination of the reddening is not as accurate as for the broad lines ratio method. However, for one object (SAXJ1353+18) both methods could be used and gave consistent results (Maiolino et al. 2000a, Vignali et al. 2000). The interesting property of these objects is that they were selected (and actually discovered) in the hard X-rays (Fiore et al. 1999) and, therefore, are not affected by some of the biases that might affect the other objects and which will be discussed later.

### 2.2 Evidence for a low $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$

On column 4 of Tab.1 we report the $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ ratio relative to the Galactic standard value of $`1.7\times 10^{22}\mathrm{mag}\mathrm{cm}^2`$ (Bohlin et al. 1978). Except for a few cases, $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ is significantly lower than the Galactic standard value, by a factor ranging from a few to $``$100. This is graphically shown in Fig.1 where the $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ ratio relative to Galactic is plotted as a function of the intrinsic X-ray luminosity. Fig.1 does not really show a correlation between the two quantities, but rather a bimodal behavior: AGNs with luminosities higher than $`10^{42}\mathrm{erg}\mathrm{s}^1`$ are systematically characterized by $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ lower than Galactic, though they show a large spread, while those few Low Luminosity AGNs (LLAGNs) in our sample, with $`\mathrm{L}_\mathrm{X}<10^{42}\mathrm{erg}\mathrm{s}^1`$, are characterized by $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ consistent with, or even higher, than Galactic. The markedly different behavior of LLAGNs with respect to the other AGNs in the sample, might reflect the fact that the physics of these objects is intrinsically different from the “classical” AGNs, as suggested by various authors (eg. Ho 1999). In the following we will focus on the other AGNs in our sample, which have luminosities more typical of classical Seyfert galaxies or of QSOs.

Admittely, our sample is not very large and is not representative of the population of obscured AGNs, given that we had to select our objects in a relatively narrow range of absorptions. In particular, the absorption must be low enough to enable us to detect the broad lines. We cannot exclude the existence of AGNs (in the normal-high luminosity range $`>10^{42}\mathrm{erg}\mathrm{s}^1`$) with an $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ consistent with Galactic. For instance, objects with $`\mathrm{N}_\mathrm{H}`$ higher than a few times $`10^{22}\mathrm{cm}^2`$ and Galactic $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ would be characterized by an extinction $`A_V`$ so high to make broad lines undetectable, possibly even in the IR, and therefore would be excluded from our sample. Nonetheless, we would like to stress the following results: 1) the existence of a population of AGNs characterized by a value of $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ significantly lower than Galactic is proven beyond any doubt; 2) with the exception of the three LLAGNs discussed above, we could not find AGNs whose $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ is consistent with the Galactic standard value.111Note that even with the bias discussed above we would have expected to find some object with Galactic $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$, at least in those cases with $`\mathrm{N}_\mathrm{H}10^{22}\mathrm{cm}^2`$ and observed in the near-IR. Determining if this is a property common to all AGNs or determining what fraction of them is characterized by this feature cannot be done with the current sample alone.

## 3 Evidences for a low AV/NH

The evidence for a low $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ discussed above might be translated into evidence for a low $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ if a standard conversion factor $`\mathrm{A}_\mathrm{V}/\mathrm{E}_{\mathrm{B}\mathrm{V}}=3.1`$ applies to these objects. Indeed, such a conversion factor might not apply and, therefore, a reduced reddening does not necessary imply a reduced absorption. However, there are other observational evidences supporting the idea that also the AV/NH ratio AGNs is generally significantly lower than the Galactic standard value.

### 3.1 The intermediate 1.8–1.9 type Seyferts

The intermediate type Seyferts discussed in the former section are only a small fraction of those available in various Seyfert samples; in particular the discussion in Sect.2 is limited to the sources for which enough information is available to derive both reddening and gaseous column density. In some of these objects the weakness of the broad H$`\beta `$ or H$`\alpha `$ lines is due to intrinsic properties of the BLR (Goodrich 1995), but more often the weakness of the broad lines lines is ascribed to dust absorption (Maiolino & Rieke 1995). In the latter case the extinction must be about $`\mathrm{A}_\mathrm{V}3`$, given that usually in these objects the faint broad lines have a flux similar to the narrow components, while in type 1 objects broad lines are about 10 times stronger than the narrow lines. On the other hand, most of the intermediate type Seyferts are characterized by an absorbing column density between $`10^{22}cm^2`$ and $`10^{23}cm^2`$ and, in some cases, even higher than $`10^{23}cm^2`$ (Risaliti et al. 1999). With a Galactic conversion factor such column densities would imply an extinction of $`\mathrm{A}_\mathrm{V}30`$, which would obviously make any broad line undetectable (unless reflected, paper II). This suggests that, at least in this class of objects, the AV/NH ratio must be about a factor of 10 lower than Galactic.

### 3.2 QSOs with hard X-ray absorption

Local Seyfert 1 galaxies are generally characterized by an X-ray spectrum with low or no cold absorption, in agreement with the extinction inferred by their optical spectra. However, at higher luminosities, there is evidence for a population of QSOs whose hard X-ray spectrum is characterized by significant cold absorption along the line of sight although their optical spectrum has prominent broad lines and, usually, also a blue continuum typical of unabsorbed AGNs.

In particular, although most of the optically selected PG QSOs do not show evidence for X-ray gaseous absorption in excess of the Galactic value, a few of them are affected by an excess absorption with $`\mathrm{N}_\mathrm{H}>10^{21}\mathrm{cm}^2`$. This is shown in the top panel of Fig. 2, where the distribution of cold absorbing columns, as inferred from the hard X-rays, is reported for the PG QSOs.

Obviously, even a small amount of dust associated with the gaseous column along the line of sight can redden the QSO continuum enough to be missed by U–B color selection criterion of the PG survey. Indeed, surveys of type 1 QSOs selected at wavelength less sensitive to dust obscuration and reddening have found a larger fraction of objects which show evidence for significant gaseous absorption. In particular, the lower two panels of Fig. 2 show the distribution of NH for the type 1 QSOs in the ASCA survey presented by Akiyama et al. (2000) and in the radio-selected type 1 AGNs presented by Sambruna et al. (1999). These distributions clearly show a larger fraction of type 1 AGNs with NH even in excess of $`10^{22}\mathrm{cm}^2`$. On the other hand, the prominent broad lines and, often, the blue continuum (although not as “blue” as for the PG QSOs) suggest that the optical absorption and reddening cannot be large. In particular, Akiyama et al. (2000) show that, although significant X-ray absorption is observed in their ASCA selected objects, the ratio between optical and hard X-ray emission is not different from soft X-ray selected QSOs (i.e. not affected by significant absorption). This indicates that optical absorption must be lower than about 1 mag222The hard X-ray flux is not significantly affected by absorption for NH values which are typical of these objects ($`\mathrm{N}_\mathrm{H}<10^{23}\mathrm{cm}^2`$, Fig. 2). and, together with the NH inferred from the X-rays, indicates an $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ about a factor of ten lower than Galactic.

The HELLAS BeppoSAX survey (Fiore et al. 1999, Comastri et al. 2000) has also found some type 1 blue QSOs with flat X-ray spectra suggesting large columns of gas along our line of sight. Although statistically less relevant than the surveys mentioned above, we could study some of these blue X-ray–absorbed QSO more in detail (eg. Maiolino et al. 2000): we found an $`\mathrm{E}_{\mathrm{B}\mathrm{V}}`$ which is a factor of $``$30 lower than expected from the X-ray absorbing NH and a Galactic $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$.

Within this context it is puzzling that, at variance with the ASCA and BeppoSAX results, Chandra has not found a large number of objects of this class so far (Fiore et al. 2000, Brandt et al. 2000b, Mushotzky et al. 2000): only a few blue type 1 QSOs with significant NH were discovered by Chandra’s surveys. Possibly this is due to the sensitivity of Chandra which peaks below 2 keV and, therefore, probably biases the results in favor of soft, little absorbed X-ray sources. XMM, whose sensitivity is much more uniform up to $`7keV`$, should tackle this issue.

### 3.3 Broad Absorption Line QSOs

Probably the most extreme case of type 1 AGNs with low $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ are the Broad Absorption Line (BAL) QSOs. Brandt et al. (2000a) and Gallagher et al. (1999) found that these objects are extremely weak both in soft X-rays and hard X-rays with respect to the optical. They suggest that this property is not due their intrisic SED, but to a large column of absorbing gas along the line of sight, probably related to the same medium seen in UV resonant absorption lines. The X-ray emission, heavily suppressed even in the 2–10 keV band, suggests that the gaseous absorbing column must be relatively high: $`10^{23}cm^2`$ or higher. Yet, the optical to UV spectrum is rather typical of normal, unabsorbed QSOs (except for the presence of the broad resonant absorption lines), implying little or no dust absorption. Combined with the large absorbing column inferred from the X-rays, this implies that these objects are characterized by an $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratio which is nearly two orders of magnitude lower than Galactic.

### 3.4 Grism selected QSOs

As mentioned above, even a small amount of dust might redden the QSO spectra enough to exclude them from surveys based on color selection criteria such as the Palomar Green survey. On the other hand, grism QSO surveys are mostly based on the detection of broad lines whose flux might be little affected by extinction in the case of a low $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$, even if a substantial amount of gas were present along the line of sight. Therefore, grism-selected QSOs are among the best suited objects to search for effects of a low $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratio. Indeed, by comparing the optical and X-ray emission of grism selected QSOs, Risaliti et al. (2000) find evidence for a significant population of QSO characterized by a low $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratio. The reader is addressed to that paper for an exhaustive discussion.

## 4 Evidences for peculiar properties of dust grains

In the former two sections we have discussed observational evidences on a reduced dust reddening and absorption towards active nuclei with respect to what would be expected from the gaseous NH, assuming a Galactic gas–to–dust ratio and extinction curve. This discrepancy can be ascribed to various effects, as discussed in detail in paper II. One possibility is that the properties of dust grains in the circumnuclear region of AGNs are different from the diffuse ISM of our Galaxy. In the following we present two observational evidences supporting this scenario.

### 4.1 The 9.7$`\mu `$m silicate feature

The silicate feature at 9.7$`\mu `$m observed in the mid-IR spectra of many Galactic sources is commonly ascribed to silicate grains with sizes smaller than $`3\mu `$m. Ground-based mid-IR studies of type 2 Seyferts claimed the detection of a deep silicate absorption feature in the spectra of many Sy2s, which was regarded as an evidence supporting the unified model. However, recent ISO spectra have shown that many of such detections were probably spoilt by the narrow bandwidth of the data, limited by the atmospheric transmission. Indeed, the presence of strong PAH features on both sides of the silicate dip prevented a reliable determination of the continuum in ground-based observations. The average ISO spectrum of a sample of Sy2s obtained by Clavel et al. (2000) does not show evidence for any silicate dip, thus questioning the case for a significant absorption feature claimed in former studies. In particular the very conservative upper limits on the EW of the silicate feature given by Clavel and collaborators for the average spectrum of Sy2s (EW$`<`$0.32$`\mu `$m) implies $`\tau _{9.7}<0.6`$.

On the other hand the mid-IR continuum produced by the active nucleus is certainly suppressed in Sy2s with respect to Sy1s. This is apparent, as discussed by Clavel et al. (2000), in the much larger equivalent width of the PAH features of Sy2s with respect to Sy1s. In particular, Clavel et al. estimate an average dust extinction in the mid-IR band and, more specifically, at 7.7$`\mu `$m of about 1.83 mag with a dispersion of $`\pm `$0.74 mag. For a Galactic dust composition this would give $`\tau _{9.7}=6.2\pm 2.5`$, i.e. one order of magnitude higher than the upper limit of 0.6 derived from the constraints on the silicate absorption feature, as discussed above.

Note that the absorptions derived above assume a uniform foreground screen model. In the reality the emitting region is extended ($``$ 10 pc) and mixed with the absorber, though most of the emission comes from the hot, inner region. As a consequence, the optical depths derived above actually gives an “equivalent” optical depth for screen absorption. Nonetheless, the relation between $`\tau _{9.7}`$ and $`\tau _{7.7}`$ does not change much with the geometry of the absorber, since we are comparing two different forms of absorption of the same continuum at nearly the same wavelength. In particular, the relation remains unchanged in the case of partial covering of the absorber (eg. patchy extinction) since the optical depth at 7.7$`\mu `$m and at 9.7$`\mu `$m are measured by means of the equivalent width of two features at nearly the same wavelength. In the case that the emitting region is mixed with the absorber the expected depth of the silicate feature is lower. If the emitting region is completely mixed with the absorber (which is an extreme case since the colder absorbing gas must be located in the outer parts with respect to the mid-IR nuclear emitting region) the “equivalent” optical depth in the silicate feature inferred from the (equivalent) $`\tau _{7.7}`$ should be $`\tau _{9.7}=3.1\pm 0.68`$ which is still much higher than the upper limit given by the non-detection of the silicate feature in abosorption. Therefore, althogh the meaning of the optical depth derived at 7.7$`\mu `$m and at 9.7$`\mu `$m depends on the geometry of the circumnuclear emitting/absorbing matter, the discrepancy between these two quantities remains significant and is not to ascribe to geometrical effects. The only concern to this regard, is that the inner hot region, which is absorbed along our line of sight by the outer colder regions, might be characterized by the silicate feature in emission, as it is observed in Sy1s (Clavel et al. 2000). Since the upper limit on $`\tau _{9.7}`$ was based on the upper limit on the equivalent width of the 9.7$`\mu `$m absorption feature assuming an intrinsic featureless continuum, the possible presence of the emission feature implies a higher upper limit on $`\tau _{9.7}`$. However, even in Sy1s the 9.7$`\mu `$m silicate emission is relatively faint (EW$`0.25\mu `$m) and, in the worst case, accounting for this emission feature would only increase the upper limit on $`\tau _{9.7}`$ by 0.35, i.e. $`\tau _{9.7}<0.95`$, which is still much lower than what was derived from the featureless absorption at 7.7$`\mu `$m ($`\tau _{9.7}=4.8`$).

This result is a very convincing indication that although dust must be absorbing the mid-IR nuclear radiation, such dust must have properties different from the Galactic diffuse interstellar medium. As discussed in paper II the most likely explanation is that dust in the circumnuclear region is predominantly composed of grains with size larger than 3$`\mu `$m, which do not contribute to the feature at 9.7$`\mu `$m.

### 4.2 The 2175Å carbon dip

The absorption feature at 2175Å observed in the diffuse interstellar medium was commonly ascribed to small graphite grains with radii $`100200`$Å. More recent studies ascribe most of the profile of this absorption feature to even smaller dust particles (PAHs, Weingartner & Draine 2000). Observing this feature in absorbed AGNs is really difficult since even a small amount of dust generally suppresses nearly completely the UV emission. On the other hand unabsorbed AGNs are generally free of any dust absorption. Yet there are a few type 1 Seyferts and QSOs whose optical and near-IR broad line ratios suggest the presence of some dust reddening along the line of sight (Lacy et al. 1982, McAlary et al. 1986, Puetter et al. 1981). Some of these objects are still relatively bright in the UV and have been observed with the HST spectrometers (FOS and STIS); these are therefore among the best suited objects to look for the carbon dip feature at 2175Å. In Fig.3 we show the UV spectra of five of these slightly reddened type 1 AGNs for which we could retrive HST archival spectra. Since we are mostly interested in the continuum shape and in the broad carbon absorption feature the spectra were smoothed to a resolution of about 1000 km/s. The thin dashed line shows the template of type 1 AGNs obtained by Francis et al. (1991) reddened with the standard Galactic extinction curve by an amount consistent with an $`\mathrm{E}_{\mathrm{B}\mathrm{V}}`$ consistent both with the reddening of the broad lines and with the shape of the UV continuum outside the carbon dip. The most important result is that the Galactic extinction curve sistematically predicts a deep feature around 2175Å which is undetected or much weaker in the observed spectra 333Note that some of the objects show a small FeII bump at 2500Å which is stronger than in the template; however, such an emission feature cannot account for the missing dip expected at 2175Å.

As discussed in paper II, this observational evidence supports the idea that small grains are depleted in the dusty medium responsible for the reddening.

## 5 Evidences for $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ higher than Galactic?

Finally, we shall discuss whether there is any evidence in other previous studies for an $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ higher than Galactic in any AGNs, besides the three LLAGNs discussed in Sect.2.

Veilleux et al. (1997) found a few AGNs whose $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ is higher than Galactic (some of which are also included in our sample). For most of these objects a high lower limit to $`\mathrm{A}_\mathrm{V}`$ was obtained by comparing a broad Pa$`\beta `$ with an upper limit on the broad component of H$`\alpha `$. Paradoxically, as acknowledged by the same authors, for the same objects the comparison between broad Pa$`\beta `$ and broad Br$`\gamma `$ gives little or no reddening. We believe that this result reflects problems with estimating the upper limits for the broad component of H$`\alpha `$ which, most probably, was underestimated. The only object in their paper for which $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ seems larger than Galactic and for which they do not use upper limits on the broad line components is NGC2992. However, the much higher quality optical and near-IR (simultaneous) spectra obtained by Gilli et al. (2000) clearly indicate an absorption significantly lower than inferred by Veilleux et al. (1997).

A significant fraction of Sy1s shows evidence for highly ionized gas along our line of sight, identified through the presence of absorption edges of OVII and OVIII in the soft X-rays, referred to as “warm absorber” (Reynolds & Fabian 1995, George et al. 1998). In some cases the X-ray spectra also show evidence for neutral gas in addition to the highly ionized gas. For some of these Sy1s with warm absorbers there is also evidence for dust reddening (Komossa & Fink 1997a,b, Leighly et al. 1997, Komossa & Bade 1998, Reynolds & Fabian 1995). The comparison between the $`\mathrm{E}_{\mathrm{B}\mathrm{V}}`$ and the column of neutral gas suggesting an $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ higher than Galactic. However, we note that there is no reason for the highly ionized gas to be dust-free if it is located at a distance larger than the sublimation radius. Indeed, the columns inferred for the highly ionized component are in agreement (or higher) with the dust reddening for a Galactic dust–to–gas ratio. Secondly, as pointed out by various authors, the ionized gas probably has multiple components at different ionization stages, which are not sampled by the OVII and OVIII edges, but which are characterized by columns high enough to account for the observed reddening (Kraemer et al. 1999, 2000, Brandt et al. 1996, Reynolds & Fabian 1995).

## 6 Conclusions

We have reported various observational evidences indicating that the dust reddening and absorption of the nuclear region of AGNs is generally much lower than the values expected from the gaseous column density measured in the X-rays, if a standard Galactic dust–to–gas ratio and extinction curve are assumed.

Quantitatively, the most convincing argument supporting a low $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ ratio is the comparison between the reddening inferred from the optical and infrared broad line ratios and the NH derived from the hard X-rays for a sample of 20 objects. With the exception of three low luminosity AGNs ($`\mathrm{L}_\mathrm{X}<10^{42}\mathrm{erg}\mathrm{s}^1`$), whose physics might differ from Seyfert– and QSO–like luminosity systems, all the objects appear characterized by an $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ ratio sistematically lower than the Galactic standard value by a factor ranging from $``$3 up to $``$100.

Also, we have presented additional evidences suggesting a reduced $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratio in AGNs. The presence of substantial gaseous absorbing column ($`\mathrm{N}_\mathrm{H}`$ a few times $`10^{22}`$) in intermediate type 1.8–1.9 Seyferts contrasts with their optical appearance, but can be reconciled by a low $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratio. A more extreme case of this effect is observed in a number of type 1 AGNs (mostly QSOs) whose X-ray spectrum shows evidence for substantial gaseous absorption despite their optical unabsorbed appearance. Two classes of type 1 AGNs which seems to be nearly systematically affected by this phenomenon (i.e. significant X-ray absorption but little, or no, optical absorption) are the Broad Absorption Line QSOs and the grism selected QSOs.

The samples used are probably affected by selection effects, which are discussed in the body of the paper, and therefore cannot be considered as representative of the whole population of AGNs. However, we can certainly state that at least a sub-population of the AGNs is characterized by a low $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ or $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ with respect to the Galactic value. Also, we could not find evidence for a Galactic standard $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ or $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ ratio in nearly any object.

We presented additional evidences indicating that the properties of dust grains in the circumnuclear region of AGNs are different with respect to the Galactic diffuse interstellar medium:
– although type 2 Seyfert nuclei appear significantly absorbed by dust in the mid–IR, their average ISO spectrum does not show evidence for the silicate absorption feature at 9.7$`\mu `$m which, instead, is expected to be very deep in case of heavy absorption;
– some type 1 Seyferts which appear affected by some reddening do not show evidence for the carbon dip at 2175Å, while according to the Galactic standard extinction curve this absorption feature should be prominent in their UV spectra.
Both these observational evidences suggest that dust in the circumnuclear region of AGNs is depleted of the small grains ($`<3\mu `$m and 100-200Å respectively) which are responsible for these absorption features. A dust grain distribution biased in favor of large grains would also make the extinction curve flatter (Laor & Draine 1993). If the bias for large grains is due to coagulation this would also explain the reduced $`\mathrm{E}_{\mathrm{B}\mathrm{V}}/\mathrm{N}_\mathrm{H}`$ and $`\mathrm{A}_\mathrm{V}/\mathrm{N}_\mathrm{H}`$ (Kim & Martin 1996). A more detailed discussion on the interpretation of these results is given in a companion paper (Maiolino et al. 2000b, paper II).

Regardless of the intepretation of these observational phenomena, these results have important consequences on the unified theories of AGNs. The finding that the dust absorption in generally significanlty lower with respect to what expected from the gaseous column density implies that for several AGNs the optical classification might be de-coupled from the X-ray classification. In particular, AGNs which appear obscured (type 2) in the X-rays might appear as relatively unobscured (type 1) in the optical.

###### Acknowledgements.

We are grateful to B. Draine for enlightening discussions during the early stages of this work. This paper also benefits of useful comments from A. Natta, M. Walmsley and from the referee R. Antonucci. This work was partially supported by the Italian Space Agency (ASI) under grant ARS–99–15 and by the Italian Ministry for University and Research (MURST) under grant Cofin98–02–32.