- About this Journal ·
- Abstracting and Indexing ·
- Advance Access ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents

Advances in High Energy Physics

Volume 2013 (2013), Article ID 627137, 15 pages

http://dx.doi.org/10.1155/2013/627137

## Viscous Coefficients of a Hot Pion Gas

Theoretical Physics Division, Variable Energy Cyclotron Centre, 1/AF, Bidhannagar, Kolkata 700064, India

Received 26 March 2013; Accepted 20 May 2013

Academic Editor: Edward Sarkisyan-Grinbaum

Copyright © 2013 Sourav Sarkar. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The steps essentially involved in the evaluation of transport coefficients in linear response theory
using Kubo formulas are to relate the defining *retarded* correlation function to the corresponding
*time-ordered* one and to evaluate the latter in the conventional perturbation expansion. Here we
evaluate the viscosities of a pion gas carrying out both the steps in the *real-time* formulation. We
also obtain the viscous coefficients by solving the relativistic transport equation in the Chapman-Enskog approximation to leading order. An in-medium cross-section is used in which spectral
modifications are introduced in the propagator of the exchanged .

#### 1. Introduction

One of the most interesting results from experiments at the Relativistic Heavy Ion Collider (RHIC) is the surprisingly large magnitude of the elliptic flow of the emitted hadrons. Viscous hydrodynamic simulations of heavy ion collisions require a rather small value of , being the coefficient of shear viscosity and the entropy density, for the theoretical interpretation of this large collective flow. The value being close to , the quantum lower bound for this quantity [1], matter produced in these collisions is believed to be almost a perfect fluid [2].

This finding has led to widespread interest in the study of nonequilibrium dynamics, especially in the microscopic evaluation of the transport coefficients of both partonic as well as hadronic forms of strongly interacting matter. In the literature one comes across basically two approaches that have been used to determine these quantities. One is the kinetic theory method in which the nonequilibrium distribution function which appears in the transport equation is expanded in terms of the gradients of the flow velocity field. The coefficients of this expansion which are related to the transport coefficients are then perturbatively determined using various approximation methods. The other approach is based on response theory in which the nonequilibrium transport coefficients are related by Kubo formulas to equilibrium correlation functions. They are then perturbatively evaluated using the techniques of thermal field theory. Alternatively, the Kubo formulas can be directly evaluated on the lattice [3] or in transport cascade simulations [4] to obtain the transport coefficients.

Thermal quantum field theory has been formulated in the imaginary as well as real-time [5–9]. For time independent quantities such as the partition function, the imaginary time formulation is well suited and stands as the only simple method of calculation. However, for time dependent quantities like two-point correlation functions, the use of this formulation requires a continuation to imaginary time and possibly back to real-time at the end. On the other hand, the real-time formulation provides a convenient framework to calculate such quantities, without requiring any such continuation at all.

A difficulty with the real-time formulation is, however, that all two-point functions take the form of matrices. But this difficulty is only apparent. Such matrices are always diagonalisable, and it is the 11 component of the diagonalised matrix that play the role of the single function in the imaginary time formulation. It is only in the calculation of this 11 component to higher order in perturbation that the matrix structure appears in a nontrivial way.

In the literature transport coefficients are evaluated using the imaginary time formulation [10–12]. Such a coefficient is defined by the *retarded* correlation function of the components of the energy-momentum tensor. As the conventional perturbation theory applies only to *time-ordered* correlation functions, it is first necessary to relate the two types of correlation functions using the Källen-Lehmann spectral representation [13–16]. We find this relation directly in real-time formulation. The time-ordered correlation function is then calculated also in the covariant real-time perturbative framework to finally obtain the viscosity coefficients of a pion gas.

We also calculate the viscous coefficients in a kinetic theory framework by solving the transport equation in the Chapman-Enskog approximation to the leading order. This approach being computationally more efficient [12] has been mostly used in the literature to obtain the viscous coefficients. The cross-section is a crucial dynamical input in these calculations. Scattering amplitudes evaluated using chiral perturbation theory [17, 18] to lowest order have been used in [19, 20], and unitarization improved estimates of the amplitudes were used in [21] to evaluate the shear viscosity. Phenomenological scattering cross-section using experimental phase shifts has been used in [20, 22–24] in view of the fact that the cross-section estimated from lowest order chiral perturbation theory is known to deviate from the experimental data beyond centre of mass energy of 500 MeV primarily due to the pole which dominates the cross-section in the energy region between 500–1000 MeV. All these approaches have used a vacuum cross section. To construct an in-medium cross-section we employ an effective Lagrangian approach which incorporates and meson exchange in scattering. Medium effects are then introduced in the propagator through one-loop self-energy diagrams [25].

