text:
# Dephasing time of disordered two-dimensional electron gas in modulated magnetic fields

## I Introduction

Dephasing is one of the key elements in the study of quantum coherent phenomena in mesoscopic systems. Coupling with environment suppresses the quantum interference of electrons. The phase breaking time, or dephasing time is the characteristic time beyond which the phase coherence is lost. Although a static magnetic field does not destroy all quantum effects, it may introduce a cutoff to the interference effect. Dephasing of disordered electron gas by coupling with a uniform magnetic field has been studied in the early stage of weak-localization theory. In the absence of spin-orbit scattering, the magnetic field suppresses the weak-localization effect and leads to a positive magnetoconductance, which has been observed(for a review, see ). The dephasing rate due to coupling with uniform field turns out to be proportional to the field amplitude. These results can also be established by qualitative considerations according to Khmelnitskii.

In recent years there has been increasing interest in hybrid semiconductor systems both for fundamental understanding and for the potentiality of enhancing the functionality of the devices. The disposition of superconducting or magnetic microstructures on the surface of heterostructure with a two-dimensional electron gas(2DEG) may produce inhomogeneous magnetic field which influences electron motion locally. Some interesting consequences of the modulated magnetic fields have already been reported. It was shown that the effect of electron-electron Umklapp scattering can be observed in a 2DEG at a GaAs/Ga1-xAlxAs interface by imposing a spatially alternating magnetic field normal to the 2DEG plane. Kubrak et al fabricated a few different types of hybrid ferromagnetic-semiconductor devices, which allow them to study how these different modulated magnetic fields influence the transport properties of 2DEG. And magnetoresistance oscillations due to internal Landau band structure of a 2D electron system in periodic magnetic field are observed. Theoretically, Peng calculated the transport properties in a parabolic channel exposed to a periodically modulated magnetic field. Gumbs and Zhang developed a magneto-transport theory for the magnetoconductivity of a square lattice in a periodically modulated magnetic field, and predicted some anomalies due to commensurability effects. Recently Matulis and Peeters studied the semiclassical magnetoresistance in weakly modulated magnetic fields. They considered the case when the field is periodic in one dimension and with zero mean. In the limit of small magnetic field amplitude, it is shown that the contribution of the magnetic modulation to magnetoresistance increases as $`B^{3/2}`$ in the diffusive limit, while increases linearly in $`B`$ in the ballistic limit.

The magnetoresistance in the diffusive limit is well described by standard weak-localization theory. It arises from the suppression of magnetic field on the cooperon propagator which represents the interference of time-reversal trajectories. Taking the magnetic field into consideration, the cooperon propagator satisfies,

$$\left[i\omega +D\left(i\frac{2e}{c}𝐀(𝐫)\right)^2\right]C(𝐫,𝐫^{};\omega )=\frac{1}{\tau }\delta (𝐫𝐫^{}),$$
(1)
which may be viewed as an imaginary-time Schr$`\ddot{\mathrm{o}}`$dinger equation with parameters properly replaced. For the uniform field situation, the solution of the corresponding Schr$`\ddot{\mathrm{o}}`$dinger equation is the well-known Landau levels. The dephasing time was found to be

$$\frac{1}{\tau _\varphi ^u}=\frac{4DeH}{\mathrm{}c},$$
(2)
where the superscript $`u`$ represents result for uniform field.

In this paper we first study the dephasing time of disordered 2DEG due to coupling with the following modulated magnetic field

$$𝐇=(0,0,H/\mathrm{cosh}^2((xx_0)/\delta )).$$
(3)
The corresponding Schr$`\ddot{\mathrm{o}}`$dinger equation has been solved exactly by Hud$`\stackrel{´}{\mathrm{a}}`$k, where it was shown that a proper transformation converts the equation into a form solved earlier by Morsh and Feshbach. Making use of this solution, we construct the cooperon propagator and calculate the dephasing time of a disordered 2DEG in this modulated field. We find in the weak inhomogeneity limit, $`\delta >>L`$, where $`L`$ is the linear size of the sample, $`\tau _\varphi ^1`$ is proportional to $`H`$, as in the uniform field case. In the strong inhomogeneity limit, it is shown that

$$\frac{1}{\tau _\varphi }=D\left(\frac{e}{\mathrm{}c}\right)^2H^2\delta ^2.$$
(4)
In the intermediate regime, a crossover between these two limits is expected to occur at $`H_c=\frac{\mathrm{}c}{4e}\delta ^2`$. The dephasing rate dependence on the magnetic field amplitude is shown to be related to the nature of corresponding single particle motion. Bound states lead to linear dependence, while nearly free motion results in quadratic. A semiclassical Monte Carlo algorithm is developed to study the dephasing time, which is of qualitative nature but efficient in uncovering the dependence of $`\tau _\varphi `$ on $`H`$ for arbitrarily complicated magnetic field modulation. Computer simulations support analytical results. The considerations are then generalized to situation where the modulated magnetic field is periodic in one direction with zero mean, and it is argued that this crossover between linear and quadratic dependence can be expected for a large class of modulated magnetic fields.

Before going to quantitative calculations, let’s see how qualitative considerations of Khmelnitskii can predict Eq.(4). Let the magnetic field be nonzero only in a stripe of width $`\delta `$. Consider a large loop of area $`D\tau `$, we are interested in the flux piercing it. When $`\delta `$ is small, the effective area where the field is nonzero is $`\delta \sqrt{D\tau }`$. Requiring the phase change due to this loop to be of order 1, $`\frac{e}{\mathrm{}c}H\delta \sqrt{D\tau _\varphi }1`$, one immediately obtains Eq.(4). In subsequent quantitative calculations, we will see that understanding of the origin of this dependence allows us to generalize the result to more complicated situations.

This paper is organized as follows. In section II, the solution of the Schr$`\ddot{\mathrm{o}}`$dinger equation by Hud$`\stackrel{´}{\mathrm{a}}`$kis briefly outlined. The construction of cooperon propagator and calculation of dephasing time is presented in Section III. Section IV contains a description of the numerical algorithm and the simulation results. The generalization to situation with modulated magnetic field periodic in one direction is presented in Sec.V. The conclusion is given in Section VI.

## II Solution of the Schr$`\ddot{\mathrm{o}}`$dinger equation in $`𝐇=(0,0,H/\mathrm{cosh}^2((xx_0)/\delta ))`$

For completeness, the solution of the Schro$`\ddot{\mathrm{o}}`$dinger equation in the modulated magnetic field $`𝐇=(0,0,H/\mathrm{cosh}^2((xx_0)/\delta ))`$ is briefly outlined in this section. Interested readers are referred to the original article of Hud$`\stackrel{´}{\mathrm{a}}`$k and the book of Morsh and Feshbach.

Under the Landau gauge

$$𝐀=(0,H\delta \mathrm{tanh}((xx_0)/\delta ),0),$$
(5)
the Schr$`\ddot{\mathrm{o}}`$dinger equation can be written as

$`{\displaystyle \frac{1}{2m}}\left[\mathrm{}^2{\displaystyle \frac{d^2}{dx^2}}+\left(p_y{\displaystyle \frac{eH\delta }{c}}\mathrm{tanh}({\displaystyle \frac{xx_0}{\delta }})\right)^2\right]\chi (x)`$ (6)
$`=E\chi (x),`$ (7)
where we have separated variables

$$\psi (x,y)=e^{\frac{i}{\mathrm{}}p_yy}\chi (x).$$
(8)
This Schr$`\ddot{\mathrm{o}}`$dinger equation describes the motion of a particle in the potential $`V(x)`$,

$$V(x)=\frac{\mathrm{}^2}{2m\delta ^2}\left(2\pi \frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}\right)^2\left(\frac{p_y\delta /\mathrm{}}{2\pi \mathrm{\Phi }/\mathrm{\Phi }_0}\mathrm{tanh}\left(\frac{xx_0}{\delta }\right)\right)^2$$
(9)
where $`\mathrm{\Phi }_0=hc/e`$ is the flux quantum, $`\mathrm{\Phi }=H\delta ^2`$ is a measure of the magnetic flux by the external field. Introducing

