## Servicios Personalizados

## Articulo

## Indicadores

- Citado por SciELO
- Accesos

## Links relacionados

- Similares en SciELO

## Compartir

## Brazilian Journal of Physics

*versión impresa* ISSN 0103-9733

### Braz. J. Phys. vol.40 no.3 São Paulo set. 2010

#### http://dx.doi.org/10.1590/S0103-97332010000300014

**Effect of electron inertial delay on Debye sheath formation**

**U. Deka ^{I,*}; C.B. Dwivedi^{II}**

^{I}Department of Physics, Sikkim Manipal Institute of Technology, Majitar, East Sikkim-737136, Bharat (India)

^{II}Retired from Centre of Plasma Physics, Tepesia, Sonapur, Kamrup-782402, Assam, Bharat

**ABSTRACT**

Present contribution deals with the role of weak but finite electron inertia on the sheath formation condition. As reported earlier this becomes effective when the ions' drift velocity exceeds the phase velocity of the acoustic wave fluctuations. Such situation has natural existence near the sheath edge. Keeping this in mind we have revisited the problem of usual Bohm sheath condition. Analytical and numerical analysis have been performed to re-derive the local condition for plasma sheath formation. It is found that the weak but finite electron inertia reduces the threshold value of ion Mach number that may be, at least in principle, of qualitative value to define the sheath edge boundary. Consideration of finite but weak equilibrium electron flow at the defined sheath edge shrinks the width of non-neutral space charge layer over which major potential drop and charge imbalance occurs. Detailed numerical analysis and results of quantitative and qualitative importance are included in the text.

**Keywords:** Debye sheath plasma, presheath plasma, sheath edge singularity, weak electron inertia, Bohm condition.

**1. INTRODUCTION**

Natural formation of a nonlinear non-neutral space charge layer in the vicinity of a boundary wall confining the plasma is well established [1-3]. This is basically an outcome of the plasma-boundary wall interaction processes. The localized non-neutral space charge layer, thus formed near the wall boundary, is equivalent to a localized non-neutral plasma of few Debye lengths order termed as Debye sheath plasma (DSP). This is well established that the existential condition for such layer, popularly known as the Debye sheath or simply sheath, is governed by the Bohm criterion [2]. However, there are many questions related to the plasma sheath structure, which are unresolved. For example, location of Debye sheath edge is still not exactly specified and no proper and firm definition and physical and mathematical modelling of the sheath edge has been evolved.

Two-layer theory of the plasma sheath, as pointed out by Riemann in his review article [3], is inadequate due to singularity problem at the presheath termination that coincides with the Debye sheath edge. The plasma sheath as a whole comprises of bulk plasma, quasineutral presheath plasma and non-neutral Debye sheath plasma. Normally, the system of bulk plasma and quasineutral presheath plasma is treated as a single unit of quasineutral presheath plasma. But in reality the presheath plasma region is a source for ambipolar acceleration of the ions as per the requirement of existence condition of the Debye sheath plasma. The natural properties of the presheath-Debye sheath transition layer is still not well understood [4-6]. Similarly, the proper theoretical, numerical and experimental descriptions of the presheath potential distribution are also not performed [7].

Of course, numerous efforts are made to resolve these issues, at least, theoretically [8]. But good experimental measurements are not numerous and it is still lacking to complement theoretical efforts [8]. In fact, the need of more versatile and improved diagnostics is required to produce spatially resolved data on plasma sheath potential profiles with more accuracy in electrostatic potential and ion flow measurements. Laser induced fluorescence techniques [9,10] may offer a good alternative diagnostics for the purpose. In spite of all these efforts the scientific debates and discussions persist about the patching problem of presheath-Debye sheath layers and universal significance of the Bohm condition [11-13].

As pointed out by Riemann in his review article [3], the Bohm condition is basically a local condition. It is derived by retaining only linear order potential terms in Poisson equation. Subsequently, the demand for monotonous solution of DSP potential dictates the Bohm condition to arise in an inequality form of *M*_{0} > 1. Here *M*_{0}, termed as the Mach number, is defined as the ratio of the drift velocity of the ion fluid at the defined sheath entrance and the phase velocity of the ion acoustic wave. The validity of Bohm condition in equality form i.e. *M*_{0} = 1 is also recognized [12] but this can not be justified by two-layer theory of Poisson equation analysis. In reality, this is well known that the presheath scale analysis of the Poisson equation results into singularity at *M*_{0}→ 1.