In Section 2 we derive the spectral representations for the retarded and time-ordered correlation functions in the real-time version of thermal field theory. We also review the formulation of the nonequilibrium density operator and obtain the expressions for the viscosities in terms of equilibrium (retarded) two-point functions. The time-ordered function is then calculated to lowest order with complete propagators in the equilibrium theory. In Section 3 we briefly recapitulate the expressions for the viscosities obtained by solving the Uehling-Uhlenbeck transport equation in the kinetic theory framework. We then evaluate the cross-section in the medium briefly discussing the one-loop self-energy due to loops evaluated in the real-time formulation discussed previously. We end with a summary in Section 4.

#### 2. Viscous Coefficients in the Linear Response Theory

##### 2.1. Real-Time Formulation

In this section we review the real-time formulation of equilibrium thermal field theory leading to the spectral representations of bosonic two-point functions [16]. This formulation begins with a comparison between the time evolution operator of quantum theory and the Boltzmann weight of statistical physics, where we introduce as a complex variable. Thus while for the time evolution operator, the times and are any two points on the real line, the Boltzmann weight involves a path from to in the complex time plane. Setting this , where is real, positive, and large, we can get the contour shown in Figure 1, lying within the region of analyticity in this plane and accommodating real-time correlation functions [6, 8].

Let a general bosonic interacting field in the Heisenberg representation be denoted by , whose subscript collects the index (or indices) denoting the field component and derivatives acting on it. Although we will call its two-point function as propagator, can be an elementary field or a composite local operator. (If denotes the pion field, it will, of course, not have any index.)

The thermal expectation value of the product may be expressed as
where for any operator denotes *equilibrium* ensemble average:
Note that we have two sums in (1), one to evaluate the trace and the other to separate the field operators. They run over a complete set of states, which we choose as eigenstates of four-momentum . Using translational invariance of the field operator:
we get
Its spatial Fourier transform is
where the times , are on the contour . We now insert unity on the left of (5) in the form
(We reserve for the variable conjugate to the real-time.) Then it may be written as
where the spectral function is given by as follows:

In just the same way, we can work out the Fourier transform of as follows: with a second spectral function given by The two spectral functions are related by the KMS relation [26, 27] in momentum space, which may be obtained simply by interchanging the dummy indices , in one of and using the energy conserving -function.

We next introduce the *difference* of the two spectral functions:
and solve this identity and the KMS relation (11) for as follows:
where is the distribution-like function
In terms of the true distribution function
it may be expressed as

With the previous ingredients, we can build the spectral representations for the two types of thermal propagators. First consider the *time-ordered* one:
Using (7), (9), and (13), we see that its spatial Fourier transform is given by [6]

As , the contour of Figure 1 simplifies, reducing essentially to two parallel lines, one the real axis and the other shifted by , points which will be denoted, respectively, by subscripts 1 and , so that [8]. The propagator then consists of four pieces, which may be put in the form of a matrix. The contour ordered s may now be converted to the usual time ordered ones. If , are both on line (the real axis), the and orderings coincide, . If they are on two different lines, the ordering is definite, . Finally if they are both on line , the two orderings are opposite, .

Back to real-time, we can work out the usual temporal Fourier transform of the components of the matrix to get where the elements of the matrix are given by [16] Using relation (16), we may rewrite (20) in terms of as follows:

The matrix and hence the propagator can be diagonalised to give where and are given by Equation (22) shows that can be obtained from any of the elements of the matrix , say . Omitting the indices , we get

Looking back at the spectral functions defined by (8) and (10), we can express them as usual four-dimensional Fourier transforms of ensemble average of the operator products, so that is the Fourier transform of that of the commutator: where the time components of and are on the real axis in the -plane. Taking the spectral function for the free scalar field we see that becomes the free propagator, .