$`F`$ $`=`$ $`2\pi {\displaystyle \frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}},`$ (10)
$`P`$ $`=`$ $`p_y\delta /\mathrm{},`$ (11)
$`z`$ $`=`$ $`(xx_0)/\delta ,`$ (12)
the potential may be written as

$$V(z)=\frac{\mathrm{}^2F^2}{2m\delta ^2}\left(\frac{P}{F}+\mathrm{tanh}z\right)^2.$$
(13)
For $`p_y0`$, this potential is an asymmetric well, with different limiting value of $`V(+\mathrm{})`$ and $`V(\mathrm{})`$. Hud$`\stackrel{´}{\mathrm{a}}`$k observed that the following transformation

$`\nu \mathrm{cosh}^22\mu `$ $`=`$ $`P^2+F^2,`$ (14)
$`\nu \mathrm{sinh}^22\mu `$ $`=`$ $`2FP,`$ (15)
can convert the Schr$`\ddot{\mathrm{o}}`$dinger equation into a form that was solved by Morsh and Feshbach earlier,

$`{\displaystyle \frac{d^2\chi (z)}{dz^2}}+(\epsilon \nu \mathrm{cosh}2\mu \nu \mathrm{sinh}2\mu \mathrm{tanh}z`$ (16)
$`+\nu {\displaystyle \frac{\mathrm{cosh}^2\mu }{\mathrm{cosh}^2z}})\chi (z)=0,`$ (17)
where $`\epsilon =2m\delta ^2E/\mathrm{}^2`$.

One may distinguish situations between $`|P|F`$ and $`|P|>F`$. In the former case, there is a discrete part in the energy spectrum as well as a continuous one, while the latter leads to only continuous part. The solution is as follows. If

$$|P|<F(F|P|)^2,$$
(18)
then the energy spectrum for the motion in the $`x`$-direction contains a discrete part given by

$`E_n(p_y)={\displaystyle \frac{p_y^2}{2m}}\left[1{\displaystyle \frac{F^2}{\left(\sqrt{F^2+\frac{1}{4}}(n+\frac{1}{2})\right)^2}}\right]`$ (19)
$`+{\displaystyle \frac{\mathrm{}^2}{2m\delta ^2}}\left[F^2\left(\sqrt{F^2+{\displaystyle \frac{1}{4}}}(n+{\displaystyle \frac{1}{2}})\right)^2\right],`$ (20)
with $`n=0,1,\mathrm{},[\mathrm{n}_{\mathrm{max}}]`$, and

$$n_{max}=\sqrt{F^2+\frac{1}{4}}\frac{1}{2}\sqrt{|P|F}.$$
(21)
The corresponding eigenfunction is

$`\chi _n(z)=N_n\mathrm{exp}(a_nz)(e^z+e^z)^{b_n}\times F(n,`$ (22)
$`2\sqrt{F^2+1/4}n,a_n+b_n+1,e^z/(e^z+e^z)),`$ (23)
where

$$a_n=\frac{|P|F}{\sqrt{F^2+\frac{1}{4}}(n+\frac{1}{2})},$$

$$b_n=\sqrt{F^2+\frac{1}{4}}(n+\frac{1}{2}).$$

Here $`F(\alpha ,\beta ,\gamma ,\delta )`$ denotes the hypergeometric function, and $`N_n`$ is a normalization constant.

Thus there is discrete part if

$$|P|<P_d(F)=F+\frac{1}{2F}\sqrt{1+\frac{1}{4F^2}}.$$
(24)
The continuous spectrum may be divided into two parts. For

$$(F|P|)^2<\epsilon <(F+|P|)^2,$$
(25)
the energy is given by

$$E(k,p_y)=\frac{1}{2m}\left(\left(\frac{\mathrm{}k}{\delta }\right)^2+(|p_y|p(F))^2\right).$$
(26)
where $`k`$ is the momentum in $`x`$ direction, $`p(F)=F\mathrm{}/\delta `$. The corresponding wavefunction is

$`\chi _{k,p_y}^A(z)=N_{k,p_y}\mathrm{exp}(az)(e^z+e^z)^b\times `$ (27)
$`F(b\gamma +{\displaystyle \frac{1}{2}},b+\gamma +{\displaystyle \frac{1}{2}},K_++1,e^z/(e^z+e^z)),`$ (28)
with

$$a=(k_++ik)/2,b=(k_+ik)/2,$$

$$\gamma =\sqrt{F^2+\frac{1}{4}},k_+=\sqrt{4F|P|k^2}.$$

This wavefunction vanishes exponentially for $`x\mathrm{}`$ when $`p_y<0`$, and bounded from above when $`x\mathrm{}`$. For $`p_y>0`$, the opposite is true. The inequality (25) is equivalent to

$$0<k^2<4|P|F.$$
(29)
If the parameters are such that the inequality (29) is replaced by

$$k^24|P|F,$$
(30)
then there is another type of continuous energy spectrum. The wavefunctions are nonvanishing but bounded as $`|x|\mathrm{}`$. They describe overbarrier motion. The energy is still given by Eq.(26), but the corresponding wavefunctions are

$`\chi _{k,p_y}^B(z)=N_{k,p_y}\mathrm{exp}(iz(k_+k)/2)(e^z+e^z)^{i\frac{k+k_+}{2}}`$ (31)
$`\times F((i(k+k_+)+12\gamma )/2,(i(k+k_+)`$ (32)
$`+1+2\gamma )/2,1ik_+,e^z/(e^z+e^z)),`$ (33)
where

$$\gamma =\sqrt{F^2+\frac{1}{4}},k_+=\sqrt{k^24F|P|}.$$

The above results are valid for $`|P|F`$. When $`|P|>F`$, an analytical continuation can be performed. And it was found that there is no discrete part in the spectrum, the wavefunctions and energy are of the same form as that for $`|P|F`$, but now valid for this region of parameters.

## III Dephasing time

In the weak inhomogeneity limit $`\delta L`$ where $`L`$ is the linear size of the sample, $`F1`$ for not very weak field. Then the potential well (given by Eq.(9)) is deep enough to host many discrete levels, which are reminiscence of the Landau levels. Recall that low energy states dominate in the cooperon propagator, the continuous part(with energies higher than the barrier height) gives negligible contribution, thus only the discrete part of the spectrum is important in the cooperon. It can be shown that in this limit, inhomogeneity brings a correction of order $`O(\frac{1}{F})`$ to the usual uniform field weak-localization magnetoresistance. And the dephasing time has the same form with the situation of a uniform field. The opposite limit, $`F1`$, which we will focus on, is more interesting.

Assuming the field described by Eq.(3) can be realized in experiment, let’s show that for realistic parameters it is possible to have the discrete part of the spectrum absent in the strong inhomogeneity limit. For $`\delta 100\mathrm{n}\mathrm{m}`$, $`F\frac{H\delta ^2}{\mathrm{\Phi }_0}10H`$ if $`H`$ is in Tesla, thus up to $`H100\mathrm{G}\mathrm{a}\mathrm{u}\mathrm{s}\mathrm{s}`$, one can take $`F`$ as a small quantity. Then we check if the discrete part of the spectrum exists in this limit. The criterion for its existence is $`P=\frac{p_y\delta }{\mathrm{}}<F^3<10^3`$. Since $`p_y>\frac{\mathrm{}}{L_\varphi }`$, one has $`P=\frac{p_y\delta }{\mathrm{}}>\frac{\delta }{L_\varphi }`$. For a system with $`L_\varphi <10^4\mathrm{nm}`$, this gives $`P=\frac{p_y\delta }{\mathrm{}}>\frac{10^2\mathrm{nm}}{10^4\mathrm{nm}}10^2`$. Thus the inequality for the existence of discrete levels does not hold in this situation. One therefore concludes there is no discrete part in the spectrum. In this limit,

$$E=\frac{\mathrm{}^2}{2m\delta ^2}[p_x^2+(|p_y|F)^2],$$
(34)
and the cooperon propagator is

