text:
# Wavefronts may move upstream in semiconductor superlattices

## I Introduction

Current instabilities in doped semiconductor superlattices (SL) have been an active subject of research during this decade. For strongly coupled SL, Bloch oscillations  and Wannier-Stark hopping  produce negative differential conductivity (NDC) at high electric fields. This may result in self-sustained oscillations of the current due to recycling of charge dipole domains as in the Gunn effect of bulk n-GaAs . For weakly coupled SL, sequential tunneling is the main mechanism of vertical transport. Under dc voltage bias conditions, stationary electric field domains may form if doping is large enough . Below a critical doping value, the existing charge inside the SL may not be able to pin domain walls, and current self-oscillations appear . These oscillations may be due to recycling of charge monopoles (domain walls) or dipoles depending on the boundary condition at the injecting contact region (in a typical n+-n-n+ configuration with the SL imbedded between highly doped regions, the doping at the emitter region is crucial) . Driven chaotic oscillations have also been predicted  and observed in experiments . Lastly, there are ways to tune the charge inside the SL (and therefore obtain stationary domains or self-oscillations) without replacing it by a different one. For example, by applying a transverse magnetic field  or by photoexciting the SL .

Transport in weakly coupled SL can be described by simple rate equation models for electron densities and average fields in the wells, . Many of the effects related above have been explained by means of a simple discrete drift model . In this model, the tunneling current between two adjacent wells, $`J_{ii+1}`$, equals the 2D electron charge density at well $`i`$ times a drift velocity, which depends on the electric field at the same well. By starting from a microscopic sequential tunneling model, it has been shown that the discrete drift model is a good approximation at low temperatures and for fields above the first plateau of the SL current-voltage characteristic . For low dc voltages on the first plateau, a discrete diffusion (which is a nonlinear function of the field) should be added. This term contains the contribution to $`J_{ii+1}`$ of the tunneling from well $`i+1`$ back to well $`i`$ (which vanishes for large enough electric fields) . In this paper we report an interesting consequence of electron diffusivity at low fields: if the current is sufficiently high, and so is the doping, a domain wall (monopole wave) which connects two domains may travel in a direction opposite to the flow direction for electrons (i.e., upstream, in the positive current direction!). This striking phenomenon is contrary to the usual situation: a monopole either moves downstream (in the direction of the flow of electrons), or it remains stationary, . We substantiate our claim both by numerical simulations of the discrete drift-diffusion model and by rigorous mathematical analysis based upon a comparison principle . Mathematical analysis yields useful bounds for critical values of current and well doping, and for monopole velocity.

There are related fields for which differential-difference equations (similar to discrete drift-diffusion models) model the systems of interest. Well-known are propagation of nerve impulses along myelinated fibers, modelled by discrete FitzHugh-Nagumo equations , motion of dislocations  and sliding charge density waves , modelled by variants of the Frenkel-Kontorova model , etc. The theory of wavefront propagation has been developed for some of these models, which are simpler than ours: convection is typically absent from them and diffusion is purely linear .

The rest of the paper is as follows. We write the drift-diffusion model with appropriate boundary conditions in Section II. There we render these equations dimensionless and explain the results of numerical simulations on a current biased infinitely long SL. Furthermore, we find by numerical simulations that our results for infinite SL may be realized in finite SL with appropriate boundary conditions under constant current bias. The theoretical analysis based on the comparison principle is presented in Section III. Section IV contains our conclusions. Finally some material of a more technical nature is relegated to the Appendices.

## II Discrete drift-diffusion model

### A Equations and boundary conditions

At low enough temperatures (much less than a typical Fermi energy of a SL well measured from the first subband, say 20 meV or 232 K), the following discrete drift-diffusion equations model sequential vertical transport in a weakly doped SL :

$`{\displaystyle \frac{\epsilon }{e}}{\displaystyle \frac{dF_i}{dt}}+{\displaystyle \frac{n_iv(F_i)}{d+w}}D(F_i){\displaystyle \frac{n_{i+1}n_i}{(d+w)^2}}=J(t),`$ (1)
$`F_iF_{i1}={\displaystyle \frac{e}{\epsilon }}(n_iN_D^w).`$ (2)
Eq. (1) is Ampère’s law establishig that the total current density, $`eJ`$, is sum of displacement and tunneling currents. The latter consists of a drift term, $`en_iv(F_i)/(d+w)`$, and a diffusion term, $`eD(F_i)(n_{i+1}n_i)/(d+w)^2`$. We have adopted the convention (usual in this field) that the current density has the same direction as the flow of electrons. Eq. (1) holds for $`i=1,\mathrm{},N1`$. Eq. (2) is the Poisson equation, and it holds for $`i=1,\mathrm{},N`$. $`n_i`$ is the 2D electron number density at well $`i`$, which is singularly concentrated on a plane located at the end of the well. $`F_i`$ is minus an average electric field on a SL period comprising the $`i`$th well and the $`i`$th barrier (well $`i`$ lies between barriers $`i1`$ and $`i`$: barriers 0 and $`N`$ separate the SL from the emitter and collector contact regions, respectively). Parameters $`\epsilon `$, $`d`$, $`w`$, and $`N_D^w`$ are well permittivity, barrier width, well width and 2D doping in the wells, respectively.

Drift velocity and diffusion coefficient are depicted in Fig. 1 for the 9nmGaAs/4nmAlAs SL of Ref. . We have obtained them from microscopic calculations presented in Ref.  (which is appropriate for these sample parameters ) by setting $`v(F)=J(N_D^w,N_D^w,F)(d+w)/N_D^w`$ and $`D(F)=(J(N_D^w,N_D^w,F)/n_{i+1})(d+w)^2`$. Here $`eJ(n_i,n_{i+1},F_i)`$ is the tunneling current between wells $`i`$ and $`i+1`$, $`J_{ii+1}`$. We assume that the tunneling current is a function of the average field at the $`i`$th SL period, $`F_i=F`$, and of the 2D electron densities at wells $`i`$ and $`i+1`$, $`n_i`$ and $`n_{i+1}`$, respectively. Notice that our model for the tunneling current,

$`eJ(n_i,n_{i+1},F_i)={\displaystyle \frac{en_iv(F_i)}{d+w}}eD(F_i){\displaystyle \frac{n_{i+1}n_i}{(d+w)^2}}`$ (3)
$`{\displaystyle \frac{en_iv^{(f)}(F_i)en_{i+1}v^{(b)}(F_i)}{d+w}},`$ (4)
is reasonable for temperatures much lower than a typical Fermi energy in the wells measured from the first subband (say 20 meV), . The tunneling current density should change sign if we reverse the electric field and exchange the electron densities at wells $`i`$ and $`i+1`$: $`J(n_i,n_{i+1},F_i)=J(n_{i+1},n_i,F_i)`$. This inversion symmetry implies

$$v^{(f)}(F)=v^{(b)}(F)\text{and}v(F)=v(F),$$

where $`v^{(b)}(F)=D(F)/(d+w)`$ and $`v^{(f)}(F)=v(F)+v^{(b)}(F)`$. See Figure 1(d).

Equations (1) and (2) should be supplemented with appropriate bias, initial and boundary conditions. Among possible bias conditions, we shall consider the extreme cases of current bias ($`J(t)`$ specified) and voltage bias:

$`(d+w){\displaystyle \underset{i=1}{\overset{N}{}}}F_i=V,`$ (5)
with specified $`V=V(t)`$. Using (5) ignores potential drops at the contact regions and at barrier 0, and it overestimates the contribution of barrier $`N`$ by a factor $`1+w/d`$ . These contributions are negligible for long SL ($`N=40`$ or larger), so that we shall adopt the simpler expression (5). Appropriate boundary conditions have been derived under the same approximations as in (1) . They are

$`{\displaystyle \frac{\epsilon }{e}}{\displaystyle \frac{dF_0}{dt}}+j_e^{(f)}(F_0){\displaystyle \frac{n_1w^{(b)}(F_0)}{d+w}}=J(t),`$ (6)
$`{\displaystyle \frac{\epsilon }{e}}{\displaystyle \frac{dF_N}{dt}}+{\displaystyle \frac{n_Nw^{(f)}(F_N)}{d+w}}=J(t),`$ (7)
where the emitter current density, $`ej_e^{(f)}(F)`$, the emitter backward velocity, $`w^{(b)}(F)`$, and the collector forward velocity, $`w^{(f)}(F)`$ are functions of the electric field depicted in Fig. 3 of Ref.  for contact regions similar to those used in experiments .

To analyze the discrete drift-diffusion model, it is convenient to render all equations dimensionless. Let $`v(F)`$ reach its first positive maximum at $`(F_M,v_M)`$. We adopt $`F_M`$, $`N_D^w`$, $`v_M`$, $`v_M(d+w)`$, $`eN_D^wv_M/(d+w)`$ and $`\epsilon F_M(d+w)/(eN_D^wv_M)`$ as the units of $`F_i`$, $`n_i`$, $`v(F)`$, $`D(F)`$, $`eJ`$ and $`t`$, respectively. For the first plateau of the 9/4 SL of Ref. , we find $`F_M=6.92`$ kV/cm, $`N_D^w=1.5\times 10^{11}`$ cm-2, $`v_M=156`$ cm/s, $`v_M(d+w)=2.03\times 10^4`$ cm$`{}_{}{}^{2}/`$s and $`eN_D^wv_M/(d+w)=2.88`$ A/cm2. The units of current and time are 0.326 mA and 2.76 ns, respectively. Then (1) to (5) become