We next consider the *retarded* thermal propagator:
where again , are on the contour (Figure 1). Noting equations (7), (9), and (12) the three-dimensional Fourier transform may immediately be written as
As before we isolate the different components with real-times and take the Fourier transform with respect to real-time. Thus for the 11 component we simply have
whose temporal Fourier transform gives
This 11 component suffices for us, but we also display the complete matrix as follows:Though we deal with matrices in real-time formulation, it is the 11 component that is physical. Equations (23) and (30) then show that we can continue the *time-ordered* two-point function into the *retarded* one by simply changing the prescription as follows:
The point to note here is that for the time-ordered propagator, it is the *diagonalised* matrix and not the matrix itself, whose 11 component can be continued in a simple way.

##### 2.2. Transport Coefficients

We now use the linear response approach to arrive at expressions of the transport coefficients as integrals of retarded Green’s functions over space. We follow the method proposed by Zubarev [28], which is excellently reviewed in [10]. Here the system is supposed to be in the hydrodynamical stage where the mean free time of the constituent particles is much shorter than the relaxation time of the whole system under consideration. Thus local equilibrium will be attained quickly, while global equilibrium will be approaching gradually. Since the system is assumed to be not far from equilibrium, we may retain only linear terms in space-time gradients of thermodynamical parameters, like temperature and velocity fields. We assume the energy-momentum of the system to be conserved as follows: The nonequilibrium density matrix operator is constructed in the Heisenberg picture, where it is independent of time as follows: Following Zubarev, we construct the operator as follows: where . Here is a Lorentz invariant quantity defining the local temperature, and is the four-velocity field of the fluid The construction (35), which smooths out the oscillating terms, resembles the one used in the formal theory of scattering [28, 29] and selects out the retarded solution.

The expression (35) is actually independent of ; the time derivative is As and are finite, the right hand side of (37) goes to zero as . Also integrating (35) by parts, we get We now consider the space integral of (38). Using the energy-momentum conservation rule (33), we integrate the second term in (38) by parts and neglect the surface integrals to get where we have abbreviated the first and second terms by and , respectively. Then the nonequilibrium statistical density matrix is given by The first term in (39) characterises local equilibrium: where the second term including the thermodynamical force describes deviation from equilibrium.

In order to expand in a series in we define the function such that the boundary conditions at and correspond to the equilibrium and nonequilibrium density matrices as follows: We then differentiate with respect to to get which can be integrated to give It can be solved iteratively. Keeping up to the first order term (linear response) and setting , we get the required result

Applying this formula to the energy-momentum tensor, we get its response to the thermodynamical forces as [10] where is the correlation function to be evaluated. As the correlation is assumed to vanish as , it can be put in terms of the conventional retarded Green’s function. Omitting indices, it is with

We now use (47) to obtain the expectation value of the viscous-shear stress part of the nonequilibrium energy momentum tensor which is given by where is the equilibrium part, is the viscous-shear stress tensor, and is the heat current. Also, with a view to separate scalar, vector and tensor processes, the quantity in (39) is expanded as with , being the sound velocity. Using now the fact that the correlation function between operators of different ranks vanishes in an isotropic medium, one can write from (47) with . Following Hosoya et al. [10], we write the correlation function as where . Assuming now that changes in the thermodynamic forces are small over the correlation length of the two-point function, the factor can be taken out of the integral giving finally: where

Again, starting with the pressure on the l.h.s of (47) and following the steps as described previously, we obtain where the bulk viscosity is given in terms of a retarded correlation function by Here with the energy density and .

Recall that denotes *equilibrium* ensemble average. From now on we will drop the subscript “0” on the correlation functions.

##### 2.3. Perturbative Evaluation

Clearly the spectral forms and their interrelations derived in Section 2.1 hold also for the two-point function appearing in (56) for the shear viscosity. We begin with four-dimensional Fourier transforms. To calculate the 11 element of the retarded two-point function we consider the corresponding time-ordered one which can be calculated perturbatively. The viscous stress tensor can be extracted from the energy momentum tensor using the formula where in which denotes the pion triplet. We take the lowest-order chiral Lagrangian given by [18] The time-ordered correlator, to leading order, is then given by Wick contractions of pion fields in which is obtained as In the so-called skeleton expansion, these contractions are expressed in terms of complete propagators (see Figure 2) to get where is given by (19) and is determined by the derivatives acting on the pion fields as follows: where the pion isospin degeneracy factor .