$`C(r,r^{};\omega )`$ (35)
$`=\left({\displaystyle \frac{\mathrm{}}{\delta }}\right)^2{\displaystyle 𝑑p_x𝑑p_y\frac{\psi _{p_x,p_y}^{}(x,y)\psi _{p_x,p_y}(x^{},y^{})}{iw\tau +\frac{D\tau }{\delta ^2}\left[p_x^2+(|p_y|F)^2\right]}},`$ (36)
where the integral should be done under the constraint

$$\frac{D\tau }{\delta ^2}\left[p_x^2+(|p_y|F)^2\right]1,$$
(37)
or

$$p_x^2+(|p_y|F)^2q^2,q^2=\frac{\delta ^2}{D\tau },$$
(38)
as it is a condition for the perturbation theory.

Before proceeding, we need to discuss the boundary condition. If we use the wavefunctions Eq.(28) and (33) in their present form, then the cooperon propagator constructed naively in Eq.(36) does not satisfy the zero amplitude boundary condition

$$C(\pm L,y)=C(x,\pm L)=0.$$
(39)
The boundary condition in $`y`$ direction can be taken care of in the same way as in the uniform field case, where the plane wave wavefunctions are replaced by their linear combinations, and $`p_y`$ takes only some discrete allowed values. The boundary condition in the $`x`$ direction can be dealt with in the same manner. We note that from Eq.(26), the eigenstates with $`k`$ and $`k`$ are degenerate, thus a linear combination of them is also an eigenstate. The new eigenstate can be constructed as

$$\psi _{p_x,p_y}^{A,B}(x)=\alpha \chi _{p_x,p_y}^{A,B}+\beta \chi _{p_x,p_y}^{A,B},$$
(40)
which are required to satisfy the boundary condition

$$\psi _{p_x,p_y}^{A,B}(L)=\psi _{p_x,p_y}^{A,B}(L)=0.$$
(41)
As in the case of uniform field, only some discrete $`p_x`$ will be allowed. Thus the continuous part of the spectrum becomes again discrete in this boundary condition. The integration over momentum will be replaced by summation over these discrete allowed values. We note that the dispersion relation is not altered in the construction Eq.(40).

Introducing the dimensionless momenta $`q_x=p_x/F`$, $`q_y=p_y/F`$, $`q^{}=q/F`$, $`P^{}=P/F`$, and going to time domain, we find

$`C(r,r^{};t,t^{})`$ (42)
$`={\displaystyle \frac{\mathrm{}^2F^2}{\delta ^2\tau }}{\displaystyle \underset{q_x^2+(|q_y|1)^2q^2}{}}\psi _{q_x,q_y}^{}(x,y)\psi _{q_x,q_y}(x^{},y^{})\times `$ (43)
$`\times \mathrm{exp}\left\{{\displaystyle \frac{DF^2}{\delta ^2}}\left[q_x^2+(|q_y|1)^2\right](t^{}t)\right\}`$ (44)
$$=\frac{\mathrm{}^2F^2}{\delta ^2\tau }\{\underset{q_x^2<4|P^{}|}{}\psi _{q_x,q_y}^A(x,y)\psi _{q_x,q_y}^A(x^{},y^{})+$$

$$+\underset{q_x^24|P^{}|}{}\psi _{q_x,q_y}^B(x,y)\psi _{q_x,q_y}^B(x^{},y^{})\}_{q_x^2+(|q_y|1)^2q^2}$$

$$\times \mathrm{exp}\left\{\frac{DF^2}{\delta ^2}\left[q_x^2+(|q_y|1)^2\right](t^{}t)\right\}.$$

We see that the magnetic field results in a characteristic time scale beyond which the cooperoon propagator or the interference effect is no longer important. Therefore

$$\frac{1}{\tau _\varphi }=\frac{DF^2}{\delta ^2},$$
(45)
which leads to Eq.(4).

Thus the dephasing rate depends quadratically on the field amplitude in the strong inhomogeneity limit. In this calculation, we see that the dependence of $`\frac{1}{\tau _\varphi }`$ on $`H`$ is a result of the nature of corresponding single particle motion. From Eq.(9), when $`F`$ is large, the potential well is deep enough to accommodate many bound states with discrete levels(which are reminiscence of Landau levels). Dominant contribution to the cooperon is from these low lying levels which results in linear dependence. In contrast, if $`F`$ is small, the potential well can be so shallow that discrete level does not appear, and the particle executes overbarrier motion which is nearly free with continuous spectrum. When this part of spectrum is dominant in the cooperon, the dephasing rate is a quadratic function of the field amplitude. The quadratic dependence implies that the electrons are more slowly dephased compared with in the uniform field case, where the dependence is linear. Physically this is because for a given diffusion time, electrons in this regime can visit a larger area than the constrained motion. Since the magnetic field is nonzero only in a limited region, the phase change accumulated during this time is smaller than that of the uniform case. So the electron has to wander for a longer time before it gets dephased. Therefore the constrained motion of single particle leads to linear dephasing rate dependence on field amplitude, and nearly free overbarrier motion to quadratic.

Eq.(4) is obtained in the limit of small $`F`$, which may be realized with tiny $`\delta `$ and small $`H`$. More realistic for experimental observation is situation with moderate $`\delta `$. In this case, one expects a crossover from the linear(uniform limit) to quadratic(inhomogeneity limit) dependence. This point can be illustrated in the following way. Consider the electron which has diffused for a time $`\tau `$. If $`D\tau <\delta ^2`$, the inhomogeneity is not noticed by the electron. If during this time the electron has already been dephased, $`\tau >\tau _\varphi ^u`$, where $`\tau _\varphi ^u`$ is the dephasing time in the uniform field(given by Eq.(2)), then the inhomogeneity is not important at all. When $`\tau `$ is beyond this time scale, the inhomogeneity effect enters. Thus, for a given $`\delta `$, there is a crossover field implied by

$$D\tau _\varphi ^u=\delta ^2,$$
(46)
which yields

$$H_c=\frac{\mathrm{}c}{4e\delta ^2}.$$
(47)
This is equivalent to saying that the inhomogeneity effect enters when $`F1`$.

On the other hand, if we can adjust the width $`\delta `$, then for a given $`H`$, there is a crossover length scale $`\delta _c`$,

$$\delta _c=\sqrt{\frac{\mathrm{}c}{4eH}},$$
(48)
which is the magnetic length for a uniform field.

In passing on, let’s mention that from Eq.(44), one can see that the magnetoconductance is proportional to $`H^2\delta ^2`$ in this small $`F`$ limit. In the large $`F`$ limit, a logarithmic dependence is expected as in the uniform field case. These can also be established by qualitative argument of Khmelnitskii.

## IV Monte Carlo simulation

A Monte Carlo algorithm has been developed to simulate the dephasing process. In order to be consistent with that implied by Eq. (1), the simulations will be semiclassical in nature. Trajectories will be used and the only quantum mechanical effect will be in the phases. In this approach, a particle performs random walk in a square lattice. The value of the perpendicular field is assigned to a dual lattice. We trace all the closed loops that are formed. Once a loop is formed, all the inner points are picked up(for technical details, see ), and we calculate and record the phase change $`\varphi _i`$ due to this loop. Then the trajectory of this loop is erased, and the particle continues the random walk. The phase accumulated this way is

$$\varphi =\underset{i}{}\varphi _i,$$
(49)
whose average is zero since the loop has equal possibility to be clockwise and counter-clockwise. Then

$`\delta \varphi ^2`$ $`=`$ $`\varphi ^2`$ (50)
$`=`$ $`{\displaystyle \underset{i}{}}\varphi _i^2.`$ (51)
The random walk stops as soon as

$$\delta \varphi ^21,$$
(52)
and we specify the time(total number of random walk steps) that has been spent to reach this as the dephasing time. An ensemble of such walks are performed and the dephasing time is averaged. In the simulations, the mean free path $`l`$ and the flux quantum $`\mathrm{\Phi }_0`$ have been set to 1, so that

$$\varphi _i=\underset{j}{}H_j,$$
(53)
where $`j`$ runs over all the inner points in the dual lattice of loop $`i`$.