Infact, this is pertinent to comment that the DSP potential confines the plasma particles electrostatically and also screens the bulk plasma from the influence of any external potential field effect. This is important to further note that any process that alters the Debye shielding behavior of plasmas may affect many plasma phenomena. There is an interesting reporting where it is mentioned that sufficiently intense laser pump wave frustrates the Debye screening in the direction of the oscillations under laser electric field [14]. As a consequence, the frequency of ion waves is increased from the ion acoustic wave frequency to ion plasma frequency, which is long understood and also observed experimentally in case of thermal, unmagnetized plasmas. Detailed discussions are found in reference [14].

Similar effect, due to weak but finite electron inertial delay, is reported to occur in an ion beam plasma system where the ions are drifting with flow velocity exceeding the phase velocity of the ion acoustic wave [15]. As a result, the fluid mode of ion acoustic wave fluctuations resonantly couples with the ion beam mode to produce relaxational growth of the acoustic fluctuations at down shifted wave frequency (Doppler effect), which other wise is not possible in the asymptotic limit of *m _{e}/m_{i }*→ 0, i.e. inertialess electrons. In the light of inertia induced frustration of the Debye shielding [15], it is therefore important to review the Bohm condition of DSP formation as because the transonic point satisfies the condition for its physical realization.

In simple fluid model approach of theoretical analysis of the *scrape-off layer, SOL, *it is emphasized that the electrons can satisfy the Boltzmann relation quite closely and yet have a substantial fluid velocity of the order of acoustic speed [16]. It is thus reasonable to judge the practical significance of the asymptotic limit of *m _{e}/m_{i }*→ 0 in deriving the result of Bohm condition for sheath formation.

A simple attempt has been made to revisit the problem of sheath formation with due consideration of finite equilibrium fluid flow of thermal electron cloud as an inertial correction. Hydrodynamical plasma model is used to carry out analytical and numerical calculations to report the impact of weak but finite electron inertial delay effect on Bohm condition. It is found that the local threshold value of the Mach number for DSP formation is slightly reduced. Here the local threshold value of *M*_{0} has usual meaning in the sense that it satisfies the condition for monotonous solution of Poisson equation, when specified at a defined Debye sheath edge, where the consideration of only linear order space charge imbalance is valid. Numerical calculations are also carried out to depict the same. The DSP width is found to be decreased in comparison to that of the ideal case of inertialess electrons.

The body of the research paper is organized in the following form. Section 2 includes physical discussions of the plasma model and basic equations. Subsection 2.1 contains the basic equations. Brief review of the asymptotic result of Bohm condition in the limit of *m _{e}/m_{i }*→ 0 is included in subsection 2.2. Subsection 2.3 incorporates the calculation of sheath formation condition with finite but weak role of electron inertia. Numerical analysis of the Poisson equation and discussions of the analysis with and without electron inertia as a correction is reported in section 3. Results and conclusions are compiled in Section 4.

**2. PLASMA MODEL AND BASIC EQUATIONS**

Let us consider a simple two-component collisionless, source free and unmagnetized non-isothermal plasma. This is confined by a boundary wall that absorbs the plasma particles through surface recombination process. The boundary wall is assumed to be non-emitting. The excess thermal flux loss of the electrons at the wall then creates a net loss of negative charge in the bulk plasma and thus causes self-consistent negative biasing to the wall. Subsequently, the ions respond dynamically to this wall potential biasing, which ultimately evolves in time to a floating potential value producing thereby a localized steady state nonlinear and non-neutral space charge layer near the wall.

The thickness of non-neutral space charge layer, thus formed, is of the order of Debye scale length. This requires the Bohm condition to be fulfilled at its entrance point termed as the Debye sheath edge supposed to lie at presheath terminal point. This is also termed as the transonic point. Before we look into the role of weak but finite electron inertial delay effect on the Bohm condition, let us first refresh about the basic physical idea of the usual DSP with inertialess electrons.

**2.1 Basic equations**

For physical discussions of the problem we use simple fluid model of plasma particles. This is a widely used mathematical model approach and yields quite satisfactory results of practical importance. Of course, there are many practical situations where the kinetic approach may be more appropriate to discern new physical ideas, which otherwise can not be possible by fluid approach. Let us start with Poisson equation that describes the electrostatic field potential distribution due to continuous charge distribution in plasmas.