To work out the integral in (64), it is more convenient to use given by (20) than by (21). Closing the contour in the upper or lower half -plane, we get where The imaginary part of arises from the factor as follows: where its real part is given by the principal value integrals.

Having obtained the real and imaginary parts of , we use relations similar to (24) to build the element of the diagonalised matrix: Finally can be continued to by a relation similar to (31) as follows: Note that in (69), (70), we retain the terms in the numerator to put it in a more convenient form. Change the signs of and in the first and second term, respectively. Noting relations like and , we get where

Returning to the expression (56) for , we now get the three-dimensional spatial integral of the retarded correlation function by setting in (59) and Fourier inverting with respect to as follows: This completes our use of the real-time formulation to get the required result. The integrals appearing in the expression for have been evaluated in [10, 11], which we describe in the following for completeness.

As shown in [10], the integral over , , and in (56) and (73) may be carried out trivially to give The dependence of is contained entirely in as follows: Changing the integration variables in (71) from , to and , we get where

It turns out that the integral over becomes undefined, if we try to evaluate with the free spectral function given by (26). As pointed out in [10], we have to take the spectral function for the complete propagator that includes the self-energy of the pion, leading to its finite width in the medium Note that this form of the spectral function trivially follows on replacing (where ) with in the free spectral function (26) which can be written as

Then becomes having double poles at for and also at . The integral over may now be evaluated by closing the contour in the upper/lower half plane to get where we retain only the leading (singular) term for small . In this approximation (76) gives

Proceeding analogously as mentioned above, the lowest order contribution to the bulk viscosity can be obtained as [10]

The width at different temperatures is known [30] from chiral perturbation theory. The quantity can also be interpreted as the collision frequency, and the inverse of which is the relaxation time . For collisions of the form this is given by (see, e.g., [25])
where is the cross-section. Note that the lowest order formulae for the shear and bulk viscosities obtained previously in the linear response approach *coincide* with the expressions which result from solving the transport equation in the relaxation-time approximation.

#### 3. Viscous Coefficients in the Kinetic Theory Approach

The kinetic theory approach is suitable for studying transport properties of dilute systems. Here one assumes that the system is characterized by a distribution function which gives the phase space probability density of the particles making up the fluid. Except during collisions, these (on-shell) particles are assumed to propagate classically with well-defined position, momenta, and energy. It is possible to obtain the nonequilibrium distribution function by solving the transport equation in the hydrodynamic regime by expanding the distribution function in a local equilibrium part along with nonequilibrium corrections. This expansion in terms of gradients of the velocity field is used to linearize the transport equation. The coefficients of expansion which are related to the transport coefficients satisfy linear integral equations. The standard method of solution involves the use of polynomial functions to reduce these integral equations to algebraic ones.

##### 3.1. Transport Coefficients at First Chapman-Enskog Order

The evolution of the phase space distribution of the pions is governed by the (transport) equation: where is the collision integral. For binary elastic collisions which we consider, this is given by [22] where the interaction rate is as follows: and . The factor comes from the indistinguishability of the initial state pions.

For small deviation from local equilibrium, we write, in the first Chapman-Enskog approximation, where the equilibrium distribution function (in the new notation) is given by with , , and representing the local temperature, flow velocity, and pion chemical potential, respectively. Putting (88) in (85), the deviation function is seen to satisfy where the linearized collision term Using the form of as given in (89) on the left side of (90) and eliminating time derivatives with the help of equilibrium thermodynamic laws, we arrive at [31] where , , , and indicates a space-like symmetric and traceless combination. In this equation, wherewith and . The functions are integrals over Bose functions [31] and are defined as , denoting the modified Bessel function of order . The left hand side of (90) is thus expressed in terms of thermodynamic forces with different tensorial ranks. In order to be a solution of this equation must also be a linear combination of the corresponding thermodynamic forces. It is typical to take as which on substitution into (92) and comparing coefficients of the (independent) thermodynamic forces on both sides yields the set of equations: ignoring the equation for which is related to thermal conductivity. These integral equations are to be solved to get the coefficients and . It now remains to link these to the viscous coefficients and . This is achieved by means of the dissipative part of the energy-momentum tensor resulting from the use of the nonequilibrium distribution function (88) in where Again, for a small deviation , close to equilibrium, so that only first-order derivatives contribute, and the dissipative tensor can be generally expressed in the form [32, 33] Comparing, we obtain the expressions of shear and bulk viscosity as follows: The coefficients and are perturbatively obtained from (96) by expanding in terms of orthogonal polynomials which reduces the integral equations to algebraic ones. After a tedious calculation using the Laguerre polynomial of 1/2 integral order, the first approximation to the shear and bulk viscosity comes out as where The integrals are given by where is the chemical potential of pions. The exponents in the Bose functions are given by