$`{\displaystyle \frac{dE_i}{dt}}+v(E_i)n_iD(E_i)(n_{i+1}n_i)=J,`$ (8)
$`E_iE_{i1}=\nu (n_i1),`$ (9)
$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}E_i=\varphi .`$ (10)
Here we have used the same symbols for dimensional and dimensionless quantities except for the electric field ($`F`$ dimensional, $`E`$ dimensionless). The parameters $`\nu =eN_D^w/(\epsilon F_M)`$ and $`\varphi =V/[F_MN(d+w)]`$ are dimensionless doping and average electric field (bias), respectively. For the 9/4 SL, $`\nu 3`$. We recall that $`i=1,\mathrm{},N1`$ in (8) and $`i=1,\mathrm{},N`$ in (9). The boundary conditions (6) and (7) become

$`{\displaystyle \frac{dE_0}{dt}}+J_e(E_0)w_e(E_0)n_1=J,`$ (11)
$`{\displaystyle \frac{dE_N}{dt}}+w_c(E_N)n_N=J,`$ (12)
where

$`J_e(E_0)={\displaystyle \frac{j_e^{(f)}(F_ME_0)(d+w)}{N_D^wv_M}},`$ (13)
$`w_e(E_0)={\displaystyle \frac{w^{(b)}(F_ME_0)}{v_M}},`$ (14)
$`w_c(E_N)={\displaystyle \frac{w^{(f)}(F_ME_N)}{v_M}}.`$ (15)
Figure 2 shows $`J_e`$, $`w_e`$ and $`w_c`$ as functions of the electric field. They are dimensionless versions of the curves plotted in Figure 3 of Ref. .

### B Numerical simulations

Simple solutions of the drift-diffusion equations (8) - (9) under constant current bias are stationary or moving monopole wavefronts connecting two electric field domains. Let us consider monopole solutions with profiles $`\{E_i\}`$ which are increasing functions of $`i`$, for they are compatible with realistic boundary conditions in which the emitter region is highly doped . We have simulated numerically on a large SL,

$`{\displaystyle \frac{dE_i}{dt}}{\displaystyle \frac{D(E_i)+v(E_i)}{\nu }}(E_{i1}E_i)`$ (16)
$`{\displaystyle \frac{D(E_i)}{\nu }}(E_{i+1}E_i)=Jv(E_i),`$ (17)
with fixed $`J`$, which is equivalent to (8) - (9). Let $`E^{(1)}(J)<E^{(2)}(J)<E^{(3)}(J)`$, be the three solutions of $`v(E)=J`$ for $`v_m<J<1`$, where $`(E_m,v_m)`$ is the minimum of $`v(E)`$ for $`E>1`$. For the 9/4 SL of Fig. 1, $`E_m=9.8571`$, $`v_m=0.02192`$. We have simulated (17) for different values of $`\nu >0`$ and of $`J(v_m,1)`$. The initial condition was chosen so that $`E_iE^{(1)}(J)`$ as $`i\mathrm{}`$, and $`E_iE^{(3)}(J)`$ as $`i\mathrm{}`$. We observed that, after a short transient, a variety of initial conditions sharing these features evolved towards either a stationary or moving monopole. For systematic numerical studies, we therefore adopted an initial step like profile, with $`E_i=E^{(1)}(J)`$ for $`i<0`$, $`E_i=E^{(3)}(J)`$ for $`i>0`$ and $`E_0=E^{(2)}(J)`$. The boundary data were taken to be $`E_N=E^{(1)}(J)`$, $`E_N=E^{(3)}(J)`$ with $`N`$ large.

Our results show that the dimensionless doping $`\nu `$ determines the type of solution of (17) which is stable. There are two important values of $`\nu `$, $`\nu _1<\nu _2`$.

* For $`0<\nu <\nu _1`$ and each fixed $`J(v_m,1)`$, only traveling monopole fronts moving downstream (to the right) were observed. For $`\nu >\nu _1`$, stationary monopoles were found. According to the arguments of Wacker et al  for the discrete drift model with $`D(E)=0`$, stationary monopoles exist for dimensionless doping larger than a critical value. An upper bound for this critical doping is

$`\nu _c=v_m{\displaystyle \frac{E_m1}{1v_m}},`$ (18)
which equals $`\nu _c=0.198`$ for our numerical example. We have found that $`\nu _1=0.16`$. This agreement with results obtained assuming $`D(E)=0`$ is not surprising: we shall prove in Section III that (18) holds as well for the model (8) - (9) with nonzero diffusivity.
* For $`\nu _1<\nu <\nu _2`$, traveling fronts moving downstream exist only if $`J(v_m,J_1(\nu ))`$, where $`J_1(\nu )<1`$ is a critical value of the current. If $`J(J_1(\nu ),1)`$, the stable solutions are steady fronts (stationary monopoles). We have found that $`\nu _2=0.33`$.
* New solutions are observed for $`\nu >\nu _2`$. As before, there are traveling fronts moving downstream if $`J(v_m,J_1(\nu ))`$, and stationary monopoles if $`J(J_1(\nu ),J_2(\nu ))`$, $`J_2(\nu )<1`$ is a new critical current. For $`J_2(\nu )<J<1`$, the stable solutions of (17) are monopoles traveling upstream (to the left). As $`\nu `$ increases, $`J_1(\nu )`$ and $`J_2(\nu )`$ approach $`v_m`$ and 1, respectively. Thus stationary solutions are found for most values of $`J`$ if $`\nu `$ is large enough.

Figure 3 depicts $`J_1(\nu )`$ and $`J_2(\nu )`$ as functions of $`\nu `$. Notice that $`J_1`$ decreases from $`J_1=1`$ to $`J_1=v_m`$ as $`\nu `$ increases from $`\nu _1`$. Similarly, $`J_2`$ decreases from $`J_2=1`$ to a minimum value $`J_20.53`$ and then increases back to $`J_2=1`$ as $`\nu `$ increases. Monopole velocity as a function of current has been depicted in Figure 4 for four different doping values, $`\nu =0.5`$, $`\nu =1`$, $`\nu =3`$ and $`\nu =10`$. For larger $`\nu `$, the interval of $`J`$ for which stationary solutions exist becomes wider again, trying to span the whole interval $`(v_m,1)`$ as $`\nu \mathrm{}`$. For very large $`\nu `$, the velocities of downstream and upstream moving monopoles become extremely small in absolute value.

Notice that if we use the complete sequential tunneling current instead of the drift-diffusion approximation (4) in Eq. (1), the situation is the same. Figure 5 depicts monopole velocity versus current for well doping corresponding to the 9/4 SL of Ref. . Results obtained with the complete sequential tunneling current or with approximation (4) (corresponding to Fig. 4 with $`\nu =3`$) are compared. Both velocity curves are similar, and their quantitative discrepancies are irrelevant in view of the uncertainties involved in a theoretical calculation of the tunneling current (typically the off-resonance current is larger than the theoretical prediction).

Once different stable monopole solutions (moving either downstream or upstream, stationary) have been identified, we raise the natural question of whether they are compatible with boundary conditions. Another series of numerical simulations was carried out to answer this. We solved numerically (17) for a current biased finite SL ($`N=40`$) with boundary conditions (11) - (15). Our results are depicted in Figure 6 for realistic doping at the contact layers. We observe that the emitter boundary condition results in the creation of a charge accumulation layer near this contact. A charge depletion layer is formed near the collector contact as a result of the corresponding boundary condition. Except for these layers, existence and configuration of monopoles moving downstream, upstream or remaining stationary, agrees with the previous simulations (corresponding to an infinitely long current-biased SL with a monopole-like initial condition).

## III Mathematical analysis of traveling monopoles and stationary solutions

In this Section, we study theoretically moving or stationary monopoles on an infinitely long, current-biased SL. Our findings will confirm the picture suggested by the numerical simulations of the previous Section for any doped weakly coupled SL. Furthermore, we shall prove stability of the different monopole solutions and find bounds for the critical values of $`\nu `$ and $`J_i`$. Our results are based upon and extend ideas first proposed by J.P. Keener for discrete FitzHugh-Nagumo equations, corresponding to signal transmision in myelinated neurons . Mathematically analogous problems arise in models of propagation of defects in crystals . These problems have the following structure,

$`{\displaystyle \frac{dE_i}{dt}}d(E_{i+1}2E_i+E_{i1})=Jv(E_i),`$ (19)
which is much simpler than (17). Here the parameter $`d>0`$ is a constant diffusion coefficient, and $`v(E)`$ a ‘cubic’ function with three branches as the electron drift velocity of Fig. 1.

For (19), there are critical values of $`J`$, $`J_1`$ and $`J_2`$, characterizing wavefront behavior . For $`J>J_2(d)`$, there exist wavefront solutions of (19) moving upstream (to the left). For $`J<J_1(d)`$, there are wavefronts moving downstream (to the right), whereas for $`J_1(d)<J<J_2(d)`$, stationary fronts exist. The width of the interval $`(J_1(d),J_2(d))`$ is an increasing function of $`d`$.

### A Propagation failure and stationary solutions

In Appendix A we state and prove a comparison principle for (17). As a consequence, if our initial field profile is monopole-like \[monotone increasing with well index, and sandwiched between $`E^{(1)}(J)`$ and $`E^{(3)}(J)`$\], so is the electric field profile for any later time $`t>0`$; see Appendix A:

$`\{E_i(0)\}\text{increasing with }i`$ (20)
$`E^{(1)}(J)<E_i(t)<E_{i+1}(t)<E^{(3)}(J),i,t>0.`$ (21)
We now obtain sufficient conditions for an initial monopole not to propagate upstream or downstream. Under these conditions, the monopole may remain stationary or move downstream or upstream, respectively. Let us start with a condition pinning the left tail of a monopole. As $`E_{i1}<E_i`$ and $`E_{i+1}<E^{(3)}(J)`$, we have:

$`{\displaystyle \frac{dE_i}{dt}}`$ $`=`$ $`{\displaystyle \frac{D(E_i)}{\nu }}(E_{i+1}E_i)`$ (22)
$`+`$ $`{\displaystyle \frac{D(E_i)+v(E_i)}{\nu }}(E_{i1}E_i)+Jv(E_i)`$ (23)
$``$ $`{\displaystyle \frac{D(E_i)}{\nu }}[E^{(3)}(J)E_i]+Jv(E_i)0,`$ (24)
provided there exist $`a_l<b_l`$ such that

$`{\displaystyle \frac{D(E)}{\nu }}[EE^{(3)}(J)]Jv(E),E(a_l,b_l),`$ (25)
and then we choose some initial field, $`E_i(0)(a_l,b_l)`$. The previous inequality then implies $`E_i(t)(E^{(1)}(J),b_l)`$ for all $`t>0`$. This in turn forbids a monopole to move upstream (to the left). We say that condition (25) pins the left tail of the monopole. Whether such $`(a_l,b_l)`$ exist, depends on the parameters $`\nu `$ and $`J`$; see Figure 7.

Let us now pin the right tail of a monopole. As $`E_{i+1}>E_i`$ and $`E_{i1}>E^{(1)}(J)`$, we have

$`{\displaystyle \frac{dE_i}{dt}}`$ $`=`$ $`{\displaystyle \frac{D(E_i)}{\nu }}(E_{i+1}E_i)`$ (26)
$`+`$ $`{\displaystyle \frac{D(E_i)+v(E_i)}{\nu }}(E_{i1}E_i)+Jv(E_i)`$ (27)
$``$ $`{\displaystyle \frac{D(E_i)+v(E_i)}{\nu }}[E^{(1)}(J)E_i]+Jv(E_i)0,`$ (28)
provided there exist $`a_r<b_r`$ such that

$`{\displaystyle \frac{D(E)+v(E)}{\nu }}[EE^{(1)}(J)]Jv(E),`$ (29)
$`E(a_r,b_r),`$ (30)
and we choose some initial field, $`E_i(0)(a_r,b_r)`$. The previous inequality then implies $`E_i(t)(a_r,E^{(3)}(J))`$ for all $`t>0`$. A monopole cannot then move downstream (to the right), and we say that its right tail is pinned. Figure 8 illustrates our arguments: for fields larger than $`a_r`$, the $`E_i`$’s tend to increase above $`a_r`$ toward $`E^{(3)}(J)`$. Then the monopole cannot move downstream. For $`E_i<b_l`$, the fields tend to $`E^{(1)}(J)`$, and the monopole cannot move upstream. As before, the existence of $`(a_r,b_r)`$ depends on the values of $`\nu `$ and $`J`$; see Figure 7.

Figures 7 show the curves $`Jv(E)`$, $`D(E)(EE^{(3)})/\nu `$ and $`[D(E)+v(E)](EE^{(1)})/\nu `$ for $`\nu =3`$ and different values of $`J`$. At $`J=0.08`$, Figure 7(a) shows that there is an interval $`(a_l,b_l)`$ as in (25), but no interval $`(a_r,b_r)`$ as in (30) exist. Then the left tail of a monopole is pinned, but its right tail is free. In this conditions, a monopole may move downstream. Figure 7(b) shows a monopole with both its left and right tails pinned for $`J=0.2`$. Then our theory implies that wavefront propagation fails and a monopole-like stationary solution is stable. Numerical simulations show that there are stationary solutions when $`J[0.09,0.53]`$. Finally, Figure 7(c) shows that, if $`J=0.34`$, the right tail of a monopole is pinned, but not its left tail. Under these conditions a monopole may move upstream. For larger $`\nu `$, the estimates become sharper. For instance, when $`\nu =10`$, (25) and (30) hold for $`J[0.05,0.45]`$. Direct numerical simulations show that stationary solutions exist for $`J[0.04,0.55]`$. Systematic use of these criteria allows us to estimate the critical doping values $`\nu _j`$ and critical current values $`J_i(\nu )`$, $`i=1,2`$ defined in the previous Section; see Appendix B. Instead of looking for $`J_1(\nu )`$ and $`J_2(\nu )`$, it is more convenient to look for their inverse functions, which we may call $`\nu _1(J)`$ and $`\nu _2^\pm (J)`$. According to Figure 3, the inverse function of $`J_2(\nu )`$ is two-valued, and its two branches are $`\nu _2^{}(J)<\nu _2^+(J)`$. We have found the following upper bounds $`\nu _{1b}(J)`$ and $`\nu _{2b}^+(J)`$ for $`\nu _1(J)`$ and $`\nu _2^+(J)`$, respectively:

$`\nu _{1b}(J)=v_m{\displaystyle \frac{E_mE^{(1)}(J)}{Jv_m}},`$ (31)
$`\nu _{2b}^+(J)=D(1){\displaystyle \frac{E^{(3)}(J)1}{1J}}.`$ (32)
If $`\nu >\nu _{1b}(J)`$, the right tail of the monopole is pinned, whereas the left tail of the monopole is pinned if $`\nu >\nu _{2b}^+(J)`$; see Appendix B. Notice that $`\nu _{1b}(J)`$ is a decreasing function of $`J`$. Therefore the critical value $`\nu _1`$ (above which there are stationary solutions) is smaller than $`\nu _{1b}(1)`$, which is exactly Wacker et al’s bound, (18). This explains why the bound (18) gives surprisingly good results even for the first plateau of the SL current-voltage characteristics \[despite having been obtained under the assumption $`D(E)0`$\] . Notice that the bound (32) is reasonable for large dopings and currents $`J1`$. See Figure 3 for a comparison between the critical curves $`J_1(\nu )`$ and $`J_2(\nu )`$ and the bounds (31) and (32).

### B Propagation: traveling fronts

Having shown that only one tail of a monopole is pinned suggests that the monopole may move in the opposite direction. Direct simulations show that this is often the case, and we will prove this now.

An upstream traveling wave solution of (17) may have the form

$`E_i(t)=w(i+ct),c>0.`$ (33)
We will look for an electric field profile $`w(z)`$, $`z=i+ct`$, which is not an exact solution of (17), but instead it satisfies

$`c{\displaystyle \frac{dw}{dz}}{\displaystyle \frac{D(w(z))+v(w(z))}{\nu }}[w(z1)w(z)]`$ (34)
$`{\displaystyle \frac{D(w(z))}{\nu }}[w(z+1)w(z)]+v(w(z))J0.`$ (35)
If this subsolution is initially below an initial field profile, i.e. $`w(i)<E_i(0)`$, for all $`i`$, then the comparison theorem of Appendix A guarantees that $`E_i(t)>w(i+ct)`$ for later times. As $`w(i+ct)`$ moves upstream, so does $`E_i(t)`$, and the electric field profile corresponds to a monopole moving upstream with velocity at least $`c`$. See Figure 9: a subsolution “pushes” the monopole upstream, whereas a supersolution (defined below) “pushes” the monopole downstream.

How do we find a reasonable subsolution? An idea is to try a piecewise continuous solution which equals $`E^{(1)}(J)`$ for $`z<z_0`$ and a larger constant $`A`$, $`A(E^{(2)}(J),E^{(3)}(J))`$ \[and therefore $`v(A)J0`$\], for $`z>z_1`$, with $`z_1>z_0`$. For $`z_0<z<z_1`$, $`w(z)`$ is an unspecified smooth increasing function with $`w(z_0)=E^{(1)}(J)`$ and $`w(z_1)=A`$. Now we shall select conveniently the numbers $`z_0`$, $`z_1`$, $`c`$ and $`A`$ so that (35) holds. Clearly, (35) holds for $`z+1<z_0`$ and for $`z1>z_1`$. Suppose that $`0<z_1z_0<1`$. Then there are five possibilities:

1. $`z<z_0`$, $`z+1>z_1`$. Then $`w(z1)=w(z)=E^{(1)}(J)`$, $`w(z+1)=A`$, which inserted in (35) yields $`D(E^{(1)})(AE^{(1)})/\nu 0`$ (obviously true).
2. $`z1<z_0`$, $`z>z_1`$. Then $`w(z1)=E^{(1)}(J)`$, $`w(z)=w(z+1)=A`$, which inserted in (35) yields

$`Jv(A){\displaystyle \frac{D(A)+v(A)}{\nu }}[AE^{(1)}(J)].`$ (36)
3. $`z_0<z1<z_1`$. Then $`E^{(1)}(J)<w(z1)<A`$ and $`w(z)=w(z+1)=A`$, which yields $`Jv(A)[D(A)+v(A)][Aw(z1)]/\nu `$. This inequality holds if (36) does.
4. $`z_0<z<z_1`$. Then $`w(z1)=E^{(1)}(J)`$, $`E^{(1)}(J)<w(z)<A`$ and $`w(z+1)=A`$. Inserting this in (35), we find

$`c{\displaystyle \frac{dw}{dz}}{\displaystyle \frac{D(w(z))+v(w(z))}{\nu }}[E^{(1)}w(z)]`$ (37)
$`+{\displaystyle \frac{D(w(z))}{\nu }}[Aw(z)]+Jv(w(z)).`$ (38)
Let us now assume that we can select $`A(E^{(2)},E^{(3)})`$ such that the right hand side of this expression is positive, say

$`{\displaystyle \frac{D(w)+v(w)}{\nu }}[E^{(1)}(J)w]+{\displaystyle \frac{D(w)(Aw)}{\nu }}`$ (39)
$`+Jv(w)\delta >0,E^{(1)}<w<A,`$ (40)
and that we choose $`c`$ so that $`cdw/dz<\delta `$. Then (35) holds.
5. $`z_0<z+1<z_1`$. Then $`w(z1)=w(z)=E^{(1)}(J)`$, and $`E^{(1)}(J)<w(z+1)<A`$, which inserted in (35) yields $`D(E^{(1)})[w(z+1)E^{(1)}]/\nu 0`$ (obviously true).

Summarizing the previous arguments, provided (36) and (40) hold, $`w(z)`$ is a subsolution obeying (35). The parameter $`A`$ can be found graphically. First of all, we depict the functions $`Jv(E)`$ and $`f_1(E;J)[D(E)+v(E)][EE^{(1)}(J)]/\nu `$. Possible values of $`A`$ are those $`E`$ for which $`Jv(E)f_1(E;J)`$. For such A, we may plot the left side of (40),

$`f_2(E;J,A){\displaystyle \frac{D(E)+v(E)}{\nu }}[E^{(1)}(J)E]`$ (41)
$`+{\displaystyle \frac{D(E)(AE)}{\nu }}+Jv(E).`$ (42)
If $`f_2(E;J,A)>0`$ for $`E(E^{(1)}(J),A)`$, then the selected value of $`A`$ allows us to construct the sought subsolution. See Figure 10 for a practical realization of this graphical construction.

We have proved rigorously that monopoles may move upstream under favorable circumstances. Our proof using subsolutions may yield a very practical additional bonus: an upper bound, $`c^{}`$, for the velocity of the monopole. Let us choose $`\delta (J,A)=`$ min$`{}_{E^{(1)}<E<A}{}^{}f_{2}^{}(E;J,A)`$, $`z_1z_0=1`$, and $`w(z)=[AE^{(1)}(J)](zz_0)`$ for $`z_0<z<z_1`$. Then $`c^{}=\delta (J,A)/[AE^{(1)}(J)]`$. In Figure 4, $`c^{}`$ is represented by a line of thick dots for doping $`\nu =3`$ corresponding to the 9/4 SL.

In a similar vein, we can construct supersolutions which push the monopole field profile to the right; see Figure 9. Now we start from a monopole profile moving downstream,

$`E_i(t)=w(ict),c>0.`$ (43)
The electric field profile $`w(z)`$, $`z=ict`$ should satisfy

$`c{\displaystyle \frac{dw}{dz}}+{\displaystyle \frac{D(w(z))+v(w(z))}{\nu }}[w(z1)w(z)]`$ (44)
$`+{\displaystyle \frac{D(w(z))}{\nu }}[w(z+1)w(z)]+Jv(w(z))0.`$ (45)
We seek a piecewise continuous supersolution which equals a constant, $`A`$, $`A(E^{(1)}(J),E^{(2)}(J))`$, for $`z<z_0`$, and $`w(z)=E^{(3)}(J)`$ for $`z>z_1`$, with $`z_1>z_0`$. For $`z_0<z<z_1`$, $`w(z)`$ is an unspecified smooth increasing function with $`w(z_0)=A`$, and $`w(z_1)=E^{(3)}(J)`$. As for subsolutions, we now select conveniently the numbers $`z_0`$, $`z_1`$, $`c`$ and $`A`$ so that (45) holds. Clearly, (45) holds for $`z+1<z_0`$ and for $`z1>z_1`$. Suppose that $`0<z_1z_0<1`$. An analysis of the remaining five possibilities yields the following criteria to hold for $`w(ict)`$ to be a supersolution:

$`Jv(A){\displaystyle \frac{D(A)}{\nu }}[E^{(3)}(J)A]f_3(A;J),`$ (46)
$`f_4(w;J,A){\displaystyle \frac{D(w)+v(w)}{\nu }}(wA)`$ (47)
$`+{\displaystyle \frac{D(w)}{\nu }}[E^{(3)}(J)w]+Jv(w)\delta ,`$ (48)
$`\text{for }AwE^{(3)}(J),`$ (49)
$`c{\displaystyle \frac{dw}{dz}}\delta .`$ (50)
Provided such $`w(ict)`$ is found, solutions $`E_i(t)`$ of (17) with $`E_i(0)<w(i+\tau )`$ will satisfy $`E_i(t)<w(ict+\tau )`$ and propagate to the right with speed larger than $`c`$. $`\tau `$ is a constant which can be conveniently chosen to keep the monopole profile below the supersolution. Figure 11 illustrates the graphical construction of the supersolution by checking that (46) and (48) hold for particular values of $`J`$ and $`\nu `$.

As in the subsolution case, an upper bound $`c^{}`$ for the monopole velocity $`c`$ is estimated by choosing $`\delta (J,A)=`$ max$`{}_{A<E<E^{(3)}}{}^{}f_{4}^{}(E;J,A)`$, $`z_1z_0=1`$, and $`w(z)=[E^{(3)}(J)A](zz_0)`$ for $`z_0<z<z_1`$. Then $`c^{}=\delta (J,A)/[E^{(3)}(J)A]`$. In Figure 4, $`c^{}`$ is represented by a line of thick dots for doping $`\nu =3`$ corresponding to the 9/4 SL.

We can now summarize the results obtained from sub and supersolutions; see Figures 10 and 11. We find reasonably good upper bounds for the absolute value of monopole velocity. Furthermore, for $`\nu =3`$, conditions (36) and (40) hold for $`J=0.6`$ and $`A=12`$, whereas conditions (46) and (48) hold for $`J=0.05`$ and $`A=3`$. Therefore, monopoles move downstream for $`J0.05`$ and they move upstream for $`J0.6`$. Direct numerical simulations show that: (i) the estimate $`J_1=0.05`$ for the first critical current can be improved to $`J_1=0.08`$; and (ii) $`J_2=0.6`$ for the second critical current can be improved to $`J_2=0.54`$.

## IV Conclusions and final comments

We have presented a theory of monopoles moving downstream or upstream on an infinitely long doped, current biased superlattice when the fields are on the first plateau of the current–voltage characteristic. This theory has been corroborated with numerical evidence, which sharpens our results. Furthermore, we have simulated a 40-well 9nmGaAs/4nmAlAs SL  under doping and contact conditions similar to experimental ones , but under constant current bias conditions. This situation is different from the usual case of voltage bias conditions. We have obtained that it is possible to observe monopole wavefronts moving upstream when the current is kept at large enough levels. Together with our theoretical bounds for critical currents and dopings, this numerical prediction could be used to set up an experiment to observe this striking phenomenon. For this purpose, we would need an initial condition corresponding to a monopole separating two electric field domains at high enough current. In an ideal world, this situation could be obtained by first fixing a low dc voltage for the 9/4 sample at a value near the top of one of the first branches of the current - voltage characteristics. Then we could switch from voltage to current bias conditions. The outcome would be a monopole moving upstream until the emitter region is reached. Presumably an idea of the field distribution corresponding to this situation could be obtained by time-resolved photoluminescence measurements .

There are technical problems that must be overcome if one wants to observe these features in real experiments: when we switch, there will always be a Faraday-like inductive pulse which will probably perturb the state of the system in an uncontrolled way. There are other possible biases we could think of. Under dc voltage bias, upstream moving monopoles are probably created for a short time during relocation experiment . In these experiments, one has a doped SL with a current-voltage characteristics correponding to multiple stationary monopole solution branches. Voltage is set at a particular value near the end of a branch, so that the field profile is that of a monopole layer connecting a low to a high field domain. Let the monopole layer be located at well $`i`$ (counted from the emitter contact). Then the voltage is suddenly and appropriately increased. After a certain time, the field profile settles to a new situation corresponding to a monopole layer centered at well $`i1`$ . This could be an indication of a monopole moving upstream, albeit for a short time. To increase this time, we could try to set a hybrid bias (between current and voltage bias) by including a finite series resistance in our external circuit. Additional theoretical and numerical work is needed to explore these possibilities.

## ACKNOWLEDGMENTS

We thank A. Amann, G. Platero and D. Sánchez for fruitful discussions and collaboration on the discrete drift-diffusion model. LLB thanks S.W. Teitsworth for a critical reading of the manuscript and helpful comments. This work was supported by the Spanish DGES through grant PB98-0142-C04-01, and by the European Union TMR contracts ERB FMRX-CT96-0033 and ERB FMRX-CT97-0157.

## A Comparison principle

The main theorem which we use to prove our results in Section III is the following comparison principle:

Theorem A.1 Let $`U_i(t)`$ and $`L_i(t)`$, $`iZ`$, be differentiable sequences such that

$`{\displaystyle \frac{dU_i}{dt}}d_1(U_i)[U_{i+1}U_i]`$ (A1)
$`d_2(U_i)[U_{i1}U_i]f(U_i)`$ (A2)
$`{\displaystyle \frac{dL_i}{dt}}d_1(L_i)[L_{i+1}L_i]`$ (A3)
$`d_2(L_i)[L_{i1}L_i]f(L_i),`$ (A4)
$`U_i(0)>L_i(0).`$ (A5)
where $`f`$, $`d_1>0`$ and $`d_2>0`$ are Lipschitz continuous functions. Then,

$$U_i(t)>L_i(t),t>0,iZ$$

In our discrete drift-diffusion model,

$`d_1(E)={\displaystyle \frac{D(E)}{\nu }},d_2(E)={\displaystyle \frac{D(E)+v(E)}{\nu }},`$ (A6)
$`f(E)=Jv(E).`$ (A7)
Proof: The proof is by contradiction. Set $`W_i(t)=U_i(t)L_i(t)`$. At $`t=0`$, $`W_i(0)>0`$ for all $`i`$. Let us assume that $`W_i`$ changes sign after a certain minimum time $`t_1>0`$, at some value of $`i`$, $`i=k`$. Thus $`W_k(t_1)=0`$ and $`dW_k/dt0`$, as $`tt_1`$. We shall show that this is contradictory. At $`t=t_1`$, there must be an index $`m`$ (equal or different from $`k`$) such that $`W_m(t_1)=0`$, while its next neighbor $`W_{m+j}(t_1)>0`$ ($`j`$ is either 1 or -1), and $`W_i(t_1)=0`$ for all indices between $`k`$ and $`m`$. For otherwise $`W_k`$ should be identically 0 for all $`k`$. Equation (A4) implies

$`{\displaystyle \frac{dW_m}{dt}}(t_1)d_1(U_m(t_1))W_{m+1}(t_1)`$ (A8)
$`+d_2(U_m(t_1))W_{m1}(t_1)>0.`$ (A9)
This contradicts the fact that $`dW_m/dt`$ should have been nonpositive as $`tt_1`$, for $`W_m(t_1)`$ to have become zero in the first place.

Corollary A.1 Any solution $`E_i(t)`$ of (17) with initial data $`E_i(0)(E^{(1)}(J),E^{(3)}(J))`$ satisfies $`E_i(t)(E^{(1)}(J),E^{(3)}(J))`$ for $`t>0`$.

Proof: Apply Theorem A.1 first with $`L_i=E^{(1)}(J)`$ and $`U_i=E_i`$, then with $`L_i=E_i`$ and $`U_i=E^{(3)}(J)`$.

Corollary A.2 If $`E_i(0)`$ is monotone increasing, that is, $`E_i(0)<E_{i+1}(0)`$, then, $`E_i(t)`$ is also monotone increasing, i.e., $`E_i(t)<E_{i+1}(t)`$ for $`t>0`$.

Proof: Apply Theorem A.1 with $`L_i=E_i(t)`$ and $`U_i=E_{i+1}(t)`$.

Remark. Strict inequalities in these theorems can be replaced by inequalities and the corresponding statements still hold. However the proofs become rather more technical and involved.

## B Bounds for critical doping values

We want to estimate the curves $`\nu _{1b}(J)`$ and $`\nu _{2b}^\pm (J)`$ defined in Section III. To estimate $`\nu _{2b}^+(J)`$, assume that $`J1`$ and $`\nu `$ is large. The left tail of a monopole is pinned if Eq. (25) holds. For large currents, (25) certainly holds if the curve corresponding to the left side of the inequality is above that of the right hand side, for $`E=1`$ (this is possible because $`D(E)`$ decreases rapidly to zero as the field increases). Setting $`E=1`$ in (25), we obtain

$$\frac{D(1)[1E^{(3)}(J)]}{\nu }>J1.$$

In turn, this implies $`\nu >\nu _{2b}^+(J)`$, defined in (32). This argument fails for the small values of $`J`$ used to draw Figure 7. We believe that quite different reasoning is needed to estimate $`\nu _{2b}^{}(J)`$.

The same argument yields our estimate $`\nu _{1b}(J)`$ of (31). For (30) to hold, the curve corresponding to the left side of the inequality should be below that of the right hand side for $`E=E_m`$. As $`D(E_m)0`$, we obtain

$$\frac{v(E_m)}{\nu }[E_mE^{(1)}(J)]<Jv(E_m),$$

which yields (31). Fig. 3 shows that the bound (31) is reasonably good for all eligible values of $`\nu `$ and $`J`$.


format: text

text:
# Enhancement of GMR due to spin-mixing in magnetic multilayers with a superconducting contact

## Abstract

We study the Giant Magnetoresistance (GMR) ratio in magnetic multiayers with a single superconducting contact in the presence of spin-mixing processes. It has been recently shown  that the GMR ratio of magnetic multilayers is strongly suppressed by the presence of a superconducting contact when spin-flipping is not allowed. In this Letter we demonstrate that the GMR ratio can be dramatically enhanced by spin-orbit interaction and/or non-collinear magnetic moments. The system is described using a tight-binding model with either $`s`$-$`p`$-$`d`$ or $`s`$-$`d`$ atomic orbitals per site.

PACS numbers: 75.70.Pa, 74.80.Dm

Hybrid nanostructures form a fascinating melting pot for studying the interplay between fundamental quantum phenomena, often revealing new and unexpected physics. One recently-recognized class of such structures, involving the coexistence of superconducting contacts and ferromagnetic domains, has led to the identification of a number of fundamental issues , several of which are currently unresolved. In this Letter, we examine one such issue, posed by experiments on giant magnetoresistance (GMR) in magnetic (M) multilayers with superconducting (S) contacts and current-perpendicular-to-the-plane (CPP). Recognizing that the sub-gap conductance of such structures is mediated by Andreev scattering, it was recently noted  that in the absence of spin-flip processes, the conductance of a metallic (i.e. diffusive) multilayer in the presence of aligned magnetic moments is almost identical to that of the multilayer when adjacent moments are anti-aligned and therefore the conventional GMR ratio should be strongly suppressed. Since large GMR ratios are observed experimentally , it is clear that even a qualitative understanding of transport in such structures must incorporate the effects of spin-mixing. The aim of this Letter is to present the first theoretical description of CPP GMR in M-multilayers with S-contacts, which incorporate spin-flip scattering. As sources of spin-mixing we consider both spin-orbit (SO) coupling and non-collinear magnetizations in adjacent magnetic layers.

The system under consideration is a disordered magnetic multilayer consisting of an alternating sequence of magnetic layers each of length $`l_M`$ and non-magnetic layers (N) of length $`l_N`$. The building block of the magnetic structure is the bilayer \[M/N\] of length $`l_B=l_N+l_M`$. The magnetic moments of even-numbered M-layers make an angle $`\theta `$ relative to those of odd-numbered M-layers. Experimentally, $`\theta `$ can be varied by applying an external magnetic field with antiparallel (AP) alignment ($`\theta =\pi `$) typically occurring at zero field and parallel (P) alignment ($`\theta =0`$) at large enough fields. The current flows perpendicular to the planes of the multilayer, which makes contact with a metallic normal lead on the left-hand side of the multilayer and a superconducting lead on the right-hand side. GMR is the drastic increase in electrical conductance $`G(\theta )`$ that occurs when the system switches from the AP to the P alignment with the conventional GMR ratio defined by: $`\rho =\frac{G(0)G(\pi )}{G(\pi )}`$.

Following , the multilayer and leads are modelled using a tight-binding Hamiltonian on a cubic lattice with hoppings to nearest neighbours. Lattice imperfections and impurities are simulated by adding to the on-site energies a random number in the range $`[\frac{W}{2},+\frac{W}{2}]`$. The on-site Hamiltonian has the following structure:

$$H=\left(\begin{array}{cccc}H^p& \mu _{xy}& \mathrm{\Delta }& 0\\ \mu _{xy}^{}& H^p& 0& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& 0& H^h& \mu _{xy}^{}\\ 0& \mathrm{\Delta }^{}& \mu _{xy}& H^h\end{array}\right)$$
(1)
where $`H^{p()}`$ is the Hamiltonian for up (down)-spin particles ($`s`$ and $`d`$ bands), $`H^{h()}=H^{p()}`$ is the Hamiltonian for up(down)-spin holes and $`\mathrm{\Delta }`$ is the superconducting order parameter. Here $`\mu _{xy}=\mu _x+i\mu _y`$, where $`\mu _{x(y)}`$ is the $`x(y)`$-component of the exchange field $`\stackrel{}{\mu }`$. Note that $`\stackrel{}{\mu }`$ is non-zero only for electrons in the $`d`$-band in the M-layers and $`\mathrm{\Delta }`$ is non-zero only in the right-hand-side superconducting lead. Within the tight-binding formulation, SO interaction can be included by adding to the Hamiltonian the following term:

$$V_{SO}=V_1\underset{i,j,\alpha ,s}{}\stackrel{}{\sigma }\stackrel{}{R}_{i,j}c_{\alpha ,i}^\sigma c_{\alpha ,j}^\sigma $$
(2)
where $`V_1`$ is a constant which determines the interaction strength, $`\stackrel{}{\sigma }`$ is a vector of Pauli matrices, $`\stackrel{}{R}_{i,j}`$ is the unit vector which connects site $`i`$ with the neighbouring site $`j`$. $`c_{\alpha ,i}^\sigma `$ is the annihilation operator for electrons of spin $`\sigma `$ in the $`\alpha `$ ($`s`$, $`d`$)-band on site $`i`$. In the presence of disorder (2) produces spin-flip scattering since it couples electrons with different spin on neighbouring sites.

In the presence of disorder, to study the largest possible sample cross sections, we consider 2 orbitals per site, which is the minimal model capable of reproducing scattering potential at the N/M interface and interband scattering . The tight-binding parameters are chosen to reproduce the GMR ratio and conductances obtained from an ab initio material specific calculation for Cu/Co multilayers .

In the presence of spin-flip scattering the two-spin-fluid approximation does not hold and the zero-temperature, zero-bias, normal-state Landauer formula takes the form:

$$G^{NN}=\frac{e^2}{h}\underset{\sigma \sigma ^{}}{}\text{Tr}\left\{t^{\sigma \sigma ^{}}t^{\sigma \sigma ^{}}\right\}$$
(3)
where $`t^{\sigma \sigma ^{}}`$ is the matrix of transmission amplitudes for injected $`\sigma ^{}`$-spin electrons in the left-hand lead into $`\sigma `$-spin electrons in the right-hand lead. When the right-hand lead is in the superconducting state, the conductance is given by

$$G^{NS}=\frac{e^2}{h}2\underset{\sigma \sigma ^{}}{}\text{Tr}\left\{r_a^{\sigma \sigma ^{}}r_a^{\sigma \sigma ^{}}\right\}$$
(4)
where $`r_a^{\sigma \sigma ^{}}`$ is the Andreev reflection matrix for injected $`\sigma ^{}`$-spin electrons in the left-hand lead to be reflected into $`\sigma `$-spin holes. In what follows, the scattering amplitudes are calculated exactly by solving the Bogoliubov-de Gennes equation using an efficient recursive Green’s function technique .

We shall now turn to the central results of this Letter, namely that in the presence of strong-enough spin-mixing, either produced by SO coupling or non-collinear moments, GMR in the presence of a S-contact approaches that of two normal contacts, with values of $`\rho `$ of the order of 100 %. First consider the effect of SO coupling within the multilayer. Fig. 1 shows the conventional GMR ratio as a function of the SO interaction strength $`V_1`$ for a disordered multilayer of 40 bilayers, with $`l_M=15`$ and $`l_N=8`$. As expected, in the NN case $`\rho `$ decreases monotonically from $`200`$ % at $`V_1=0`$, to zero at large $`V_1`$ ($`0.17`$ eV). In contrast for the NS case, $`\rho `$ initially increases with increasing $`V_1`$, eventually joining the NN curve at $`V_10.08`$ eV. As a second source of spin-mixing, consider the effect of non-collinear magnetic moments when $`V_1=0`$. Fig. 2 shows the $`\theta `$-dependence of the GMR ratio defined as $`\rho (\theta )=\frac{G(\theta )G(\pi )}{G(\pi )}`$. Whereas in the NN case $`\rho (\theta )`$ decreases monotonically with increasing $`\theta `$, in the NS case $`\rho (\theta )`$ exhibits a pronounced maximum around $`\theta =\pi /8`$.

To understand these results, first consider the case of non-collinear moments. In the presence of two normal-metallic contacts, the conductance $`G(\theta )`$ has been theoretically studied in  where it is predicted that $`G(\theta )G(\pi )`$ tends monotonically to zero as $`\theta `$ varies from $`0`$ to $`\pi `$. In addition, the dependence of the resistance on the angle $`\theta `$ has been experimentally found  to contain a term proportional to cos$`{}_{}{}^{2}(\theta /2)`$ and a second term proportional to cos$`{}_{}{}^{4}(\theta /2)`$. In the presence of a S-contact, where $`G(0)G(\pi )`$, this behaviour is drastically changed by the presence of an extremum which occurs at some intermediate angle $`\theta _c`$, the value of which depends on the interplay between competing effects. Since $`\theta (H)`$ is a function of the applied magnetic field $`H`$, the presence of a S-contact introduces a new characteristic field $`H_c`$ for which $`\theta (H_c)=\theta _c`$. For a disordered multilayer of 22 bilayers, with $`l_M=30`$, $`l_N=16`$, the insert in Fig. 2 shows the $`\theta `$-dependence of the conductance divided by the number of open channels in the normal lead. As expected $`G^{NN}(\theta )`$ is a monotonic function of $`\theta `$, whereas $`G^{NS}(\theta )`$ possesses an extremum at $`\theta _c\pi /8`$. To understand why the extremum is a maximum, recall that for $`\theta =0`$ or $`\theta =\pi `$, when spin is conserved, current flows when a right-going (spin $`\sigma `$) electron passes through the multilayer, Andreev reflects as a left-going (spin $`\sigma `$) hole, which retraverses the multilayer. A M-layer whose moment is aligned with the spin of the incident electron is anti-aligned with the spin of the outgoing hole and consequently the number of aligned and anti-aligned M-layers encountered by a given quasi-particle is the same for both $`\theta =0`$ and $`\theta =\pi `$ (only the order differs). When the elastic mean free path is comparable with the total multilayer length, the resistance of traversed layers add in series and therefore, apart from small differences due to interference effects , $`G^{NS}(0)G^{NS}(\pi )`$. Furthermore, since a quasi-particle must necessarily traverse regions in which it is a minority spin, both $`G^{NS}(0)`$ and $`G^{NS}(\pi )`$ are low-conductance states. In contrast, as $`\theta `$ increases from zero, this conductance bottleneck is removed, because an Andreev reflected minority hole can spin-convert to a majority hole, thereby avoiding anti-aligned moments on its return journey. Of course this initial increase in $`G^{NS}(\theta )`$ is eventually overcome by the usual GMR effect which decreases $`G^{NS}(\theta )`$ as $`\theta \pi `$, thereby producing an overall maximum.

In the absence of disorder, the nature of the extremum is determined by interface scattering and band structure. To illustrate this consider a clean multilayer which is perfectly periodic and therefore the variation of the conductance with $`\theta `$ arises from tuning of the ballistic spin-filtering by the structure. Fig. 3 shows the conductance divided by the number of open channels in the left-hand-side normal lead. As expected, $`G^{NN}(\theta )`$ is a monotonic function of $`\theta `$, whereas $`G^{NS}(\theta )`$ exhibits a minimum around $`\pi /2`$ and then increases. (In this case, translational invariance in the transverse direction allowed us to use a full ab initio, $`spd`$ Hamiltonian to obtain the results of Fig. 3.) In the NN case, the dependence of the multilayer resistance on $`\theta `$ predicted by our model is in good agreement with experiment . In Ref.  the ratio between the resistance at a given $`\theta `$ and the resistance with AP alignment has been found to fit the following function:

$$\frac{R(\theta )}{R(\pi )}=1a\mathrm{cos}^2(\theta /2)+b\mathrm{cos}^4(\theta /2)$$
(5)
where $`a`$ and $`b`$ are fitting constants. In Fig. 4 we show the plot of such a ratio for the disordered multilayer considered above, along with the best fit to function (5). In addition we also checked that this ratio cannot be fitted with the same accuracy assuming a pure dependence on cos$`{}_{}{}^{2}(\theta /2)`$ (i.e. with $`b=0`$). For $`G^{NS}(\theta )`$ however, no such analytic results currently exist.

Let us now turn attention to the effect of SO coupling. Figs. 5a and 5b show the conductances as a function of the SO strength $`V_1`$ for, respectively, the NN and the NS case. In the NN case (Fig. 5a) $`G_P^{NN}`$ decreases as $`V_1`$ increases and eventually joins the curve for $`G_{AP}^{NN}`$. This can be understood in terms of the heuristic model presented in Ref. , because, as the SO strength increases, the average length required for a spin to flip (spin relaxation length $`\lambda _{sf}`$) gets shorter. Therefore in the P alignment, an injected majority electron travels through the multilayer for a length $`\lambda _{sf}`$ before being scattered into a minority spin, thereby producing a decrease in the conductance. This suggests that the value of $`V_1`$ for which $`G_P^{NN}G_{AP}^{NN}`$ corresponds to a spin relaxation length $`\lambda _{sf}`$ close to the period $`l_B`$ of the multilayer. We have carried out a range of simulations which show that this value of $`V_1`$ does not depend on the overall length of the multilayer, but decreases with increasing $`l_B`$. As expected, the conductance with AP alignment does not change significantly with $`V_1`$.

In the NS case (Fig. 5b) the conductance in the P aligned state rapidly increases with $`V_1`$, reaching a maximum and thereafter decreases, eventually joining the curve for the AP configuration. Clearly the enhancement in $`G_P^{NS}`$ is produced by the onset of spin-flip scattering. The abrupt increase is understandable, since even a relatively small probability for spin flipping opens a highly conductive “channel” if the spin-flip events take place in the vicinity of the interface. As one can see in the insert of Fig. 1 for larger values of $`V_1`$, the conductance $`G_P^{NS}`$ joins $`G_P^{NN}`$ and together they decrease thereafter. As in the normal case $`G_{AP}^{NS}`$ depends weakly on $`V_1`$. The value $`V_10.08`$, at which $`G^{NS}`$ is maximum, corresponds to a spin relaxation length close to the total length of the multilayer and, as expected, separate simulations show that this value of $`V_1`$ decreases with increasing total length. Similarly the value of $`V_1`$ at which the GMR ratio (of Fig. 1) vanishes corresponds to a spin-relaxation length of the order the bilayer thickness $`l_B`$ and is independent of the total length of the system.

In conclusion, we have demonstrated that spin-mixing plays a crucial rôle in determining both the qualitative and quantitative features of GMR in magnetic multilayers with a S-contact. In contrast with the normal case, where spin-mixing suppresses GMR, we find that the GMR ratio can be dramatically enhanced by the presence of spin-orbit interactions and/or non-collinear magnetic moments. In experiments carried out to-date, the presence of large spin-orbit scattering  presumably masks the mechanism shown in Fig. 2, which is predicted to be a generic feature in the absence of other spin-mixing processes. This suggests that lighter metals and superconductors would be more appropriate for observing the new extrema predicted in this Letter. Finally we note that for the future it would be of interest to examine spin-mixing in non-diffusive NS structures such as clean spin-valves , where the GMR ratio can be non-zero or negative, even in the absence of spin-flip processes.


text:
# Three-dimensional vortex dynamics in Bose-Einstein condensates

## Abstract

We simulate in the mean-field limit the effects of rotationally stirring a three-dimensional trapped Bose-Einstein condensate with a Gaussian laser beam. A single vortex cycling regime is found for a range of trap geometries, and is well described as coherent cycling between the ground and the first excited vortex states. The critical angular speed of stirring for vortex formation is quantitatively predicted by a simple model. We report preliminary results for the collisions of vortex lines, in which sections may be exchanged.

The formation and properties of vortices in Bose-Einstein condensates are attracting intense current interest. A number of schemes for producing vortices in trapped condensates have been proposed theoretically, and the first experimental creation of a topological vortex in a two-component condensate has been reported , in agreement with a numerical model . Simulations of flow of a condensate past a spherical object have shown vortex ring solutions , and numerical studies of vortex formation in rotating traps have shown the formation of vortex arrays . Very recently the first report has been made of the observation of vortices in a single-component condensate .

In an earlier paper , we showed that a rotationally stirred Bose-Einstein condensate, simulated in two dimensions, exhibits a simple single vortex cycling regime, a behaviour we described as nonlinear Rabi cycling. In the present paper we demonstrate that the single vortex cycling regime is also found in full three-dimensional simulations of condensate stirring (as shown in Fig. 1), in a parameter regime accessible to current experiments. The single vortex regime is studied in a variety of trap geometries (oblate to prolate) as is the dynamical stability of central vortices. The line character of a vortex in three dimensions also gives rise to some qualitatively new behaviour that has no analogue in two dimensions. We show that in a spherical trap, stirring of the condensate can produce two vortex lines which collide and exchange sections.

As in the earlier work, our treatment is based on the time-dependent Gross-Pitaevskii (GP) equation,

$$i\frac{\psi (𝐫,t)}{t}=\left[^2+V(𝐫,t)+C|\psi (𝐫,t)|^2\right]\psi (𝐫,t),$$
(1)
in dimensionless units . The total external potential $`V(𝐫,t)`$ is given by $`(x^2+y^2+\lambda ^2z^2)/4+W(𝐫,t)`$, where the first term represents a cylindrically symmetric trap with trap anisotropy parameter $`\lambda `$. The second term $`W(𝐫,t)`$ is the contribution of the stirrer, a far-blue detuned Gaussian laser beam cylindrically symmetric about a line parallel to the $`z`$ axis,

$`W(𝐫,t)=W_s(t)\mathrm{exp}\left[\left({\displaystyle \frac{|𝝆𝝆_s(t)|}{w_s/2}}\right)^2\right],`$ (2)
where $`𝝆`$ is the projection of $`𝐫`$ into the $`z=0`$ plane, and $`𝝆_s`$ is the centre of the Gaussian stirrer in that plane.

The condensate is stirred in a manner similar to our earlier work , with the stirrer moving in a circle of radius $`\rho _s`$ at a constant angular velocity $`\omega _f`$. Here, however, we begin with the ground state of the harmonic trap (with no stirrer present), and minimise transient effects by increasing the stirrer amplitude $`W_s`$ linearly from zero to a final value of $`W_0`$ between $`t=0`$ and $`t=\pi `$. With an appropriate choice of parameters, as in Fig. 1, a single vortex with positive circulation is drawn close to the centre of the condensate, and for continuous stirring cycles in and out. This is exactly the behaviour found previously in a two-dimensional system , and in fact the parameter choice for the single vortex cycling behaviour in three dimensions is guided by the same considerations outlined there, as we discuss below. Figure 2, where we plot the phase of the wavefunction, shows that the central feature of Fig. 1 is a vortex of positive circulation. In Fig. 3 the evolution of the angular momentum expectation value (which in our dimensionless units is 0 for the ground state and 1 for the central vortex state) is presented and confirms that for continuous stirring (thin line), the vortex cycles in and out of the condensate.

In contrast to the two-dimensional case, the three-dimensional simulation can be unambiguously related to a realistic current experimental scenario. The $`C=1000`$ case shown in Figs. 13 corresponds to $`N=1.8\times 10^4`$ atoms of 87Rb in the $`|F=2,m_f=2`$ hyperfine state (for which we use an s-wave scattering length $`a=5.29`$ nm) in a time orbiting potential (TOP) trap with radial trap frequency $`\omega _r=2\pi \times 15`$ Hz and $`\lambda =\omega _z/\omega _r=\sqrt{8}`$. For this example, our computational units of time and length  correspond to $`10.6`$ ms and $`1.97`$ $`\mu `$m respectively, and the Thomas-Fermi diameter of the condensate in the radial direction of 11.7 corresponds to about 23 $`\mu `$m. The stirrer Gaussian $`1/e`$ diameter $`w_s=4`$ corresponds to a Gaussian $`1/e^2`$ beam diameter of approximately 11 $`\mu `$m. The required intensity of the stirring laser is detuning dependent, but can be easily calculated from the light shift potential $`\mathrm{}\mathrm{\Omega }^2/4\mathrm{\Delta }`$, where $`\mathrm{\Omega }`$ and $`\mathrm{\Delta }`$ are the Rabi frequency and detuning respectively of the atom-laser interaction. For example, if the laser beam is 50 nm blue-detuned from the 780 nm atomic transition, $`\mathrm{\Delta }=1.65\times 10^{14}`$ s-1, and a laser power of 300 nW is required for the 3 nK height stirrer we describe. The first angular momentum peak in Fig. 3 occurs at 0.26 s, which is much less than the lifetime of the condensate. The occurrence of the single vortex cycling behaviour in this experimentally accessible parameter regime suggests a method of preparing a central vortex state; stir the condensate as we have described, then remove the stirrer at the time when the angular momentum peak is expected. We illustrate the implementation of this method in Fig. 3 with the thick curve, for which the stirrer is linearly withdrawn between $`t=7\pi `$ and $`t=8\pi `$. Once the stirrer is removed, the vortex remains in the condensate, circling about and close to the centre, and $`L_z`$ is constant.

The vortex cycling behaviour illustrated by the example in Figs. 13 indicates that the two-state model described in our earlier work (for the two-dimensional case ) also has application in three-dimensions. Here we briefly outline the modifications required to extend the two-state model into three dimensions. We represent the condensate by an axisymmetric component $`\varphi _s`$, and a vortex-like component $`\varphi _ve^{i\theta }`$ which has a circulation about the $`z`$ axis, where $`\varphi _s`$ and $`\varphi _v`$ are real and nonnegative, so

$$\psi (𝐫,t)=a_s(t)\varphi _s(\rho ,z,n_v)+a_v(t)\varphi _v(\rho ,z,n_v)e^{i\theta },$$
(3)
with $`\rho =\sqrt{x^2+y^2}`$. We choose $`\varphi _s`$ and $`\varphi _v`$ to satisfy the following time-independent coupled equations,

$`\mu _s\varphi _s`$ $`=`$ $`[{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{d}{d\rho }}\left(\rho {\displaystyle \frac{d}{d\rho }}\right){\displaystyle \frac{d^2}{dz^2}}`$ (6)
$`+{\displaystyle \frac{\rho ^2+\lambda ^2z^2}{4}}+C(n_s\varphi _s^2+2n_v\varphi _v^2)]\varphi _s,`$
$`\mu _v\varphi _v`$ $`=`$ $`[{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{d}{d\rho }}\left(\rho {\displaystyle \frac{d}{d\rho }}\right){\displaystyle \frac{d^2}{dz^2}}+{\displaystyle \frac{1}{\rho ^2}}`$ (8)
$`+{\displaystyle \frac{\rho ^2+\lambda ^2z^2}{4}}+C(n_v\varphi _v^2+2n_s\varphi _s^2)]\varphi _v,`$
where the eigenvalues $`\mu _s`$ and $`\mu _v`$ are found by solving Eqs. (8) for a given choice of the vortex fraction $`n_v=|a_v|^2`$. The quantity $`n_s=|a_s|^2=1n_v`$. We note that when $`n_v=0`$, Eqs. (3) reduce to Eq. (6), the time-independent GP equation for the ground state, while if $`n_v=1`$, Eqs. (3) reduce to Eq. (8), the time-independent GP equation for the first excited vortex. The resulting nonlinear Rabi equations for the two-state system \[obtained by substituting Eq. (3) into Eq. (1)\] are then identical in form to Eqs. (5–6) of , but with the wavefunctions $`\varphi _s`$ and $`\varphi _v`$ now defined over $`\rho `$ and $`z`$ instead of $`r`$.

The underlying assumption of our two-state model is that the ground state (or more generally, $`\varphi _s`$) is coupled primarily to the vortex state of lowest energy with axis parallel to the stirrer axis. In two dimensions, there are a limited number of excited condensate states in the appropriate energy range that can be coupled to the symmetric state by the stirring potential, and we have previously shown  that the two state model works well within a specified validity range. In three-dimensional cylindrically symmetric traps, many more condensate eigenstates exist, and although symmetry considerations limit the possible coupling from the symmetric state, a two state model places more constraint on the description of the system than in the two-dimensional case. We do not aim therefore to provide a detailed representation of the stirring behaviour with our two state model, instead we use it primarily to provide conceptual understanding of the cycling behaviour. Nevertheless we also find it gives a quantitatively accurate prediction of the critical frequency of rotation $`\omega _c`$, the stirring frequency which causes a single vortex to cycle right to the centre of the condensate. We have examined the single vortex cycling behaviour using the full GP solution, for a variety of trap geometries including the oblate ($`\lambda =\sqrt{8}`$), spherical ($`\lambda =1`$), and prolate ($`\lambda =1/3`$) cases, with a condensate of $`C=1000`$, and a stirring potential of the same size and position as in Fig. 1. We find that the critical rotational frequency for single vortex cycling is accurately predicted by the two state model as the difference in energy between the ground and vortex states of the condensate (modified by the finite size of the stirrer, see Eq. (7) of ). As in , a sufficiently large Rabi frequency is required to distort the energy barrier between ground and vortex state enough to permit cycling. When the critical frequency of rotation is exceeded, additional vortices penetrate the condensate.

It is worth noting that as the $`z`$ component of the trapping potential become weaker (that is, as $`\lambda `$ becomes smaller), then in addition to the vortex state, the stirring potential increasingly excites the dipole (centre of mass) oscillation of the condensate. This can be understood in terms of the Ehrenfest theorem: for our choice of parameters the mean force $`W`$ of the stirrer increases as the condensate radial size decreases because the condensate is increasingly concentrated near the steepest part of the stirring potential. The centre of mass motion contributes angular momentum $`L_z_{COM}`$ to the condensate angular momentum, where

$$L_z_{COM}=xP_yyP_x.$$
(9)
We find that $`L_zL_z_{COM}`$ (which is a measure of the angular momentum caused by the presence of the vortex) for all three trap geometries cycles between 0 and 1 when $`\omega _f`$ is slightly less than $`\omega _c`$ (the single vortex cycling regime) and oscillates about 1 when $`\omega _f`$ just exceeds $`\omega _c`$ (and multiple vortices penetrate the condensate). $`L_z_{COM}`$ is small for the $`\lambda =\sqrt{8}`$ case shown in Fig. 3.

If the single stirrer used in Fig. 1 is replaced by two stirrers of half the height but on opposite sides of the $`z`$ axis, no vortex cycling occurs, and $`L_z`$ exhibits only small oscillations of amplitude less than $`0.04`$. This is easily understood from our two state model, because the Rabi frequency $`\mathrm{\Omega }`$ (see Eq. (6b) of ) is zero for this stirrer configuration, as it is for any stirrer configuration of even symmetry. It is worth pointing out that our scenario of coherent cycling should be distinguished from vortex production by damping in a trap with a symmetric rotating distortion, such as the recently reported experiment , or numerical damping schemes such as imaginary-time propagation .

Our three-dimensional simulations also allow us to test (in the mean-field limit) the stability of pure vortex eigenstates. Using as an initial state a central vortex (an eigenstate of the time-independent GP equation) for $`C=1000`$, a perturbation was applied by linearly inserting and withdrawing (between $`t=0`$ and $`t=\pi `$) a stirrer of the same size as in Fig. 1 but at a fixed position in the laboratory frame. We find that for a range of trap geometries ($`\lambda =\sqrt{8}`$, $`1`$, or $`1/3`$), the $`L_z=1`$ vortex is dynamically stable, and remains near the centre of the condensate. We can understand this result in terms of Rabi cycling. The stationary stirrer ($`\omega _f=0`$) provides a potential which is far-detuned from resonance ($`\omega _f=\omega _c`$), and thus there is a low probability of transfer out of the initial (vortex) state. Indeed if the stirrer remains fixed in the condensate after insertion and is not withdrawn, we find that the vortex still remains near the centre of the condensate and follows a small closed path displaced slightly away from the stirrer. By contrast, for these geometries, the $`L_z=2`$ central vortex eigenstate immediately dissociates into two unit vortices when perturbed. We note that these stability results refer to dynamical stability in the GP limit (that is, $`T=0`$), as opposed to thermodynamic stability .

In addition to confirming the basic features of the two state model of vortex cycling, our three-dimensional simulations also reveal some new phenomena. The vortices produced in three dimensions by stirring are line vortices. In the $`\lambda =\sqrt{8}`$ configuration, the lines are nearly straight and parallel to the stirrer axis, and when two lines of the same circulation form, we find they repel each other. For a spherical trap ($`\lambda =1`$) however, vortex lines along any direction are degenerate (in the absence of the stirrer) with the vortex eigenstate about the $`z`$ axis. Accordingly, the vortex lines can become more curved than in the oblate case, as we find in our simulations below the critical frequency, where the vortex that forms can have appreciable curvature. Above the critical frequency, when multiple vortices are drawn into the condensate, there is a striking consequence of this potential for vortex lines to deform. In a collision between two such curved vortex lines, sections of the lines can be exchanged, as we show in the collision sequence in Fig. 4. The condensate is stirred somewhat above the critical frequency in a spherical trap, and by $`t=6.79`$ two vortex lines (both of positive sense ) have formed. The figure shows only the vortex lines (detected numerically), which have been shaded light and dark to distinguish one from the other. In Fig. 4(a) the dark vortex line approaches the light vortex line from behind, and because of their mutual repulsion, the light line has developed a bulge. In Fig. 4(b) the dark line is now close to the light line, causing the bulge to become recurved, with two small portions now nearly antiparallel to the dark line. In Fig. 4(c) the two lines have completed a collision and exchanged their central sections.

In conclusion, we have shown that the single vortex cycling regime we discovered previously in simulations of two-dimensional condensates is also present in three-dimensional condensates, for a range of trap geometries. We have related our simulations to current experimental configurations, and have demonstrated how the cycling behaviour could be used to generate a central vortex state. We have investigated the dynamical stability of central vortex states, and have also given initial evidence of the rich dynamical behaviour of vortex lines.

This work was supported by the Marsden Fund of New Zealand under contract PVT902.


text:
# Extended X-ray emission in the radio loud galaxy 3C382

## 1 Introduction

Soft X-ray excess emission above a simple extrapolation of the hard energy spectrum is found in a considerable number of AGNs, mainly radio-quiet sources (see Mushotzky, Done, & Pounds 1993 and references therein). There is a growing body of evidence for its spectral ubiquity below $``$2 keV or so in Seyferts (e.g. Pounds et al. 1994) as in highly luminous quasars (e.g. Saxton et al. 1993). Most previous studies before ROSAT converge toward the idea that the soft excess is a rather common feature among radio-quiet quasars, whereas it is almost absent in their radio-loud counterparts. ROSAT/PSPC data of radio loud sources have shown that a soft excess component is also present in radio loud sources (e.g. Buehler et al. 1995; Prieto 1996; Siebert et al. 1998).

Possible interpretations for the soft X-ray excess in AGNs include thermal emission from the inner regions of an accretion disk, scattering by highly ionized material in its vicinity (Pounds et al. 1986; Ross & Fabian 1993), or thermal emission due to shock-heated gas in the close vicinity of the nucleus (Viegas & Contini 1994). The poor spatial resolution of the ROSAT/PSPC makes difficult the separation between possible components of the observed emission. Indeed, the large PSPC resolution beam, $``$ 25 arcsec at 1 keV, makes plausible that an important part of the observed emission to be due to an extended gas component surrounding the AGN. In the particular case of radio-loud galaxies which are characterized by large radio sizes, a hot surrounding medium becomes a necessary component for providing the working surface for the radio emission. In the analysis of the 3CRR sample by Prieto, a first attempt to fit the PSPC spectra of sources with extended emission –mostly in Fanaroff & Riley (1974) type I sources (FRI)– with a single power-law leaded to extreme step spectral index, the reason being due to the dominant contribution of the gaseous medium in which those sources usually reside. In the case of FRII, a single power-law fit provided a fair representation of the PSPC spectrum but with average spectral index about -1.1, and so above the extrapolation of the canonical hard-energy spectrum into the soft X-rays. Clustering of galaxies about FRII sources is less common than in FRI, in particular at low redshift; yet, FRII could contain their own extended gaseous atmosphere which may directly translate into a steepening of the PSPC spectrum. This component however may prove to be elusive with present X-ray instrumentation.

This paper presents deep ROSAT/HRI observations of the powerful X-ray radio-loud source 3C 382. This is one of the few nearby broad-line galaxies ($`z=`$0.0578) that show extremely bright and broad permitted lines (FWZI$`>25000kms1`$; Tadhunter, Perez & Fosbury 1986) and a strong continuum, with with a X-ray luminosity in 0.2-2.4 keV band of $`L_\mathrm{x}7.10^{44}`$ erg s-1 (Prieto 1996), and a radio power at 178 MHz of $`L_{178\mathrm{M}\mathrm{H}\mathrm{z}}3.10^{33}`$ erg s$`{}_{}{}^{1}Hz^1`$ (Laing et al. 1983).

EXOSAT monitoring of the source (1983–1985) tightly constrains the high-energy (above 2 keV) spectral index of 3C382 to $`\alpha =0.7\pm 0.1`$ (Ghosh & Soundararajaperumal 1992). This is also confirmed by more recent ASCA data (Wozniak et al. 1998). However, the ROSAT/PSPC spectral analysis of the source shows compatible with a power-law model with spectral index $`\alpha =1.2\pm 0.3`$, and absorbed by a column density, $`N(H)=0.78\times 10^{21}cm^2`$, that is in agreement with the Galactic value. Thus, 3C382 shows a soft excess emission below $``$ 2keV (Prieto 1996). Independently, the presence of a soft excess is also inferred from the analysis of the EXOSAT (Ghosh & Soundararajaperumal), ASCA (Wozniak et al.) and Ginga (Kaastra et al 1991) data.

Extended soft X-ray emission associated with this source is detected in the ROSAT/HRI data. On the basis of that new component a re-evaluation of the PSPC spectrum is presented.

Throughout this work $`H_0`$ = 50 km s-1 Mpc-1. 1 arcsecond corresponds to $``$ 1.7 kpc at the source.

## 2 Analysis of HRI data

The HRI observations of 3C 382 were conducted in 1996 October and 1997 April (WG900720H and WG900720H-1 datasets respectively) The corresponding total accepted times were 4514 s and 13310 s, respectively. The counts are integrated from channels 1 to 8 which enclose most of the energy accumulated in the ROSAT band, yielding for both dataset count rates of $``$1 cts s-1.

The nominal resolution of the ROSAT/HRI is $``$ 5–6 arcsec (FWHM). Residual errors in the ROSAT aspect solution are known to give rise to elongated images (David et al. 1996), the shape of the surface brightness profile dramatically departing from the expected point response function (PRF). In the case of very bright sources, improvement of the HRI spatial resolution becomes feasible by using speckle interferometric techniques such as the “shift-and-add” method. If one constructs images over short time intervals, the PRF becomes symmetric and therefore the elongation in the image appears as an apparent residual motion of the X-ray source in the sky. Such residual motion can be corrected for by de-speckling.

3C 382 is particularly suitable for that technique as it shows very bright in the ROSAT band, with $``$ 2 cts s-1 in the PSPC and $`1`$ cts s-1 in the HRI.

The procedure used follows the same criteria and approach originally presented in Schmitt, Güdel & Predehl (1994). Basically, each event file is divided into time bins of typically 50 s. For each bin, the apparent X-ray position of the source (i.e., RA and $`\delta `$) is determined as a function of the observing time. A spline function is fitted to these data points, all recorded photons being then corrected with the appropriate time-dependent correction in RA and $`\delta `$.

To validate the correction, new measurements of the source centroid are repeated on the corrected event file. The correction is considered as satisfactory if the new centroid positions cluster during the time period of the observation about an average constant value. The uncertainty in the final source centroid is about 2.5 arcsec in RA and 1.5 arcsec in $`\delta `$.

The X-ray spatial analysis was then performed on the corrected HRI event files. Because of the much higher statistic of the April event file, reliable results from the de-speckle procedure are only found from that dataset. Thus, the following analysis focuses only on this dataset.

## 3 The extended X-ray component

An X-ray contour image of 3C 382 obtained after the de-speckle procedure is shown in Figure 1a. The image is background-subtracted and smoothed with a Gaussian filter of FWHM $``$12 arcsec. The background level is estimated from different regions in the image within the central 5 arcmin. The emission is dominated by a central peak component, and a surrounding halo slightly more asymmetric towards the North-East side. Some of the morphological features are also apparent in the shortest exposed 1996 October event file -not shown. Overall, the X-ray emission extends out to about 100 arcsec from the center, $``$170 kpc at the distance of the source. This is about the size of the radio structure at 8.3 GHz (Black et al. 1992). The slight North-East asymmetry in the X-ray image is virtually coinciding with the direction of the radio axis and with that of the faint filamentary regions and bright companion galaxy seen in the HST/WFPC2 image (Fig. 1b).

The corresponding surface-brightness profile is presented in Figures 2. For comparison, the new improved ROSAT/HRI PSF profile by Predhel (1998) —currently upgrade in the EXSAS package— which includes the effect of the ROSAT mirror scattering and has been fit to Capella and Sirius, is shown superimposed. To illustrate the effect of the de-speckle procedure on the data, the surface brightness profile as derived from the original data is shown is Fig 2a, and that after de-speckling in Fig. 2b. The departure from the ROSAT/HRI PRF is clearly evidenced in Fig. 2a, where larger residuals are seen all over the profile. The de-speckle procedure produced a sharper profile which virtually fills the core of the PRF (Fig. 2b). Yet, some systematic residuals still remain at radius beyond $``$10–15 arcsec from the center, which we assume to be related to an extended emission component. This component is roughly 10% of the total emission.

In an attempt to get a better characterization of the observed profile, a combination of an unresolved component represented by the HRI/PSF, and an extended component represented by a $`\beta `$-model are fit to the data. The functional form of the $`\beta `$-model follows King’s approximation (1972) for gas confined in an isothermal sphere: $`S(r)(1+(r/r_c)^2)^{3\beta +1/2}`$, with S(r) the surface brightness a radius r, $`r_c`$ the core radius and $`\beta `$ the slope parameter. Given the uncertainty inherent with this approximation for the extended component, a simple addition of the PSF and the beta model is fit to the data.

Figure 2c shows the resulting fit corresponding to the composite model: PRF plus a $`\beta `$ model. Fig 2c shows a clear fit improvement to the observed profile as the systematic residual trend seen in Fig 2b is now removed; yet, the composite model still overpredicts the observed emission at radius beyond 10 arcsec. Overall, the $`\beta `$-model does not provide a statistically acceptable fit with reduced $`\chi ^2`$ exceeding unity. The main reason for that may reside in the validity of the beta-model for the specific case of the extended gas emission in 3C382. Besides, there is the fact that the dominant contribution of the emission is the unresolved component, which makes the modeling of the extended diffuse emission difficult. Also, there is the rather asymmetric morphology of the emission which may clearly depart from the simple model used here.

Nevertheless, the fit results from the composite model leaded to a tight range for both the core radius and the $`\beta `$ value. The minimum reduced $`\chi ^2`$, $`<1.52>`$, was obtained for a core radius $`r_c`$ 20-30 arcsec (30 - 50 kpc at the distance of 3C 382) and $`\beta 0.70.8`$. A much simpler composite model represented by the PRF plus a Gaussian model for the extended component, yielded FWHM for the extended component in the 30 - 40 arcsec range.

### 3.1 Re-evaluating the PSPC spectrum

A single power-law model with spectral index -1.2 provides a fair fit of the PSPC spectrum of 3C382, with the derived N(H) in good agreement with the Galactic value. However, the detection of an extended component in the HRI data along with the reported evidence for soft excess emission in this source from EXOSAT and ASCA data besides PSPC (cf. sect. 1) prompt to a re-evaluation of the PSPC data. Because the extended nature of the new component, its association with thermal emission arises as the most natural explanation.

Due to the short energy range and relatively low spectral resolution of the PSPC, a complex analysis of the PSPC spectrum involving a combination of several models is difficult. However, the evidence for extended emission and the fact that the hard X-ray spectral index of this source remains constant about the 0.7 value are additional inputs that can be used for forcing the PSPC fit with a more complex modeling.

Accordingly, a combination of a power-law, with spectral index fixed at $`\alpha `$= –0.7, and a thermal model are fitted to the PSPC data. In this composite model, $`N`$(H) is fixed at the Galactic value following Prieto’s results (1996). Thus, the only parameters that vary freely are the normalizations of the respective thermal and power-law components and the temperature of the thermal component.

The new composite fit (Fig. 3) provides a fair representation of the PSPC spectrum, with a reduced $`\chi ^2`$1.2 and constrained values for all the free parameters. A gas temperature of $`0.6_{0.1}^{+0.4}`$ keV is derived. Alternative fits, letting free the index of the power-law spectrum or the N(H), lead to larger $`\chi ^2`$ and unconstrained fit parameters. Fitting a single bremsstrahlung model produces also an acceptable fit but the derived N(H) value is found lower than the Galactic value.

The results from the composite model can be compared with those derived from the single power-law model and from the single thermal model in Table 1. Given the PSPC spectral limitations together with the additional complexity of the composite model, the derived fit values may be subjected to larger uncertainties. In this sense, the temperature could be largely affected whereas integrated fluxes are the least dependent on the adopted model. The total flux contribution from the thermal and power-law components in the composite model compares with that derived from the single power-law model. The reduced $`\chi ^2`$ is slighter better in the case of a simple power-law model; but we consider the difference marginal. The lowest reduced $`\chi ^2`$ is obtained with the thermal model; yet the derived N(H) in this case is inferior than the Galactic value.

### 3.2 How do the PSPC fluxes compare with those derived from the HRI data?

Fluxes estimate from the HRI data for the extended and unresolved component are primarily derived on the basis of HRI spatial analysis (§3). The total number of counts within the unresolved and extended component leads to count rates of $``$ 1 cts s-1 and about 0.1 cts s-1 respectively.

The lack of spectral resolution of the HRI hampers any direct spectral modeling of the X-ray data. Yet, the availability of the PSPC data allows us to use the same model and fit parameters as derived from the PSPC fit (§3.1) to convert the HRI counts to fluxes. These are given in Table 1

The unresolved HRI component, modeled with a power-law index $`\alpha `$= –0.7, yields a flux about a factor 2 larger than that measured by the PSPC. 3C 382 is known to be a variable source, with reported maximum-to-minimum variations of up to 120% as measured by EXOSAT (Ghosh & Soundararajaperumal 1992). Regarding the ROSAT observations, previous HRI observations of the source in 1992 March revealed a drop in the total number of counts by a factor 1.8 with respect to the values measured in the present observations taken 5 years later. This drop is thus compatible with the still lower flux measured by the PSPC in 1990. Thus, the difference found between the HRI and PSPC fluxes for the unresolved component appears compatible with the observed variability level of the source.

Regarding the extended component, bremsstrahlung models with temperatures within the range 0.6-3 keV lead to corresponding HRI fluxes $`106.\times 10^{12}`$ erg cm$`{}_{}{}^{2}s_{}^{1}`$ respectively, i.e., a factor 2 to 3.5 smaller than that derived from the PSPC for the same component. We note however that as neither the Gaussian approximation nor the $`\beta `$-model provide a statically acceptable fit to that component, the corresponding number of HRI counts for the extended component was not derived from the model but from the difference between the total integrated number of counts and that integrated within the unresolved component represented by the PRF. Given the simplicity of the method, in particular taking into account that the dominant contribution of the total emission is the unresolved component, the derived PSPC and HRI fluxes can be considered consistent with each other within the order of magnitude. Besides, it is also possible that the emission is more extended than what is detected –at least beyond the radio structure–, but the surface brightness there is too low for being detected with the HRI.

To summarize, the inferred luminosities for the extended and unresolved components estimated from the HRI data appear compatible in order of magnitude with the respective thermal and power-law luminosities derived from PSPC spectrum of the source.

## 4 Discussion

Extended X-ray emission associated with the very bright, nearby radio-loud galaxy 3C 382 is detected. The analysis of the ROSAT/HRI data shows that about 10% of the total 0.2 -2.4 keV emission is compatible with the presence of an extended component. Assuming that component to be due to hot gas emitting via bremsstrahlung, and allowing for the contribution of such thermal component into the PSPC spectrum, it is found that the non-thermal component of 3C 382 emission becomes consistent with the extrapolation of the “well established” 3C 382 high-energy power-law spectrum —above 2 keV— into the soft X-ray regime.

The temperature of the gas component as formally derived from the PSPC fit is $`0.6_{0.1}^{+0.4}`$ keV. This low temperature contrasts with the high luminosity of the gas component, $`3\times 10^{44}ergsec^1`$. Taken together the spectral limitation of the PSPC and the complexity of the model fit, the uncertainty in the temperature could however be larger. There is furthermore the possibility that a temperature gradient dominates the gas emission — a cooling flow process. In this case, the PSPC spectrum may be mostly sampling the central gas region, where the coolest and more dense gas is located. There are however other factors that could also have lead to an overestimation of the gas luminosity. The analysis by Markevitch (1998) on clusters with strong cooling flows indicates moderate temperature increase of up to 20% but luminosity decrease of up to 40% for the cluster gas after excising the cooling flow regions. On the other hand, if cooling by metals is considered – for sake of simplicity, pure bremsstrahlung has been assumed – the luminosity of the gas could decrease by about 40%, assuming a Raymond-Smith model with metal abundance Z=0.35.

Still, the derived thermal luminosity in 3C382 is about two order of magnitude larger than that found in isolated, normal elliptical galaxies (Canizares, Fabbiano & Trincheri 1987), and in low-power radio galaxies (Worral & Birkinshaw, 1994) but it is in the range found in powerful radio sources (Worrall et al. 1994; O’Dea et al. 1996; Hardcastle et al. 1999; Crawford et al. 1999). Also, the estimated core radius for 3C382, $``$50 kpc, is within the range found by Crawford et al. and Hardcastle et al. in their respective samples of 3CRR radio-loud sources.

Large X-ray halos are often seen in FRI sources, those being associated with the cluster environment in which they often reside (e.g. M 87, Perseus, 3C465). These halos largely dominate the ROSAT emission from these sources. Besides the outstanding case of Cygnus A, evidence for clustering is less obvious in classical double radio sources, particularly at low redshift (cf. Hill & Lilly 1991; Miller et al. 1999). Unambiguous extended X-ray emission in powerful FRII radio galaxies and quasars has mostly being found in sources with redshift larger than 0.1 (Hardcastle and Worall, 1999; Crawford et al 1999; O’Dea et al. 1996); yet, a few low-redshift FRII are reported to show extended X-ray emission (cf. Hardcastle and Worral). In most of these cases, the large X-ray luminosities are found compatible with thermal emission from a moderately rich cluster environment.

Comparing with Cyg A, the archetypal double radio source at z= 0.0574, 3C382 is also one of the few very bright double sources at low redshift with extended X-ray emission. Contrarily to Cyg A which presents an optical narrow line spectrum, 3C382 presents an extreme, in width and strength, broad permitted line spectrum. If this difference is interpreted as due to obscuration of the AGN region in Cyg A, it may explain why the dominant X-ray feature in Cyg A is emission from a hot diffuse gas –the AGN component is obscured at X-ray waves– whereas in 3C382, the unresolved X-ray nuclear component –presumably associated with the AGN– dominates the total X-ray emission, making more difficult the detection of any extended gas component.

The X-ray luminosity of 3C382, of about $`10^{44}ergs^1`$, compares with that of rich Abell clusters (this is also the case of Cyg A; yet Cyg A is at the center of a poor cluster of galaxies). Longair & Seldner (1979) derived however a rather poor environment in the vicinity of 3C 382 on the basis of their cross-correlation analysis between the radio position and galaxy counts. HST/WFPC2 images of 3C382 collected in parallel mode show an elliptical galaxy with a very bright unresolved nucleus and a halo very smooth (Martel et al. 1999). Yet, within the 2.5 arcminutes field of view (Fig 1b), several small galaxies can easily be distinguished in the 300 seconds exposure; the WFPC2 images also show a bright galaxy at 85 arcsec Northeast from 3C382, presenting two at least extended gaseosus tails of material in the direction of 3C382; two additional difuse regions located close to 3C382 and in the direction of the bright galaxy are also apparent. Judging from the HST images, 3C382 may be residing in a relatively poor cluster environment; also, it may be in interaction with that gas-rich galaxy companion. Such interaction could have brought plenty of gas into 3C382.

The luminosity of the halo component in 3C382 would imply a large mass of gas, of about $`10^{11}`$Mo, assuming it concentrated in a sphere of about 50 kpc –the estimate core radius derived from the HRI spatial analysis– and a temperature in the 0.6 -1 keV range. An alternative to the cluster environment is that 3C382 it may consist of a self-contained gravitational potential deep enough to restrain such large amount of gas. This could also be the case of the radio-louds 3C48 and 3C273, for which extended X-ray emission is found but the evidence for a cluster environment from optical images is minor (Crawford et al. 1999). Evidence for a massive dark halo in 3C382 comes from the velocity measurements on the extended ionized gas surrounding this galaxy. Tadhunter et al. (1986) detected ionized gas up to 25 kpc from the galaxy center. The gas follow a a rotation curve which extend flat up to those distances with velocities of about 400 km/s relative to the systemic velocity. Assuming a spherical potential, the estimated mass within a 25 kpc radius would be $`8\times 10^{11}Mo`$. This is about the gravitating mass needed to keep the X-ray gas binded to the galaxy. Following Fornan, Jones and Tucker (1985) formalism, the total gravitating mass within a 50 kpc radius is estimated between $`815\times 10^{11}Mo`$ for gas temperatures between 0.6 and 1 keV, which is in the order of magnitude of the mass derived from the ionized gas kinematics.

The results so far derived show compatible with a cooling flow process being dominating the gas emission. If the extended gas emission is modeled as that of a uniform sphere of hot gas emitting via bremsstrahlung, for a maximum radius of about 170 kpc (the size of the radio structure) and a temperature in the 0.6 - 1 keV range, the implied density would be $`simeq6\times 10^3cm^3`$. This yields a cooling time of about $`4\times 10^9yr`$, considerable smaller than the Hubble time. Thus, a cooling flow process could be operating in 3C 382.

A better characterization of the extended X-ray emission in 3C382 would demand much larger spatial resolution but also deeper observations.

Acknowledgments:

It is a pleasure to thank Peter Predehl for his support regarding the de-speckling procedure of the HRI data and Günther Hasinger and Paddy Leahy for critical reading of early versions of the manuscript.

## References

Black, A. et al. 1992, MNRAS 256, 186

Buehler, P., Courvoisier, T., Staubert, R., Brunner, H., & Lamer, G. 1995, A&A, 295, 309

David, L.P. et al. 1998, The ROSAT High Resolution Imager (HRI) Calibration Report: $`http://`$hea-www.harvard.edu./rosat/rsdc\_www/HRI\_CAL\_REPORT/hri.html