The normalized form of the Poisson equation can be written as,

Here the normalized differential operator ∇ is defined as ∇ = *d/d*ξ and ξ = x/λ_{De} is a normalized one dimensional space coordinate in plane geometry approximation. The notation ψ = –*e*Φ/*T _{e}* (with Φ denoting for unnormalized plasma sheath potential) corresponds to the normalized plasma sheath potential. The functional notations

*N*(ψ) =

_{e}*n*(ψ)/

_{e}*n*

_{0}and

*N*(ψ)/

_{i}= n_{i}*n*

_{0}, respectively denote for the normalized electron and ion fluid densities and

*n*

_{0}is the bulk plasma density. The notations

*e*and

*T*denote for the electronic Coulomb charge and electron temperature (in eV) respectively and λ

_{e}_{De}is the Debye length.

Let us now expand the unknown functions *N _{e}*(ψ) and

*N*(ψ), respectively, in the form of Taylor series expansion around the local plasma potential of arbitrary value ψ

_{i}_{0}, defined at some arbitrary position ξ = ξ

_{0}of presheath-Debye sheath interface in the transonic layer. Thus,

Where,

From asymptotic theory of electron dynamics, we know that the inertialess electrons obey Boltzmann distribution of density as *N _{e}*(ψ) = exp{–(ψ – ψ

_{0})} in the limit of

*m*→ 0. This yields

_{e}/m_{i }*a*

_{0e}= –1 and

*b*

_{0e}= +1. This implies that for Boltzmannian electrons the ratio of density and potential gradients in space coordinate are equal in magnitude but differ in sign depending on the even (positive) or odd (negative) natures of the derivatives. Similarly, the ion density distribution in normalized form can be written as [17],

This leads to *a*_{0i} = –1/, *b*_{0i} = 3/. Here, the different derivatives of ion density and potential depend on *M*_{0}, and these not only differ in sign but also differ in magnitude from each other.

**2.2 Debye sheath condition for inertialess electrons**

Let us now look at the Poisson equation (1) that can be simplified to read as,

On integrating this equation once one gets,

Now, from the Poisson eq. (5) and energy integral eq. (6), it is clear that the first term i.e. leading order term on r.h.s. is a linear term and yields monotonous solution if and only if *M*_{0} > 1, otherwise an oscillatory solution is obtained. This is where the logic of locality condition of the Bohm criterion lies. It is derived by considering the linear order potential terms only in the electron and ion density distributions to assess the nature (monotonous or non-monotonous) of spatial variation of charge imbalance in close vicinity of the transonic point.

Furthermore, on presheath scale normalization the same energy integral eq. (6) can be written as,

where, η = *x/L* = ξ(λ* _{De}/L*). From this scaling transformation, it is obvious to notice that on presheath scale, the Debye sheath scale length becomes a singular point as because the two-scale theory demands the validity of asymptotic limit of λ

_{De}

*/L*→ 0. Similarly, first order derivatives of plasma potential on presheath and Debye sheath scales are correlated as

*d*ψ/

*d*η = (

*d*ψ/

*d*ξ)(

*d*ξ/

*d*η) = (

*L*/λ

*)(*

_{De}*d*ψ/

*d*ξ). Again, it is clear to notice that the potential gradient on Debye sheath scale becomes a sharp gradient on presheath scale. Now, near the presheath termination point coinciding with the Debye sheath edge location at ξ = ξ

_{0}, where the plasma potential ψ ~ ψ

_{0}, the linear order term in eq. (7) demands that

*M*

_{0}= 1 as because

*a*

_{0e}

*= a*

_{0i}has to hold good under plasma approximation.

This is important to point out that a finite but small amount of charge imbalance may still exist at *M*_{0} = 1 due to lowest order nonlinear terms in eqs. (6) and (7) both. It is thus natural to find out the nature of local nonlinear potential solution of eqs. (6) and (7) keeping both the linear and nonlinear terms. Let us define the terms (1 – 1/) = *a*_{0} and (1/3)(–1+3/) = *b*_{0}.

For simple academic interest of physical significance of lowest order nonlinearity, eqs. (6) and (7) can be shown to admit simple soliton solutions of homogenous equilibrium plasmas of sheath and presheath scales. On the Debye sheath scale, the following soliton solution of eq. (6) can be derived;