Several features make this simulation of qualitative nature. First, the criterion Eq.(52) is a qualitative one. Secondly, in the simulation, the impurities are assumed to be on regular lattice sites, with the mean free path as lattice constant. Thirdly, sometimes we choose the field not necessarily the original one, but qualitatively the same. However, essential physics is not lost despite of its qualitative nature, and this simulation is very efficient in uncovering dependence of the dephasing time on the field amplitude, which is particularly appropriate for situations where the field modulation is complicated so that it is difficult to make progress with analytical approaches.

As a check of the algorithm, we have simulated the uniform field case. Excellent linearity in the plot of dephasing rate against $`H`$ is observed, in agreement with Eq.(2). For comparison with the modulated field we studied in the previous section, the simulation is performed with the following field

$$H(x)=\{\begin{array}{cc}H,\hfill & |x|\delta ,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}$$
(54)
which is qualitatively the same as Eq.(1). The results plotted in Fig.(1) show a clear linear dependence for large $`\delta `$, and a crossover from linear to quadratic dependence for moderate $`\delta `$. The boundary is set at $`L=5000l`$ in this simulation. Trajectories which touch the boundary are excluded. The results are obtained by averaging over $`10^4`$ random walks. When the inhomogeneity effect enters, fluctuations become significant in the result. This is because there are extreme trajectories that wander for a long time in the zero field region. However, by recording not only the average but also the standard deviation, we are able to get an estimation $`2.0\pm 0.1`$ for the exponent in the inhomogeneous limit, in agreement with Eq.(4). Thus the numerics support analytical results.

## V Generalization to modulated magnetic field periodic in one direction with zero mean

The qualitative feature of the results, linear dependence in weak inhomogeneity limit and quadratic in strong inhomogeneity limit, may be more general. Consider the situation with magnetic field modulated periodically in one direction with zero mean. While the weak inhomogeneity limit is easy to understand, we shall focus on the opposite limit to discuss possible quadratic dependence of $`1/\tau _\varphi `$ on $`H`$.

Assume the magnetic-field profile is such that nearest neighbor stripes of width $`\delta `$ have the same magnitude but opposite sign. Consider the phase change of a closed loop due to this magnetic field. Since the field is periodic with zero mean, what is important for the magnetic flux is the ratio of linear size of the loop with the width of one stripe, $`\sqrt{D\tau _\varphi }/\delta `$. If this ratio is an even integer, then the net flux is zero. If it is an odd integer, then the flux is equal to that piercing through a single stripe. In general this ratio may fluctuates, however it is clear that the average flux is proportional to $`\sqrt{D\tau _\varphi }H\delta `$. Demanding the phase change due to this flux to be 1, one obtains Eq.(4), with probably an additional numerical coefficient.

For a general magnetic field periodic in one direction, the corresponding Schr$`\ddot{\mathrm{o}}`$dinger equation is often difficult to solve analytically. For example, for $`𝐇=(0,0,H\mathrm{sin}(x/\delta ))`$, one can choose $`𝐀=(0,H\delta \mathrm{cos}(x/\delta ),0)`$. Then the resulting Schr$`\ddot{\mathrm{o}}`$dinger equation is a Whittaker-Hill equation whose solution may be reduced to three-term recurrence relations, but an explicit analytical result for the spectrum is not known. However some insight can be gained by analyzing the structure of energy spectrum. In general, if the magnetic field is $`(0,0,H(x))`$, in the Landau gauge $`\mathrm{A}=(0,A(x),0)`$. Separating variables in the wavefunction, the electron motion in $`x`$ direction is described by a Schr$`\ddot{\mathrm{o}}`$dinger equation with potential

$$V(x,p_y)=\left[p_y\frac{e}{c}A(x)\right]^2,$$
(55)
which is also a periodic function. For given $`p_y`$, this potential is bounded from above as that described by Eq.(9). The profile of this potential is a series of potential wells joining with each other. The energy spectra of electrons in a family of such modulated magnetic fields are calculated by Ibrahim and Peeters. For electron energy less than the barrier height, the wavefunctions in neighboring wells overlap and spread. So the discrete levels in single wells now form minibands due to the periodic structure. When $`H`$ and $`\delta `$ are large, the barrier can be quite high, and states in these minibands dominate in the cooperon propagator. According to the experience in dealing with the previous case, the dephasing rate is expected to be a linear function of $`H`$. When $`H`$ and $`\delta `$ are small, however, the potential well is not deep enough to support these minibands, and dominant contribution is expected from nearly free overbarrier motion, which leads to a quadratic dependence.

To examine these considerations, computer simulations are performed for the following field