Arnaud, K. A. 1993, ApJ 412, 479

Canizares, C. R., Fabbiano, G. & Trincheri. G 1987, ApJ 312. 503

Carrili, C. & Barthel, P. 1996, A&ARv 7,1

Crawford, C.S. & Fabian, A. 1995, MNRAS 273, 827

Crawford, C.S., Lehmann, I., Fabian, A., Bremer, M.N. & Hasinger, G. 1999, MNRAS 308, 1159

Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31p

Ghosh, K. K., & Soundararajaperumal, S. 1992, ApJ, 389, 179

Hardcastle, M. J. & Worral, D. 1999, MNRAS 309, 969

Hill G. J., & Lilly, S. J. 1991, ApJ, 367, 1

Kaastra J.S, Kunieda & Awaki H. 1991, A&A 242, 27

King, I. 1972, ApJ 174, L123

Laing, R. A., Riley, J. M., & Longair, M. S. 1983, MNRAS, 204, 151

Longair, M.S. & Seldner, M. 1979, MNRAS 189, 433

Markewitch, M., 1998, ApJ 504, 27

Martel et al., 1999, ApJSS 122, 81

Miller, N.A. et al, 1999, AJ 118, 1988

Mushotzky, R. F., Done, C., & Pounds, K. A. 1993, ARA&A, 31, 717