and the width of the soliton will be, 2/ The soliton is formed for *M*_{0} > 1. Similarly, on the presheath scale the following solution of eq. (7) can be derived;

such that the width of the soliton will be, 2λ_{De}/.

It is now possible to argue that the linear order local solutions of the Poisson eqs. (6) and (7) on Debye sheath and presheath scales yield monotonic potential variations for *M*_{0} > 1. However, the consideration of lowest order nonlinearity yields a localized soliton solutions for *M*_{0} > 1. It is thus reasonable to comment that the potential variation from presheath scale to Debye sheath scale may undergo an asymmetric solitary type potential transition at transonic point [18]. On presheath scale calculation, the same soliton is found to have sharp width due to scale transformation. The free energy for the electron inertia induced excitation mechanism of the nonlinear normal mode of the ion acoustic wave near the transonic zone over a finite extension is because of the transonic flow itself during the transient time [19-22].

Thus the entire DSP scale potential distribution could be treated as an equivalent to the appropriate piling up of localized solitary pulses with successively increasing amplitudes such that the global appearance looks like a smooth profile. The minimum possible amplitude of the solitary pulse at the main entrance of the Debye sheath must be equal to that of KdV soliton. Of course, these speculative ideas can be tested by simulation work. Such possibilities are mentioned in simulation work [23]. Let us now switch over to the main part of the problem.

**2.3 Debye sheath condition with inertial delay effect**

In this section, we will consider a weak but finite deviation from the exact hydrostatic balancing of the thermal electron fluid in equilibrium plasma sheath potential. For this purpose the electron continuity equation is solved to yield *N _{e}V_{e} = M_{e}*

_{0}, for weak but finite non-zero bulk flow of thermal electrons whose zeroth order density distribution is described by Boltzmannian distribution. Now, the finite but weak bulk fluid flow of thermal electrons can be estimated as,

The normalized momentum equation for the thermal electron fluid is given as below

Now, once integrating eq. (11) and then substituting the value of *V _{e}* from eq. (10), one can derive the modified electron density distribution as follows,

This is the expression for electron density distribution with finite but weak electron inertial drag effect. Now, the coefficients *a*_{0e}* , b*_{0e} can be modified to become as = –(1+2(*m _{e}/m_{i}*)) , = (1+4) and

*a*

_{0i}

*, b*

_{0i}will remain the same. Hence, the local Bohm condition derived for inertialess electrons will be modified as,

This further implies that

Since the correction term introduced by the finite but weak electron inertial drag effect to the Debye shielding is small, one can binomially expand the numerator and thus find that,

This implies that the electron inertia induced frustration of the Debye shielding reduces the Mach threshold (*M*_{0} < 1) for the DSP formation by an amount of δ*M*_{0} ~ (*m _{e}/m_{i}*). Although, this correction is very nominal but it may have qualitative value so far as the idea of wave turbulence hypothesis of transonic zone is concerned [24]. Then a natural question may arise about, how this may be possible?

For example, let us assume that the sheath edge is susceptible to weakly dispersive mode of acoustic fluctuations due to electron inertia induced instability [15] as one of the possible turbulence mechanisms. This has less threshold value than that of the dispersionless acoustic speed and the reduction in excitation threshold is proportional to for weak dispersion effect [15]. Let us further argue that the reduction in Mach number threshold for Debye sheath formation i.e. δ*M*_{0} measures the source strength of charge imbalance near the sheath edge. Now, the local balancing of weak dispersion and weak nonlinearity dictates to infer that the weak turbulence scale size λ_{turb} ≈ 2πλ_{De}(*m _{i}/m_{e}*)

^{1/2}/

*M*

_{e}_{0}.

As elaborated in section 2.2., the coefficients of linear and nonlinear terms of eqs. (6) and (7) can be rederived by introducing the corrections of electron inertia induced modification in Boltzmannian distribution to read as, .

Finally, one can directly write down the soliton solutions on sheath and presheath scales as already mentioned. Some minor quantitative changes in amplitude and width will be obtained for any given Mach value. These solutions correspond to weak nonlinearity of potential perturbations over a quasineutral homogenous equilibrium plasma background and are given just for an academic interest to apprise the significance of lowest order nonlinearity.

**3. NUMERICAL ANALYSIS AND DISCUSSIONS**