and the functions represent The relative angle is defined by .

Note that the differential cross-section which appears in the denominator is the dynamical input in the expressions for and . It is this quantity that we turn to in the next section.

##### 3.2. The Cross-Section with Medium Effects

The strong interaction dynamics of the pions enters the collision integrals through the cross-section. In Figure 3, we show the cross-section as a function of the centre of mass energy of scattering. The different curves are explained in the following. The filled squares referred to as experiment is a widely used resonance saturation parametrization [23, 34] of isoscalar and isovector phase shifts obtained from various empirical data involving the system. The isospin averaged differential cross-section is given by where The widths are given by and with and .

To get a handle on the dynamics, we now evaluate the cross-section involving and meson exchange processes using the interaction Lagrangian as follows: where and . In the matrix elements corresponding to -channel and exchange diagrams which appear for total isospin and 0, respectively, we introduce a decay width in the corresponding propagator. We get [25] The differential cross-section is then obtained from , where the isospin averaged amplitude is given by .

The integrated cross-section, after ignoring the contribution is shown by the dotted line (indicated by “vacuum”) in Figure 3 and is seen to agree reasonably well with the experimental cross-section up to a centre of mass energy of about 1 GeV beyond which the theoretical estimate gives higher values. We hence use the experimental cross-section beyond this energy.

After this normalisation to data, we now turn to the in-medium cross-section by introducing the effective propagator for the in the previous expressions for the matrix elements. This is obtained in terms of the self-energy by solving the Dyson equation and is given by where is the vacuum propagator for the meson and is the self-energy function obtained from one-loop diagrams shown in Figure 4. The standard procedure [35] to solve this equation in the medium is to decompose the self-energy into transverse and longitudinal components. For the case at hand the difference between these components is found to be small and is hence ignored. We work with the polarization averaged self-energy function defined as where The in-medium propagator is then written as The scattering, decay, and regeneration processes which cause a gain or loss of mesons in the medium are responsible for the imaginary part of its self-energy. The real part on the other hand modifies the position of the pole of the spectral function.

As discussed in Section 2.1, in the real-time formulation of thermal field theory the self-energy assumes a matrix structure of which the 11 component is given by where is the 11 component of the scalar propagator given by . It turns out that the self-energy function mentioned above can be obtained in terms of the 11 component through the relations [35, 36] Tensor structures associated with the two vertices and the vector propagator are included in and are available in [37] where the interactions were taken from chiral perturbation theory. It is easy to perform the integral over using suitable contours to obtain where is the Bose distribution function with arguments and . Note that this expression is a generalized form for the in-medium self-energy obtained by Weldon [38]. The subscript on in (116) corresponds to its values for , , , and , respectively. It is easy to read off the real and imaginary parts from (116). The angular integration can be carried out using the -functions in each of the four terms in the imaginary part which define the kinematically allowed regions in and where scattering, decay, and regeneration processes occur in the medium leading to the loss or gain of mesons [37]. The vector mesons , , and which appear in the loop have negative -parity and have substantial and decay widths [39]. The (polarization averaged) self-energies containing these unstable particles in the loop graphs have thus been folded with their spectral functions as follows: with . The contributions from the loops with heavy mesons may then be considered as a multipion contribution to the self-energy.