O’Dea, C. P., Worrall, D., Baum, S., & Stanghellini, C. 1996, AJ, 111, 92

Pounds K.A., Nandra, K., Fink, H., Makino, F. 1994, MNRAS 267, 193

Pounds, K. A., Turner, T. J., & Warwick, R. S. 1986, MNRAS, 221, 7p

Predehl, P. 1998, in preparation

Prieto, M. A. 1996, MNRAS, 282, 421

Ross, R. R., & Fabian, A. C. 1993, MNRAS, 261, 74

Saxton, R. D., Turner, M. J. L., Williams, O. R., Stewart, G. C., Ohashi, T., & Kii, T. 1993, MNRAS, 262, 63

Siebert, J. et al. 1998, MNRAS in press

Schmitt, J.H.M.M., Güdel, M. & Predehl, P. 1994, A&A 287, 843

Tadhunter, C., Perez E. & Fosbury, R. 1986, MNRAS 219,555

Viegas, S. M., & Contini, M. 1994, ApJ, 428, 113

Worral, D. & Birkinshaw, M. 1994, ApJ 427, 134

Worral, D., Lawrence, C. R., Pearson, T., & Readhead, A. C. S. 1994, ApJ, 420, L17

Wozniak P.R. et al. 1998, MNRAS 299, 449