For numerical appreciation of the analytical results and discussions in previous sections, following nonlinear Poisson equation is exactly solved by Runge-Kutta method [25],

Numerical solution of this equation has been carried out to depict the impact of electron inertia induced frustration of Debye shielding on the usual Bohm condition and Debye sheath potential distribution.

This is to clarify that the potential (ψ_{0}) at the assumed sheath edge position (ξ_{0}) is considered to be zero. Finite but small non-zero electric field (potential gradient) value of the order of 10^{–4} is considered at the assumed sheath edge position for numerical calculation to proceed. These values are considered for numerical solution of eq. (16) as an initial value problem. The spatial evolution of the potential is allowed to proceed till the floating boundary condition is reached. The minimum threshold of *M*_{0}, to ensure the monotonous plasma potential solution to exist, comes out to be around 0.9999 by numerical estimation under the limit of *m _{e}/m_{i }*→ 0. Now, the with the inclusion of finite but weak electron inertia (eq. 16) under the same specified initial conditions, the values of minimum threshold of

*M*

_{0}is slightly reduced for monotonous plasma sheath potential to occur which are shown in the the potential profiles in figure 1. Corresponding profiles of the electron-ion densities and ion flow velocity are shown in figures 2 and 3 respectively.

The change in minimum threshold of *M*_{0} is more clearly visible in figure 3. This reduction depends on the arbitrarily chosen value of equilibrium electron flow Mach number, *M _{e}*

_{0}at the considered sheath edge point for given mass ratio,

*m*. The chosen value of

_{e}/m_{i}*M*

_{e}_{0}should satisfy the condition

*M*

_{e}_{0}<< (

*m*)

_{i}/m_{e}^{1/2}[16]. These profiles are plotted for hydrogen and argon plasmas. This is to note that the non-neutral region with significant potential drop and charge imbalance is of the order of 20λ

*leaving a good portion of the profiles of quasineutral region.*

_{De} A few more numerical plots of the normalized potential (figure 4), density (figure 5) and ion flow velocity (figure 6) profiles are shown for a typical value of *M*_{0} = 1, for argon plasma with different arbitrary choices of the values of *M _{e}*

_{0}= 0,1.0,1.5,2.5 and 3.0. From these plots it is noticed that the width of the non-neutral space charge layer is reduced for finite but higher non-zero values of

*M*

_{e}_{0}. Moreover, the extension of quasineutral region is also reduced for finite but higher non-zero values of

*M*

_{e}_{0}. This implies that the sheath edge location shifts towards the wall by consideration of weak but finite electron inertial flow at the defined entrance point of sheath edge boundary.

This is to point out that the choice of the initial values of the requisite physical variables is arbitrarily considered to specify the entry point of the sheath edge boundary subject to the condition of existence of a monotonous potential solution. However, there seems to be a quite good portion of the profiles where the quasineutrality holds good even on the Debye sheath scale analysis of the Poisson eq. (16) with and without the role of electron inertia. Now, a pertinent question comes into our mind to think of how to appropriately define the sheath edge and its location? Could we qualify the sheath edge as a diffuse boundary with finite extension over which the quasineutrality holds good to a quite satisfactory level of charge imbalance? Furthermore, could we say that the initial conditions for monotonous potential profiles are the sheath formation conditions? These are the nontrivial issues to be resolved by long term research program on plasma sheath physics by proper consideration of wave turbulence energy in the mathematical formulation of plasma sheath analysis.

This is also seen that localized non-neutral space charge layer is formed even for *M*_{0} __<__ 1 on Debye sheath scale analysis of the Poisson equation. This is not surprising as because the nonlinear terms contribute to the finite level of charge imbalance to occur near the boundary wall. This is to remind the readers that in all these calculations the wall boundary potential is specified by imposing the floating wall conditions corresponding to inertialess Boltzmannian electrons.

**4. RESULTS AND CONCLUSIONS**

From numerical analysis of eq. (16) and discussions of the obtained numerical plots of the plasma sheath potential profile, one finds that the consideration of finite but weak electron inertia affects the threshold values of *M*_{0} for different values of *M _{e}*

_{0}and ionic mass of the chosen plasma systems. The non-neutral region of the potential profiles where the major potential drop and corresponding charge imbalance occurs shrinks for higher non-zero values of

*M*

_{e}_{0}. This means that the consideration of finite but weak electron inertia effect seems to play an important role to specify and define the nature of sheath edge boundary.