$$H(x)=\{\begin{array}{cc}H,\hfill & 2n\delta x<(2n+1)\delta ,\hfill \\ H,\hfill & (2n1)\delta x<2n\delta .\hfill \end{array}$$
(56)
The results are shown in Fig.2. The curves show much resemblance to that in Fig.1, and a crossover from linear to quadratic dependence is evident. There are also some differences which will be noted here. First, for a given $`\delta `$, the crossover field $`H_c`$ is smaller in the field (56) than in (54), which suggests when $`H`$ is small, electron in the field (56) is more quickly dephased than in the field (54). This may look strange at first sight, since the field (56) has zero mean, one may expect the electrons in this field are more slowly dephased. However, this result can be understood by noticing what enters is the variance of the phase instead of the average, and the field (56) is nonzero everywhere in the plane. Another difference is in the field (56) there is broader crossover region. This is associated with the fact that, in the periodic situation, states with energy lower than the barrier height are not really bound states. They share some feature of the plane waves according to Bloch’s theorem. As discussed previously, this feature has the tendency of leading to a quadratic dependence. When $`\delta `$ is moderate, the barrier is not very high for many values of $`p_y`$. Then these states can be dominant in the cooperon, resulting in a broad crossover region.

Numerical experiments are also performed for situations with some other magnetic field modulations, for example, that with circular symmetry, or periodic in both directions. Similar crossover behaviors have been observed. These simulations as well as the qualitative understanding of the origin of the dependences, suggest that the crossover from a linear to quadratic dependence in the dephasing rate on field amplitude can be quite general. For situations where the corresponding single particle potential is bounded from above, the particle motion can exhibit both bound and free nature, leading to a crossover from linear to quadratic dependence.

## VI Conclusion

We have studied the dephasing time of disordered two-dimensional electron gas in modulated magnetic field $`𝐇=(0,0,H/\mathrm{cosh}^2((xx_0)/\delta ))`$. It is shown that in the weak inhomogeneity limit, $`\delta >>L`$, where $`L`$ is the linear size of the sample, $`\tau _\varphi ^1`$ is proportional to $`H`$. This happens when the bound states with discrete spectrum of the corresponding Schr$`\ddot{\mathrm{o}}`$dinger equation dominate in the cooperon propagator. While in the strong inhomogeneity limit, the dependence is quadratic $`\tau _\varphi ^1=D\left(\frac{e}{\mathrm{}c}\right)^2H^2\delta ^2`$. In this case, the nearly free overbarrier motion gives dominant contribution to cooperon. In the intermediate regime, a crossover between these two limiting situations occurs at $`H_c=\frac{\mathrm{}c}{4e}\delta ^2`$. A semiclassical Monte Carlo algorithm has been developed to study the dephasing time, which is of qualitative nature but efficient in uncovering the dependence of $`\tau _\varphi `$ on $`H`$ for arbitrarily complicated magnetic field modulations. Computer simulations support analytical results. These considerations are generalized to the situation with magnetic field modulated periodically in one direction with zero mean, where a similar crossover is observed. We believe this crossover between linear and quadratic dependence can be expected for a large class of modulated magnetic fields where the corresponding single particle potential is bounded from above, so that the motion exhibits both bound and free nature depending on parameters.

Acknowledgements

I am grateful to Vladimir Kravtsov for helpful discussions. Discussions on the simulation algorithm with Dr. Z. F. Huang, Dr. F. Ricci-Tersenghi, and a communication with Prof. K. Kremer and Prof. K. Binder are acknowledged.


format: text

text:
# Are the Narrow Line Regions in Active Galaxies Dusty and Radiation Pressure Dominated?

## 1 Introduction

Optical spectroscopy of the ionized gas around galactic nuclei provides strong constraints on the excitation mechanisms. Line ratio diagrams of active galactic nuclei (AGN), such as those of Veilleux & Osterbrock (1987), allow us to distinguish three major classes; nuclei excited by starbursts, and the two types excited by a bona-fide active nucleus: the Seyfert narrow line regions (NLRs) and extended narrow line regions (ENLRs), and the Low Excitation Nuclear Emission Line Regions (LINERs). As far as the Seyfert galaxies are concerned, the compilations of Veilleux & Osterbrock (1987), Veilleux(1991a,b,c), and Véron-Cetty & Véron (2000) show that the emission line ratios are remarkably uniform within a given object, and very similar (with typically less than 0.5 dex variation) from one object to another.

The standard paradigm proposes that the NLR and the ENLR are excited by photons originating at or near a compact nuclear source (see, e.g. Osterbrock (1989)) having a smooth featureless power-law, or broken power-law EUV ionizing spectrum. Within this model, the clustering of the observed line ratios within such a restricted domain of parameter space presents a problem, since in such models the excitation is in large measure controlled by the ionization parameter, $`U`$, and this is free to vary over a wide range. The implications of this is that there should exist photoionized regions with very high excitation corresponding to high $`U`$, but objects such as these are missing in the Veilleux & Osterbrock (1987) diagnostic diagrams. Instead, the constraints derived from the line ratio diagnostic diagrams imply that $`U`$ has a similar value in all NLRs. Modelers have been therefore forced to make the arbitrary (and possibly unphysical) assumption that the gas density in the ionized clouds must fall exactly as the inverse square of the distance from the nucleus.

One rather successful variation on the photoionization models is to assume that the power-law EUV spectrum passes through a highly-ionized screen of diffuse gas, before being absorbed in the vicinity of dense (much higher pressure) clouds (Binette et al. (1996), Binette et al. (1997)). Although such models produce generally good results when compared with observed ratios, these have the difficulty that the geometry of the ionized gas is ill-defined and necessarily becomes a new free parameter. These models can also be questioned on physical grounds, since the assumed pressure difference between the cloud and inter-cloud components is not explained.

Another variation is to consider the integrated emission from a whole family of photoionized clouds of different densities (Komossa & Schulz (1997); Ferguson et al. (1997)). These models can be made to fit the observed spectra quite well, but the geometry of the clouds and their relative contributions to the emission remain free parameters and are poorly constrained.

In any of these photoionized models, the (often quite violent) internal dynamics of the ENLR are not addressed. These include clear evidence of non-gravitational motion, and evidence for outflow at velocities up to (or even in excess of 1000 km s-1) (eg Pedlar et al. (1989), Allen et al. (1999)). It was results such as this that led Wilson & Willis (1980) to suggest that the nucleus ejects radio components that interact with ambient gas and replenish the high kinetic energy and ionization of the NLR. This led to the development of photoionization by fast radiative shocks, which generate a strong internal photoionizing radiation field (Dopita & Sutherland (1995), (1996)). This model is successful in reproducing the strong lines observed in many AGN (see review by Dopita (2000)). In addition, it can explain other features which simple photoionization models fail to address. For example, where the NLRs are spatially resolved, there often are found to be strong correlations between radio power and either line luminosity de Bruyn & Wilson (1987) or line width Wilson & Willis (1980). Such correlations exist not only for Seyfert galaxies, but persist up to much more luminous classes of radio galaxies. However, the fast shock model has its limitations. It has difficulty reproducing the very strong coronal lines of \[Fe vii\], \[Fe x\], \[S viii\], \[S xii\], \[Si vi\] or \[Si vii\] (Wilson & Raymond (1999)) which imply photoionization with very high local ionization parameters. These models also cannot reproduce the large outflow velocities inferred in the \[O iii\] line in some objects such as NGC 1068 (e.g. Kraemer & Crenshaw (2000); Cecil et al. (2001)). In shock models, most of the \[O iii\] emission is expected to be generated in the photoionized precursor of the shock, which theoretically should be the undisturbed gas in the galaxy. Furthermore, the NLR lines of (probably) the majority of Seyfert galaxies are quite narrow and consistent with line broadening due to galactic rotation rather than being produced in shocked outflows.

In this paper we argue that some of these properties can be explained in terms of radiation pressure dominated dusty photo-ablating clouds. In the past, dust has been considered in photoionization models both through its contribution via photoelectric heating (e.g. Bottorff et al. (1998)) and as shielding component thanks to its large opacity in the FUV. However, the effects of dust photoelectric heating and of the radiation pressure on dust on both the ionization and hydrostatic structure of the NLR clouds has not been examined in a systematic fashion before. In the next section we detail these effects and show how these provide a natural explanation as to why the inferred NLR ionization parameter seems to be confined to the range $`2.5<\mathrm{log}U<2`$. In §3 we describe the basics of our dusty, radiation pressure dominated models. These are compared with other types of photoionization models, allowing us to separately understand the effects of photoelectric heating and of radiation pressure on the emergent spectrum. The predictions of all these models are compared with observations of NLR of Seyferts on some key emission line diagnostic diagrams. The concluding remarks are given in §4. In following papers we will investigate further the effects of different parameters on our radiation pressure dominated photoionization models and present a full hydrodynamical model exploring such structures.

## 2 A New Paradigm: Dusty, Radiation-pressure Dominated Photoionization

The effects of radiation pressure from line absorption have previously been considered in photoionization models of the emission line regions of AGN (Elitzur & Ferland (1986)). This is important in determining the internal density structure at large ionization parameter. The effect of dust within the ENLR has also been considered previously (Netzer & Laor (1993)). However, the radiation pressure exerted by the ionizing source acting upon dust within the photoionized gas has usually not been accounted for by the standard photoionization models, with a few notable exceptions (e.g. Baldwin et al. (1991)). The standard photoionization modelling assumes an isochoric (or constant density) structure within the ionized slab instead of an isobaric (constant pressure) structure. Dust is important not only because of the radiation pressure acting upon it, but also because of the photoelectric heating which it produces in the photoionized plasma, and in the way it competes with the gas in absorbing the ionizing photons.

However, all of this critically depends on whether dust can survive the either the strong EUV radiation field, or the relatively energetic conditions characterizing the photoionized region. In the ionized plasma of the ENLR, grains are electrically charged, locking the motion of the dust to that of the ionized gas. Both the ionized gas and the dust is likely to be streaming at roughly sonic velocities from ionization fronts located around photo-evaporating dense clouds located in the “ionization cones” of the AGN. The problem of dust survival is therefore reduced to the question of whether it can survive for a time equal to its passage through the high-emissivity region near the cloud surface as discussed in §3. Our detailed photoionization models show that the average gas temperature is $`10^4`$K, at which the thermal sputtering of grains is negligible and incapable of destroying the grains in so short a timescale. This is in contrast with the case of radiative shocks in which gas temperatures reach $`10^510^6`$K and grain destruction by sputtering occurs on very short timescales (Tielens, McKee, Seab, & Hollenbach (1994)). Similarly, our models show that grain temperatures are too low for thermal evaporation of grains by quantum heating provided that the dust grains are large enough ($`3`$ nm ). A third potential dust destruction mechanism is that of Coulomb explosion, in which the grain charging becomes large enough that the Coulomb repulsion between the atoms composing the grain exceeds the atomic bond strength. According to our models this condition is met for grains smaller than about 10 nm when the ionization parameter is of order unity. In the models described here we have simply set the minimum grain radius equal to 10 nm to allow for the destruction of the smaller grains. We will return to the detailed physics of grain destruction processes in these environments in a future paper.

In NLR regions observed by Oliva et al. (1994), Moorwood et al. (1996), and Kraemer & Crenshaw (2000) the coronal lines appear remarkably strong as compared to H$`\beta `$. The existence of such coronal gas requires that both the ionization parameter and the radiation pressure are very large (Binette et al. (1997)). Elements such as Fe and Si are normally largely sequestered in the grains. The strength of the coronal lines of these elements requires that in the coronal gas, at least, the dust has been largely destroyed. Our model, as depicted in Figure 1 is therefore one of a dusty, radiation-pressure dominated region surrounding a photoevaporating molecular cloud, which in turn is surrounded by a coronal halo within which the dust has been largely destroyed. In this paper, we will consider only the emission spectrum generated by the dusty dense inner region near the stagnation point in the flow around the cloud. This will be the highest emissivity region in the flow, and it will be the region in which the radiation pressure gradient is matched by the gas pressure gradient, allowing us to construct static rather than dynamic photoionization models.

Consider $`\mathrm{\Xi }_0`$, the ratio of the available pressure of ionizing radiation, $`\mathrm{\Phi }`$, to the gas pressure at the irradiated photoionized surface of a photoevaporating cloud111This definition implies a factor 2.3 lower value than that of Krolik et al. (1981) :

$$\mathrm{\Xi }_0=\frac{P_{\mathrm{rad}}}{P_{\mathrm{gas}}^0}\frac{\mathrm{\Phi }_0/c}{2.3nkT_e}U_0\frac{\epsilon _{\mathrm{EUV}}}{2.3kT_e}$$
(1)
where $`k`$ is Boltzmann’s constant, $`c`$ the speed of light, $`n`$ and $`T_e`$ are the hydrogen density and the electron temperature at the irradiated cloud surface (indicated by a subscript zero elsewhere). We define $`S_{}=\mathrm{\Phi }_0/\epsilon _{\mathrm{EUV}}=\frac{\phi _0(\nu )}{h\nu }𝑑\nu `$ to be the ionizing photon flux (cm-2s-1) impinging at the surface such that $`U_0=S_{}/nc`$ is the canonical ionization parameter. $`\epsilon _{\mathrm{EUV}}`$ ($`=h\nu _{\mathrm{EUV}}`$) is the local average EUV photon energy.

In hydrostatic equilibrium, the gas pressure gradient must match the local radiative volume force exerted by photon absorption, and this force can be separated into its two components, the force acting on the gas through photoionization, and the component due to dust absorption;

$`F_{\mathrm{rad}}(x)=`$ $`{\displaystyle \underset{m=1}{}}{\displaystyle \underset{i=0}{}}n_i(X_m^{+i}){\displaystyle _{\nu _m^i}^{\mathrm{}}}{\displaystyle \frac{\phi (\nu )}{c}}a_\nu (X_m^{+i})𝑑\nu `$
$`+n_H{\displaystyle _{\nu _m^k}^{\mathrm{}}}\left[\kappa _{\mathrm{abs}}\left(\nu \right)+\kappa _{\mathrm{sca}}\left(\nu \right)\left(1g\left(\nu \right)\right)\right]{\displaystyle \frac{\phi (\nu )}{c}}𝑑\nu `$
where $`\phi (\nu )`$ is the local radiation flux at frequency $`\nu `$, $`n_i(X_m^{+i})`$ the number density of ion $`+i`$ of atomic species $`m`$, $`a_\nu (X_m^{+i})`$ the corresponding photoionization cross section with threshold $`\nu _m^i`$. In the dust absorption term, $`\kappa _{\mathrm{abs}}\left(\nu \right)`$ and $`\kappa _{\mathrm{sca}}\left(\nu \right)`$ is the dust absorption cross section and scattering cross section respectively, normalized to the hydrogen density and $`g\left(\nu \right)=\mathrm{cos}\theta _\nu `$ is the mean scattering angle of the dust grains at this frequency. The boundary condition is given by $`P_0`$, the pressure at the irradiated face.

We can readily estimate a critical value of $`U_0`$ above which dust absorption dominates over the photoelectric absorption. For an equilibrium photoionized plasma, the absorption of ionizing photons in a slab is simply equal to the local recombination rate in the plasma:

$$dS_{}/dx=\alpha _Bn^2$$

where $`\alpha _B`$ is the recombination rate to excited states of hydrogen. The absorption of ionizing photons by dust is given by:

$$dS_{}/dx=\kappa nS_{}$$

where $`\kappa `$ is the effective dust opacity of order $`10^{21}\mathrm{cm}^2`$ per H atom. Dust absorption becomes relatively more important as the strength of the ionizing field increases, and dust becomes the dominant opacity in the plasma when:

$$U_0>\frac{\alpha _B}{c\kappa }$$
(2)
For a typical value of $`\alpha _B=2\times 10^{13}\mathrm{cm}^3\mathrm{s}^1`$, the critical ionization parameter is $`0.007`$.

This value is closely related to that at which radiation pressure starts to dominate either the gas pressure or the dynamical acceleration of the plasma. If $`\epsilon _{\mathrm{EUV}}`$ is the mean energy of the ionizing photons absorbed by gas, and $`\epsilon _{\mathrm{FUV}}`$ is the mean energy of the ionizing absorbed by dust, equation (2) can be approximated by:

$`F_{\mathrm{rad}}(x)`$ $`=`$ $`{\displaystyle \frac{1}{c}}\alpha _Bn^2\epsilon _{\mathrm{EUV}}+{\displaystyle \frac{1}{c}}S_{}\epsilon _{\mathrm{FUV}}\kappa n`$ (3)
$`=`$ $`{\displaystyle \frac{n^2\epsilon _{\mathrm{EUV}}}{c}}\left[\alpha _B+{\displaystyle \frac{\xi }{\psi }}U_0\kappa c\right]`$
where $`\epsilon _{\mathrm{EUV}}/\epsilon _{\mathrm{FUV}}=\psi `$ , and the effective ionization parameter for the dust–absorbed photons $`U_{\mathrm{FUV}}=\xi U`$. Generally speaking, both $`\xi `$ and $`\psi `$ will be greater than unity. Clearly, dust becomes the dominant term in this equation when:

$$U_0>\frac{\psi \alpha _B}{\xi c\kappa }$$
(4)
this expression is the same as equation (2) except for the ratio $`\psi /\xi `$ , which is of order unity. In the typical ionization field expected around an AGN, $`\phi _0(\nu )\nu ^{1.4}`$and with a ‘standard’ dust model with solar abundances, a Mathis, Rumpl & Nordsieck (1977) (MRN) size distribution $`n(a)a^{3.5}`$ and a lower size cutoff of 10 nm , we estimate that this condition gives a critical ionization parameter $`U_00.01.`$

If the ionization parameter into a static plane-parallel slab of ionized gas is increased above this critical value, then the radiation pressure gradient imposed as the EUV radiation field is absorbed through the ionized slab will induce a corresponding local gas pressure gradient. Thus, if the initial (outer layer) pressure of the ionized gas is $`P_0`$ and the radiation pressure in the EUV field is $`P_{rad}`$ then the gas pressure close to the ionization front is $`P_0+P_{rad}.`$ Thus, at high ionization parameter, the pressure in the ionized gas close to the ionization front is determined by the externally imposed ionization parameter, $`U_0,`$ and the local ionization parameter becomes independent of the initial ionization parameter. In this way, the emission line spectrum in the low- and intermediate- ionization species will become effectively independent of the external ionization parameter. This is the core of our explanation of why the emission line spectra of NLR are so similar one to another.

While $`U_0`$ defines the ionization parameter at the face of the cloud, one can define a mean ionization parameter $`\overline{U}=S_{}/\overline{n}c=U_0n_0/\overline{n}`$ where $`\overline{n}`$ is simply the mean density, that is an emissivity averaged density using a weight $`n_e^2`$. Binette et al. (1997) showed that $`\overline{U}`$ tends asymptotically towards a constant value in the limit of high $`U_0`$.

## 3 Results of Models

We have constructed static, plane-parallel, radiation-pressure dominated, dusty photoionization models using the code MAPPINGS III, which includes the implementation of eqn. (2) to evaluate the local radiation pressure gradient and hence the local density gradient. MAPPINGS III treats several of the properties of dust such as photoelectric heating, as described previously by Dopita & Sutherland (2000). The dust model includes both silicaceous and carbonaceous grains with a MRN grain size distribution and uses the absorption and scattering data of Laor & Draine (1993) (using the more recent smoothed UV data for silicates)222Available at http://www.astro.princeton.edu/$``$draine/dust/dust.diel.html. Here we have chosen the simplest “standard” values for the input parameters. For example, the models all used the Anders & Grevasse (1989) solar abundance set, with the depletion factors for each element given in Table 1. The input ionizing spectrum for the following models was a standard powerlaw continuum with $`\phi _0(\nu )\nu ^{1.4}`$.The effect of varying these and other input parameters will be investigated in a following paper (Groves et al. (2001) in preparation).

To facilitate comparison, and to separate the effects of internal density structure, radiation pressure and photoelectric heating by the grains, we have run four families of models with increasing ionization parameter $`U_0`$ (defined at the outer boundary of the ionized slab). In order of increasing physical sophistication (and reality) these are:

1. Isochoric (constant density) dust-free models with heavy element abundances in the gas phase which match those of the models including dust.

2. Isobaric dust-free models. Again these have gas-phase heavy element abundances matching the models with dust. These models are isobaric in the sense that the product $`nT`$ is kept constant, but the radiation pressure is not explicitly included.

3. Isochoric dusty models.

4. Isobaric dusty models without radiation pressure. These provide a dusty comparison to the complete model which shows the difference that radiation pressure makes.

5. Complete isobaric dusty models, which include the effect of radiation pressure on the density structure of the ionized gas.

To avoid differential effects due to different electron density in the low-excitation cases, in all cases the initial density was adjusted so that the density of the region emitting in the \[S ii\] lines was maintained constant at about 550 cm-3 to within 10%. This density was chosen to be typical of the electron densities inferred from the \[S ii\] $`\lambda \lambda 6717/6731`$ ratio in the NLRs of Seyfert Galaxies. The models were terminated at a point where the gas temperature dropped rapidly to below 5000K and the emission in the optical forbidden lines is quenched.

Generally in the literature we encounter only the first or the third of these sets of models. The profound difference in the internal structure of the isobaric dusty models including radiation pressure compared with a “standard” isochoric model is emphasized in Fig. 2 which shows the internal density and temperature structure produced in the dusty, radiation pressure dominated slab model with ionization parameter of $`U_0=0.5.`$ In this model, the radiation pressure in the EUV radiation field is roughly 30 times the initial gas pressure. At the “front” surface of the cloud the opacity of the cloud is entirely dominated by dust, and the photoelectric heating by dust is very large. As a consequence, the electron temperature in this zone is very high, of order 25000K. The electron temperature decreases strongly throughout the model towards the ionization front, and this steepens the density gradient in the model.

This strong density gradient which is set up between the high ionization and low ionization regions can provide a physical explanation of the pressure difference between these regions which was required in the $`A/M_\mathrm{I}`$ models of Binette et al. (1996). As expected, the pressure in the low-ionization region scales as the radiation pressure at high $`U_0.`$ Thus, the local ionization parameter characterizing the low and intermediate ionization region of the model becomes independent of $`U_0.`$

The emission line ratios produced by the various classes of model can be quite different one from another. We have chosen to display here only those which provide a good separation between the different classes of model, and so have a diagnostic utility. Any of the of line emission ratios used by Veilleux & Osterbrock (1987) will give such a separation.

In figure 3 we plot one of these diagnostic diagrams, the \[S ii\] $`\lambda \lambda 6717+6731`$Å/H$`\alpha `$ ratio $`vs.`$ the \[O i ii\] $`\lambda 5007`$Å/H$`\beta `$ ratio. The observational points are from Allen et al. (1999). For the dust-free case, both the isochoric nor the isobaric sequence can reproduce the observed data, though the fit would be improved with a flatter spectral index. The closest fit is found in the region $`\mathrm{log}U_02.8`$. As has been shown by many other authors, only a very restricted range of ionization parameter provides a reasonable fit.

For the dusty models, both the isochoric and isobaric sequences approach the region of the observed NLR at high ionization parameter. The dusty, radiation pressure affected isobaric model shows very slow variation in the line ratios with $`U`$ at high ionization parameter and these line ratios are in the region of interest. The dusty, isochoric models actually begin to head away from the observations as we increase in ionization parameter.

Figure 4 is a close up of the region of interest, showing only isobaric models. In addition to the two isobaric models described above, a dusty, isobaric model without radiation pressure is included for comparison. This figure illustrates the self-regulating effect of the radiation pressure in producing nearly constant line ratios for $`U_0>10^2`$. While the radiation pressure free case does reproduce the observed data quite well, it does so only for a small range of ionization parameter, $`2.3<\mathrm{log}U_0<1.3`$. By $`\mathrm{log}U_01.0`$ the model has left the region in which the observations lie. The complete dusty models, which include the effects of radiation pressure, not only lie in the region of the observed data but also display line ratios which stagnate with respect to changes in $`U`$, over a broad range of values. The implication of this is that, provided $`U_010^2`$ , a cloud of given density could be moved radially by a large amount without producing any appreciable change in its excitation state. Such a behavior is frequently seen, and is interpreted in the “classical” photoionization models as indicating a constant ionization parameter.

The conclusions derived from figure 3 and 4 are reinforced when we replace the \[O i ii\] line by the higher-excitation \[Ne v\] line. This is shown in figure 5. Now the range of acceptable models is appreciably reduced. From inspection of both Figures 3 and 5 we see that the dusty models that include radiation pressure provide acceptable solutions in the range $`1.6<\mathrm{log}U_0<0.`$ The acceptable range of solutions for the dusty isochoric models is much more restricted; $`2.0<\mathrm{log}U_0<1.2,`$ and there are no really acceptable solutions for any dust-free model. This illustrates a major advantage of the radiation pressure dominated model in that it encompasses a wider range of ionization species than any simple isochoric or isobaric model. The reason for this is that as density increases inward, one is effectively spanning a very wide range in the local ionization parameter. Furthermore, the extra heating provided by dust photoheating scales with $`U_0`$ and the higher electron temperatures enhance the intensities of the high ionization species in a way that the dust free models cannot account for. This effect is further enhanced by dust absorption which strongly reduces the luminosity of the hydrogen recombination lines in the high ionization zone. At the same time, it enhances the He ii $`\lambda 4686`$Å/H$`\beta `$ ratio to values as high as 0.5-0.6 which is in better agreement with observation, and similar to the value obtained in the matter-bounded models of Binette et al. (1996).

If we look only at line ratios which are sensitive to the excitation, then dust-free models might seem to provide an acceptable fit. An example is shown in figure 6 which plots the \[O i ii\]/\[O i\] line ratio against the \[Ne v\]/\[Ne iii\] ratio. Here similar values of the ionization parameter produce very similar values of the \[Ne v\]/\[Ne iii\] ratio, while the \[O i ii\]/\[O i\] line ratio is more strongly affected by the temperature in the \[O i\] - emitting zone. Note that, again the dusty models including radiation pressure provide acceptable solutions in the range $`1.6<\mathrm{log}U_0<0.`$

The effect of photoelectric heating is best investigated in terms of the electron temperature in the \[O i ii\] - emitting region. This is shown in figure 7 which plots the temperature sensitive \[O i ii\] $`\lambda 4363`$Å/5007Åratio against \[O i ii\] $`\lambda 5007`$Å/H$`\beta .`$ For clarity we have omitted the isochoric models, since these largely overlap with the isobaric cases. The observed points are drawn from both Allen et al. (1999) and Tadhunter, Robinson & Morganti (1989). The dust-free models are clearly much cooler than the dusty models and fail to reproduce the observed \[O i ii\] temperatures for any ionization parameter. This is the famous “temperature problem” of NLR, and which was discussed in detail by Tadhunter, Robinson & Morganti (1989). Photoelectric heating by dust, proposed as the solution to a similar problem occurring in planetary nebulae by Dopita & Sutherland (2000) seems to provide enough additional heating in NLR to bring the predicted electron temperatures into the observed range. Once again, for dusty isobaric radiation-pressure dominated models, acceptable models for the \[O i ii\] $`\lambda 4363`$Å/5007Å ratio lie in the range $`1.6<\mathrm{log}U_0<0.`$

A final major effect of the addition of dust in these models is to limit the H$`\beta `$ surface brightness as a consequence of the absorption of ionizing photons by dust. Since dust absorption strongly dominates over the gas absorption in the regime of high ionization parameter, this limits the ionized column (Str$`\ddot{\text{o}}`$mgren column) to a maximum value. For a fixed density or pressure in the region emitting the \[S II\] lines, the surface brightness of the ionized zone in H$`\beta `$ saturates at a maximal value at high $`U_0`$. By contrast, at low $`U_0,`$ dust absorption is negligible compared with gas absorption, and the H$`\beta `$ flux scales with the ionization parameter. This behavior is shown in figure 8.

If the \[S II\] density is lowered, then, at at a given $`U_0`$ the H$`\beta `$ flux simply scales with the density in the low-$`U_0`$ regime. However, since the dust and the gas absorption terms have a different density dependence, at a given $`U_0`$ in the high-$`U_0`$ regime, a lowering of the density in the \[S II\] zone leads to a greater fraction of the incident radiation field being absorbed by dust. The result of this is that the saturation value of the H$`\beta `$ surface brightness is reduced by a greater factor than simply the ratio of the densities. The roll-off between the low-$`U_0`$ and high-$`U_0`$, regimes must therefore occur at a slightly lower value of $`U_0`$ than in the high density case.

The plane-parallel models presented here strictly apply only along the (stagnation) line which lies directly towards the ionizing source on the front face of a photoevaporating cloud. Towards the sides of the cloud, lateral motions will be developed in the ionized plasma by the transverse gas pressure gradient. Around the cloud edges, the radiation pressure force acts radially outwards from the central engine. According to equation (2), this drives an accelerating flow radially away from the central AGN. Such a radiatively-driven outflow was considered by Binette (1998) in the context of the coronal gas. In the case of NGC 1068, new HST STIS data presented by Cecil et al. (2001) shows \[O iii\] tails of clouds about 100 pc out from the nucleus being accelerated up to 2000 km s-1 or greater, over a distance of about 0.15 arc sec (10 pc). This corresponds to an acceleration of order $`6\times 10^4`$cm s-2. Let us now examine whether this acceleration could be due to radiation pressure.

According to the observations of Kraemer & Crenshaw (2000), the density given by the \[S ii\] lines in NGC 1068 is about 2500 cm$`^3.`$ In the high ionization parameter models, this density is determined entirely by the radiative pressure, and such a pressure is obtained from a central source giving a surface luminosity of $`1500`$ erg cm2 s-1 at the outer boundary of the photoablating cloud. This figure is in excellent agreement with that inferred from the bolometric FIR luminosity (Dopita et al. (1998)) $`1.5\times 10^{11}`$ $`L_{}`$. This gives a flux at a distance of 70 pc from the nucleus (typical for these clouds) of about 1700 erg cm2 s$`^1.`$

If these radiative forces go into accelerating the gas in the ionized tails of the NGC1068 clouds, then from equation (3) we would require a flux of order $`1.5\times 10^4`$ erg cm-2 s-1 to provide the observed acceleration. The inferred flux is a factor ten lower than this, and this would produce a final outflow velocity of order 700km s-1. The only way around this problem would be to increase the force per unit mass of gas by increasing opacity of the dust. This could be done through grain shattering, which according to Jones, Tielens, & Hollenbach (1996) will occur whenever grain-grain collision velocities are greater than the sound speed in the grain material, about 2 km s-1. This could be achieved with either subsonic turbulence in the gas, or else subsonic radiatively-driven grain drift velocities.

A potential mechanism for grain shattering is as follows. Close to the ionization front, and within it, the net grain charge is negative, being dominated by collisions. When the grains are charged their motion is strongly coupled to the gas via Coulomb interactions. As the grains are advected with the gas flow from the ionization front, the local radiation field becomes stronger and the charge on the grains becomes smaller as photoelectric charging becomes more important. Eventually, the net charge on an individual grain becomes zero as the collisional and photoelectric charging exactly balance. At this point the grain can decouple from the gas flow and radiative pressure forces push it back towards the ionization front. The terminal drift velocity that the grain can reach is determined by the gas density, $`\rho `$, and the radiation field intensity, $`F`$; $`\text{v}_{\text{drift}}=(F/c\rho )^{1/2}`$. For the parameters computed for NGC 1068 this could amount to $`200`$ $`\text{km}\text{s}^1`$. However, long before this velocity is reached, the grain will become charged once more, and the grain motion is again locked to the gas motion. Thus, when their charge becomes zero, the grains will tend to oscillate about the position where this happens until another (charged) grain strikes and shatters them. The shattered fragments have their equilibrium point further out in the ionized flow, and of course they are much more vulnerable to the other grain destruction mechanisms we have mentioned. In our models for NGC 1068 grains are completely shattered in the region which produces \[Ne v\] and \[Fe vii\].

The idea that the ionized tails of the cloudlets in NGC 1068 are driven by radiation pressure acting on the shattered dust seems a reasonable one, but it is clear that the full solution of this problem will require a hydrodynamic model which explicitly includes the grain physics and destruction processes. We intend to present such a model in a future paper.

If dust is indeed destroyed at some point in the outflow, the gas will continue to flow outward from the cloud, but radiation pressure effects become relatively ineffective in determining the dynamics. In this case the cloud will become surrounded by a highly ionized coronal region with an equilibrium electron temperature of order $`8\times 10^4`$K. On the basis of continuity of pressure, we estimate that in the case of NGC1068, the density of this region is $`1020`$ cm-3 and therefore that the effective ionization parameter is of order unity, or somewhat greater. This would be sufficient to explain the extreme degree of ionization observed in this coronal region (Kraemer & Crenshaw (2000)).

Finally, let us examine if the hypothesis that dust can be accelerated in the cometary tail of the cloud is also consistent with the hypothesis that dust can also be destroyed in entering the coronal region. From the acceleration in the tail ($`6\times 10^4`$cm s-2) and the length of the tail (10pc) we can estimate the transit time of material through the tail; $`10^4`$yr. Along the stagnation line in the direction of the source, material flows off the cloud at about the sound speed, 20 km s-1. However, the photoionization model for NGC 1068 shows that the thickness of the ionized zone is about $`10^{18}`$cm. Therefore, the dusty plasma flows through this region in about $`1.6\times 10^4`$yr. These timescale are quite comparable, which strongly suggests that the length of the \[O iii\] tails and the production of the coronal region in NGC 1068 are both moderated by the dust destruction, occurring on a timescale of order $`10^4`$yr.

## 4 Conclusions

In this paper we propose that the NLR and ENLR of Seyfert galaxies can be self-consistently modelled by radiation-pressure dominated dusty photoionized regions surrounding photoablating dense clouds. By using the code MAPPINGS III, we have shown that these display a strong density gradient increasing towards the ionization front, are characterized by strong photoelectric heating by the dust and in which the dominant absorption of the ionizing continuum is due to the dust. In these models, the gas pressure close to the ionization front is determined by the radiation field. The models provide:

* A line spectrum consistent with the observations, and which remains approximately invariant against changes in $`U_0`$ that span a range larger than 2 orders of magnitude.
* A higher electron temperature than can be achieved in standard photoionization models, and which could help explain the long standing “temperature problem” of NLR and ENLR.
* An H$`\beta `$ surface brightness which reaches a maximum value at high values of the initial ionization parameter $`U_0`$.

In addition, they may also offer an explanation for:

* The co-existence of the low-ionization zones with a extremely high ionization parameter coronal region.
* Strong radial radiative acceleration of dusty gas which may explain the broad and blue-shifted \[O iii\] profiles seen in many Seyfert galaxies.

We will systematically explore the properties of such models in future papers.

The work of LB was supported by the CONACyT grant 32139-E and was carried out while a Visitor at the RSAA supported by a RSAA visitors grant. M. Dopita acknowledges the support of the Australian National University and the Australian Research Council (ARC) through his ARC Australian Federation Fellowship, and also under ARC Discovery project DP0208445. The authors would like to thank Ernesto Oliva for his many helpful comments and for pointing out an error that the authors missed.