Figure captions

Figure 1a:
HRI contour image of 3C 382 extracted from channels 1 to 8. It is background-subtracted and smoothed with a Gaussian filter with a FWHM $``$12 arcsec. Contours are $`10^3ctss^1arcmin^2\times `$ (2.8, 4.4, 6, 7.6, 9.2, 10, 16, 40,4000); the first contour is about the 2$`\sigma `$ level measured on the background-subtracted image. 1 sky-pixel is 0.5 arcsec.

Figure 1b:
Broad band HST/WFPC2 image of 3C382 (data set u27l6w01-2) at 7200 A an equivalent exposure of 300 seconds. 3C382 is filling the upper left CCD. The bright companion galaxy (upper right CCD) is North-East of 3C382. The image covers 2.5 arcminutes.

Figure 2a:
HRI surface brightness profile of 3C 382 before de-speckling. Points with error bars (Poissonian noise) are the data. Dashed line is the HRI PRF; the continuum line is the PRF plus background. The background level is measured in an annulus between 3 to 5 arcmin from the center (flat part of the profile). The residuals between the data and the PRF+background are shown below. Note the large departure from a point source profile due to incorrect attitude reconstruction of the ROSAT data (see text).

Figure 2b:
HRI surface brightness profile of 3C 382 after de-speckling. The interpretation of the figure is as in Fig. 2a.

Figure 2c: The same as in Fig 2b but in this case the continuum line is the model fit to the surface brightness profile using a combination of an unresolved component represented by the HRI PRF (dash) and an extended component represented by a $`\beta `$ model (dots). The reduced $`\chi ^2`$ of the fit is 1.7-2 over 35 degrees of freedom (d.o.f.).

Figure 3:
Observed and best-model fit to the PSPC spectrum of 3C 382. PSPC data are from Prieto (1996). The dotted line represents the thermal 0.6-keV component; the dashed line, the power-law $`\alpha =0.7`$ component. Residuals represent 1-$`\sigma `$ errors. The first two bins corresponding to energies $`0.1keV`$ are not included in the fit.