The in-medium cross-section is now obtained by using the full -propagator (113) in place of the usual vacuum propagator in the scattering amplitudes. The long dashed line in Figure 3 shows a suppression of the peak when only the loop is considered. This effect is magnified when the loops (solid line indicated by multipion) are taken into account and is also accompanied by a small shift in the peak position. Extension to the case of finite baryon density can be done using the spectral function computed in [40] where an extensive list of baryon (and antibaryon) loops are considered along with the mesons. A similar modification of the cross-section for a hot and dense system was seen also in [41].

We plot versus in Figure 5 obtained in the Chapman-Enskog approximation showing the effect of the in-medium propagation in the pion gas [25]. We observe 10% change at MeV due to medium effects compared to the vacuum when all the loops in the self-energy are considered. The effect reduces with temperature to less than at 100 MeV.

We noted in Section 2 that the lowest order result for in the response theory framework coincides with that obtained in the relaxation time approximation which is in fact the simplest way to linearize the transport equation. Here one assumes that goes over to the equilibrium distribution as a result of collisions, and this takes place over a relaxation time which is the inverse of the collision frequency defined in (84). The right hand side of (85) is then given by which subsequently leads to the expressions (82) and (83) for the shear and bulk viscosities [42]. In Figure 6 we show the temperature dependence of in the relaxation time approximation. The values in this case are lower than that obtained in the Chapman-Enskog method, though the effect of the medium is larger. In addition to the fact that the expressions for the viscosities are quite different in two approaches, the difference in the numerical values obtained in the two cases also depends significantly on the energy dependence of the cross-section [43].

In Figure 7 we show the numerical results for the bulk viscosity of a pion gas as function of . It is seen from an analysis of the left hand side of the transport equation that, while the shear viscosity depends on elastic processes, bulk viscosity is sensitive to number changing processes. However, in heavy ion collision experiments matter is known to undergo early chemical freezeout. Number changing (inelastic) processes having much larger relaxation times go out of equilibrium at this point and a temperature dependent chemical potential results for each species so as to conserve the number corresponding to the measured particle ratios. We hence use a temperature dependent pion chemical potential taken from [44] in this case. It is interesting to observe that decreases with in contrast to which increases. The trend followed by is similar to the findings of [45]. Additional discussions concerning the temperature dependence of viscosities for a chemically frozen pion gas are available in [31].

#### 4. Summary and Conclusion

To summarize, we have calculated the shear viscosity coefficient of a pion gas in the real-time version of thermal field theory. It is simpler to the imaginary version in that we do not have to continue to imaginary time at any stage of the calculation. As an element in the theory of linear response, a transport coefficient is defined in terms of a retarded thermal two-point function of the components of the energy-momentum tensor. We derive Källen-Lehmann representation for any (bosonic) two-point function of both time-ordered and retarded types to get the relation between them. Once this relation is obtained, we can calculate the retarded function in the Feynman-Dyson framework of the perturbation theory.

Clearly the method is not restricted to transport coefficients. Any linear response leads to a retarded two-point function, which can be calculated in this way. Also quadratic response formulae have been derived in the real-time formulation [46].

We have also evaluated the viscous coefficients in the kinetic theory approach to leading order in the Chapman-Enskog expansion. Here we have incorporated an in-medium cross-section and found a significant effect in the temperature dependence of the shear viscosity.

The viscous coefficients and their temperature dependence could affect the quantitative estimates of signals of heavy ion collisions particularly where hydrodynamic simulations are involved. For example, it has been argued in [47] that corrections to the freeze-out distribution due to bulk viscosity can be significant. As a result the hydrodynamic description of the spectra and elliptic flow of hadrons could be improved by including a realistic temperature dependence of the viscous coefficients. Such an evaluation essentially requires the consideration of a multicomponent gas preferably containing nucleonic degrees of freedom, so that extensions to finite baryon chemical potential can be made. Work in this direction is in progress.

#### Acknowledgments

The author gratefully acknowledges the contribution from his collaborators S. Mallik, S. Ghosh, and S. Mitra to various topics presented here.

#### References

- P. Kovtun, D. T. Son, and A. O. Starinets, “Viscosity in strongly interacting quantum field theories from black hole physics,”
*Physical Review Letters*, vol. 94, no. 11, Article ID 111601, 4 pages, 2005. View at Publisher ·