By this simple calculation we wish to conclude that more satisfactory definition and property of the sheath edge boundary is to be evolved as because the monotonous plasma potential profiles are obtained even for *M*_{0} __<__ 1. Hence the need for more consistent calculations of sheath-presheath scale analysis of Poisson equation in totality is realized to further review the universal significance of the usual Bohm condition. However, the role of electron inertia cannot be underestimated in any such endeavors of plasma sheath analysis. Since the sheath edge boundary is quite susceptible to the onset condition of acoustic wave turbulence [15, 18], the role of wave turbulence energy in any mathematical formulation of the plasma sheath analysis can not be ignored.

**Acknowledgement**

The first author acknowledges the M.S.T, D.S.T., Govt. of India for providing the fund under the O.Y.S. Fast Track Scheme for the Young Scientist vide sanction letter No. SR/FTP/PS-69/2007.

[1] L. Tonks, and I. Langmuir, Phys. Rev. **34**, 876 (1929). [ Links ]

[2] D. Bohm, The Characteristics of Electrical Discharges in Magnetic Fields, edited by A. Guthrie and R. Wakerling (New York: McGraw-Hill) (1949). [ Links ]

[3] K. -U. Riemann, J. Phys. D Appl. Phys. **24**, 493 (1991). [ Links ]

[4] R. N. Franklin, and J. Snell, Phys. Plasmas **8**, 643 (2001). [ Links ]

[5] R. N. Franklin, IEEE Trans.Plasmam Sci. **30**, 352 (2002). [ Links ]

[6] I. D. Kaganovich, Phys. Plasmas **9**, 4788 (2002). [ Links ]

[7] M. Goswami, and H. Ramachandran, Phys. Plasmas **6**, 4522 (1999). [ Links ]

[8] L. Oksuz, and N. Hershkowwitz, Phys. Rev. Lett. **89**, 145001-1 (2002). [ Links ]

[9] M. J. Goeckner, and J. Goree, J. Vac. Sci. Technol. **A7**, 977 (1989). [ Links ]

[10] M. J. Goeckner, J. Goree, and T. E. Sheridan, Phys. Fluids **B4**, 1663 (1992). [ Links ]

[11] H.-B. Valentini, Phys. Plasmas **3**, 1459 (1996). [ Links ]

[12] K.-U. Riemann, and P. Meyer, Phys. Plasmas **3**, 4751 (1996). [ Links ]

[13] X. P. Chen, Phys. Plasmas **5**, 804 (1998). [ Links ]

[14] R. P. Drake, and R. S. Marjoribanks, Phys. Plasmas **9**, 267 (2002). [ Links ]

[15] C. B. Dwivedi, and R. Prakash, J. Appl. Phys. **90**, 3200 (2001). [ Links ]

[16] P. C. Stangeby, The Plasma Boundary of Magnetic Fusion Devices (London: Institute of Physics Publishing Ltd) p 33 (2000). [ Links ]

[17] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, vol 1: Plasma Physics, 2nd edn. (New York: Plenum Press) ch 8 (1984). [ Links ]

[18] U. Deka, R. Prakash, A. Sarma, P. K. Karmakar, and C. B. Dwivedi, Phys. Scr. **69**, 303 (2004). [ Links ]

[19] P. K. Karmakar, U. Deka, and C. B. Dwivedi, Phys. Plasmas **12**, 032105 (2005). [ Links ]

[20] P. K. Karmakar, U. Deka, and C. B. Dwivedi, Phys. Plasmas **13**, 104702 (2006). [ Links ]

[21] P. K. Karmakar, Pramana-J. Phys. **68**, 631 (2007). [ Links ]

[22] P. K. Karmakar, J. Phys.: Conf. Ser. **208**, 012059 (2010). [ Links ]

[23] Y. Matsunaga, T. Matori, and T. Kato, Proc. of ICPP and 25^{th} EPS Conf. on Contr. Fusion and Plasma Physics, Praha, p 2403 (1998). [ Links ]

[24] U. Deka, C. B. Dwivedi, and H. Ramachandran, Phys. Scr. **73**, 87 (2006). [ Links ]

[25] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran, (New Delhi: Cambridge University Press, 1st Indian ed.) ch 16 (1992). [ Links ]

(Received on 19 March, 2010)

* Electronic address: udeka1@rediffmail.com