text:
# Giant diffusion and coherent transport in tilted periodic inhomogeneous systems.

## I Introduction

In recent times there has been a renewed interest in the study of transport properties of Brownian particles moving in periodic potential Ris , with special emphasis on coherent transport and giant diffusion Main ; March . The phenomenon of coherent or optimal transport is complimentary to the enhanced diffusion, wherein one is mainly concerned with transport currents with minimal dispersion or diffusion Fnl . Compared to free diffusion coefficient(DC, $`D=k_BT/\gamma `$), DC is suppressed in the presence of periodic potential. However, in a nonequilibrium case i.e., in the presence of bias, it has been recently shown that DC can be made arbitrarily large (giant diffusion) compared to the bare diffusion, in the presence of periodic potential Main . This enhancement at low temperature takes place near the critical threshold ( at which deterministically running solution sets in). The reason for this enhancement has been attributed to the existence of instability between locked to running solution. In some cases enhancement by fourteen order of magnitude has been predicted, so that diffusion can be observed on a macroscopic scale at room temperature Main . This enhancement decreases as we move away from the critical field in either direction. Exact result for DC in arbitrary potentials has been obtained in term of quadratures. In special cases an elegant simplification of quadrature have been carried out. Near the critical tilt, scaling behavior of DC for weak thermal noise has been obtained and different universality classes have been identified Main . Approximate expression for DC in terms of mobility has been obtained earlier which deviates from the exact results near the critical threshold March .

In a related development, study of coherent or optimal transport has been reported recently Fnl . Coherent or optimal transport of Brownian particles refer to the case of large mean velocity accompanied with minimal diffusion. This can be quantified by dimensionless Péclet number (ratio of mean velocity to DC). The transport is most coherent when this number is maximum. The particle motion is mainly determined by two time scales; noise driven escape from potential minima over the barrier along the bias, followed by the relaxation into the next minima. The former depends strongly on temperature and the later weakly on the noise strength and has a small variance. It is possible to obtain coherent transport in the parameter regime at which the traversal time across the two consecutive minima in a washboard potential is dominated by the relaxational time. At optimal noise intensity certain regularity of the particle motion is expected which accounts for the maximum of Péclet number. In some cases (molecular separation devices) for higher reliability, one requires higher currents but with less dispersion (or DC) Geier . This effect of coherent transport is related to the phenomenon of coherence resonance Cr  observed in excitable systems Fnl .

In the present work we study both the mentioned phenomena in a space dependent frictional medium. For this we have considered a simple minimal model where the potential is sinusoidal and the friction coefficient is also periodic (sinusoidal) with the same period, but shifted in phase. Frictional inhomogeneities are not uncommon in nature. Here we mention a few. Brownian motion in confined geometries, porous media experience space dependent friction Conf . Particles diffusing close to surface have space dependent friction coefficient Conf ; Surf . It is believed that molecular motor proteins move close along the periodic structure of microtubules and will therefore experience a position dependent friction Spc\_mob . Inhomogeneities in mobility ( or friction) occurs naturally in super lattice structures and Josephson junctions Falco .

Frictional inhomogeneity changes the dynamics of the diffusing particle non-trivially, thereby affecting the passage times in different regions of the potential. However, it does not effect the equilibrium distribution. Thus thermally activated escape rates and relaxational rates within a given spatial period are affected significantly. This in turn has been shown to give rise to several counter-intuitive phenomena. Some of them are stochastic resonance in absence of external periodic drive Sr , noise induced stability in washboard potential Dan . Single and multiple current reversals in adiabatic Danphys  and nonadiabatic rocked system (thermal ratchets Ratchet ) respectively have also been reported Dan2 . In these ratchet systems the magnitude of efficiency of energy transduction in finite frequency regime may be more than the efficiency in the adiabatic regime, i.e, quasistatic operation may not be efficient for conversion of input energy into mechanical work Dan3 . All these above features are absent in the corresponding homogeneous case for the same simple potential.

In our present work we show that frictional inhomogeneities can give rise to additional new features in a tilted periodic potential. The observed giant enhancement of DC near the critical tilt can be controlled (enhanced or suppressed) in a space dependent frictional medium by suitably choosing the phase difference between the potential and the frictional profile. The most surprising feature is the noise induced suppression of diffusion, leading to coherent transport. Our results are based on analytical expressions for Péclet number and DC in term of moments of first passage times.

In Section II we present our model and derive the expression for moments of first passage times. Using these, DC and Péclet number have been defined. In section III.1 we analyze the nature of DC as function of system parameters, such as the applied external force and temperature. Section III.2 is devoted to the study of coherent or optimal transport in different regimes of parameter space. Finally in section IV we present the summary of our findings.

## II Model

We consider an overdamped Brownian particle moving in a symmetric one dimensional periodic potential $`V_0(x)`$ with space dependent friction coefficient $`\eta (x)`$ under the influence of constant external tilt $`F`$ at temperature $`T`$. For simplicity we take $`V_0(x)=\mathrm{sin}(x)\text{ and }\eta (x)=\eta _0(1\lambda \mathrm{sin}(x+2\pi \varphi )),\text{ where }|\lambda |<1`$. $`\varphi `$ determines the relative phase shift between friction coefficient and potential. The correct Langevin equation for such systems has been derived earlier from microscopic considerations Langvn  and is given by

$$\dot{x}=\frac{(V_0^{}(x)F)}{\eta (x)}k_BT\frac{\eta ^{}(x)}{(\eta (x))^2}+\sqrt{\frac{k_BT}{\eta (x)}}\xi (t),$$
(1)
where $`\xi (t)`$ is a zero mean Gaussian white noise with correlation $`\xi (t)\xi (t^{})=2\delta (tt^{})`$. It should be noted that the above equation involves a multiplicative noise with an additional temperature dependent drift term which turns out to be essential for the system to approach correct thermal equilibrium state in absence of nonequilibrium forces Sancho . The quantity of central interest in this work is the effective diffusion coefficient $`D`$ given by

$$D=\underset{t>\mathrm{}}{lim}\frac{x^2(t)x(t)^2}{2t}.$$
(2)
In the absence of potential, $`D=\frac{k_BT}{\eta }`$ (uniform $`\eta `$), is the usual Einstein relation. An exact analytical expression for $`D`$ and average current $`J`$ in terms of the moments of first passage time have been recently given Main ; Fnl . If the n-th moment of the first passage time from an arbitrary point $`x_0\text{ to }b`$ is given by $`T_n(x_0b)=t^n(x_0b)`$, then

$$D=\frac{L^2}{2}\frac{T_2(x_0x_0+L)[T_1(x_0x_0+L)]^2}{[T_1(x_0x_0+L)]^3}$$
(3)
For our problem (1), the moments of first passage time follow closed recursion relation Main ; Gard .

$$T_n(x_0b)=\frac{n}{D_0}_{x_0}^b𝑑x\widehat{\eta }(x)e^{\frac{V(x)}{k_BT}}_{\mathrm{}}^x𝑑ye^{\frac{V(y)}{k_BT}}T_{n1}(yb),$$
(4)
with $`T_0(yb)=1`$. Here $`V(x)=V_0(x)Fx,D_0=k_BT/\eta _0\text{ and }\widehat{\eta }=\eta (x)/\eta _0`$. $`\eta _0`$ also happens to be the average value of friction coefficient over a period.

By using Eq. (3) and Eq. (4), and some straight forward algebra, we get

$$D=D_0\frac{_{x_0}^{x_0+L}𝑑xI_\pm (x)I_+(x)I_{}(x)}{[_{x_0}^{x_0+L}𝑑xI_\pm (x)]^3},$$
(5)
where,

$$I_+(x)=\frac{1}{D_0}\widehat{\eta }(x)e^{\frac{V(x)}{k_BT}}_{xL}^x𝑑ye^{\frac{V(y)}{k_BT}},$$
(6a)
$$I_{}(x)=\frac{1}{D_0}e^{\frac{V(x)}{k_BT}}_x^{x+L}𝑑y\widehat{\eta }(y)e^{\frac{V(y)}{k_BT}}.$$
(6b)
The average current density $`J`$ for this system has been derived earlier Dan  which in term of Eqs. (6) is given by

$$J=L\frac{1\mathrm{exp}(LF/kT)}{_{x_0}^{x_0+L}\frac{dx}{L}I_\pm (x)}.$$
(7)
The above expressions go over to the results obtained earlier for the case of space independent friction ($`\eta (x)=\eta _0,\lambda =0`$) Main . It should be noted that Eqs. (6) are applicable when $`\eta (x)\text{ and }V(x)`$ have the same periodicity. Otherwise they have to be modified appropriately. We obtain results for DC by numerically integrating Eqs. (6) using a globally adaptive scheme based on Gauss-Kronrod rules. For the special case of $`V_0(x)=0`$ we get

$$D=D_0(1+\frac{\lambda ^2}{2}\frac{1}{[(\frac{k_BT}{F})^2+1]}).$$
(8)
We would like to mention here that in the absence of potential, DC explicitly depends on system inhomogeneities ( via $`\lambda `$), however, steady current is independent of $`\lambda `$ for the same case Dan . $`F=0`$ is the equilibrium situation and as expected $`D=D_0`$, which corroborates with the fact that frictional inhomogeneities cannot affect the equilibrium properties of the system. In the high temperature regime, $`D=D_0`$ as anticipated. For asymptotically large field, DC saturates to a $`\lambda `$ dependent value. This is solely attributed to space dependent friction. This somewhat surprising result also appears in the dependence of current on $`\lambda `$ in the presence of potential and high field limit Danphys ; Dan2 .

In our subsequent analysis all the physical quantities are expressed in term of dimensionless units, DC is scaled with respect to $`D_0`$ or $`V_0/\eta _0`$ and $`T`$ is scaled with respect to $`V_0`$, where $`V_0`$ is half the potential barrier height (which is one). Period of the potential, $`L=2\pi `$.

## III Results and Discussions

### III.1 Diffusion Coefficient

Though the system response to the applied field is generally given by the stationary current density $`J=L<v>`$, but this directed motion (or average position of the particle) is accompanied by dispersion due to the inherent stochastic nature of the transport. It has been shown previously that in a homogeneous medium this dispersion or diffusion becomes very large (giant enhancement of DC) near the critical tilt. This enhancement in DC can be order of magnitude larger than the bare DC in the absence of potential. We make a systematic study of this phenomena in the presence of system inhomogeneities. Though our parameter space is large, we restrict to a narrow relevant domain where we observe effects which are surprising, and arise due to system inhomogeneities.

In fig. 1 we plot DC as function of external tilt $`F`$ for different values of temperature $`T`$ ($`\lambda =0.9\text{ and }\varphi =1.6\pi `$). It can be seen from the figure that DC exhibits a maxima as function of $`F`$. However, quantitative details depend sensitively on system parameters such as $`\varphi ,T\text{ and }\lambda `$. It can readily noticed from the curves A and B that DC has been enhanced by more than factor $`2`$ as compared with the homogeneous case. On lowering the temperature the relative enhancement of DC still increases. DC can also be suppressed by properly tuning $`\varphi `$. For $`\varphi =0.5\pi `$, DC near the critical field is suppressed by almost a factor $`2`$ as compared to the homogeneous case ( shown in the inset). Thus one can enhance or suppress DC near critical field by appropriately choosing the system parameters. It should be noted that in the case of enhancement, friction coefficient is lower on the immediate left side of the barrier and higher on the opposite side, where the relaxational motion takes place. When the frictional profile is opposite to the above case, suppression in DC occurs. When the phase difference between the two frictional profiles differ by $`\pi `$, then in one case enhancement and in the other suppression of DC can be observed as compared with the homogeneous case. Thus $`\varphi `$ controls the fluctuations of first passage times, hence DC and current.

Unlike the behavior in the homogeneous case ($`\lambda =0`$) where for small values of $`T`$ this peak value occur exactly at $`F=1`$ (critical field), here the peak position is very sensitive to $`\varphi `$ and can be shifted to either side of the bare critical field. Higher the temperature larger is the deviation from the critical threshold.

Fig. 2 shows this behavior, where we have plotted the force at which diffusion peaks ($`F_p`$) as function of $`\varphi `$ at $`\lambda =0.9\text{ for }T=0.1\text{ and }0.01`$. For $`T=0.1`$ the peak can occur at as low as $`F=0.8`$. The fact that $`F=0.8`$ is away from critical threshold, hence enhancement in DC in this regime has to be attributed to space dependent friction. This is a clear example of system inhomogeneity affecting the dynamical evolution of the particle nontrivially. This will be discussed at the end of this section to explain many of our observations. We would also like to emphasize that critical threshold is not altered at temperature $`T=0`$ for our present case as shown earlier Dan .

Since with increasing tilt the barrier to forward motion decreases (thereby reducing the effect of exponential suppression of DC in the presence of periodic potential), therefore it is natural to expect that $`D/D_0`$ will increase with increasing $`F`$ (for $`F<\text{ barrier height }`$). This is amply reflected in ref. Main , which corresponds to our $`\lambda =0`$ case. In the presence of space dependent friction ($`\lambda 0`$), $`D/D_0`$ can show a minimum (as shown in fig. 3) with increasing force ( for $`F<\text{ critical field }`$). This is surprising given the fact that current increases monotonically with increasing field (which we have checked separately) as expected. To observe this phenomena one has to properly choose the parameters. This is akin to the phenomena of coherent transport, wherein, increase in current is accompanied by decrease in DC. This we have discussed in detail in the later section as function of $`T`$, however, it is observed here as function of $`F`$ also.

Next we proceed to present the most interesting consequence of space dependent friction coefficient in our simple model. Unlike the expected result where the diffusion coefficient increases with temperature, here the diffusion can be suppressed by increasing the temperature. Fig. 4 highlights this. Here we have plotted $`D`$ as function of $`T`$. $`D`$ is scaled with $`V_0/\eta _0`$ as mentioned in the Sec. II. For the inset we have taken $`\varphi =1.36\pi \text{ and }1.44\pi \text{ at }F=0.98\text{ and }\lambda =0.9`$. For homogeneous case ($`\lambda =0`$), the minima is absent and $`D`$ increases monotonically with $`T`$. However, minima is clearly seen for the case $`\lambda 0`$. The observed suppression occurs when $`F`$ is close to the critical field. In can be clearly seen from the figure that as we go away from the critical field the minima in DC shifts to higher values of temperature and importantly it becomes shallower. Below certain value of $`F`$, minima and hence the suppression of DC disappears. The existence of the suppression in DC is attributed to the competition between two time scales. First, noise driven escape over potential barrier from the minima along the bias and the second time scale being the relaxation into the next potential well from the barrier top Fnl . It has been argued before that the second time scale is weakly dependent on noise strength and has a small variance as opposed to the first one. It is obvious that in transport processes when the second time scale dominates over the first it is expected to result in suppression of DC as function of noise strength (see for details ref. Fnl ). In the high temperature regime large thermal noise leads to large variance in the second time scale as the random motion of the particle both along and against the bias becomes equally important. Thus the DC increases as expected for higher temperatures, hence minima in DC. It is to be noted that for the case where $`\varphi =1.4\pi `$, the friction coefficient $`\eta (x)`$ is smaller near the barrier heights and moreover $`F`$ being close to critical threshold makes Arrhenius barrier crossing time scale (first time scale) smaller. As opposed to this, $`\eta (x)`$ is higher between the barrier height and the next potential minima along the bias thus slowing down the relaxation motion to the next minima. This naturally enhances the dominance of the relaxation time scale over the barrier crossing rate. This qualitatively explains our observed behaviour.

### III.2  Péclet number and coherent motion

By coherent motion we mean large particle current with minimal diffusion. This property can be quantified by the dimensionless Péclet number $`Pe`$ defined as Fnl

$$Pe=\frac{L\dot{x}}{D},$$
(9)
where $`L\dot{x}`$ is the average current density $`J`$. The expression for the current density is given in Eq. (7). We make use of expressions (5 and 7) to calculate $`Pe`$. The parameter values at which $`Pe`$ shows maxima correspond to the most coherent motion for that particular model. Higher the $`Pe`$ more coherent is the transport. It should be noted that the $`Pe`$ can show maxima though neither $`J`$ nor $`D`$ may show extrema. In fig. 5, we plot the current $`J`$, diffusion coefficient $`D`$ and the Péclet number $`Pe`$ (from top to bottom) for both space dependent (left column of fig. 5)and space independent friction (right column of fig. 5) cases. The average friction coefficient over a period in inhomogeneous case equals to the value taken for the homogeneous case. As pointed out earlier, though in the homogeneous case $`J\text{ and }D`$ is monotonic increasing function of $`T`$, the péclet number shows a maxima, the maximum value being very small compared to the space dependent friction case. For the chosen parameter values in the space dependent friction case, DC shows a minima with $`T`$ while the current increases monotonically. Moreover the magnitude of this current is larger than the corresponding current in the homogeneous friction case. This is indeed the most favorable condition of transport where increasing current is accompanied by decreasing DC. This is aptly reflected in the $`Pe`$, which shows enhancement (coherent motion or optimal transport) by an order of magnitude as compared with corresponding uniform friction case. The Péclet number is sensitively dependent on the phase factor $`\varphi `$ and it can also be suppressed (diffusion dominates the transport) which we have not reported here. Hence we can control the degree of coherent motion.

## IV  Conclusions

We have thus shown that both giant diffusion and coherent transport in a tilted periodic potential is sensitive to the frictional inhomogeneities of the medium. To analyze this problem we have taken a simple sinusoidal potential and periodic frictional profile with same periodicity but with a phase difference. Depending on the system parameters the value of DC near the critical threshold and Péclet number (indicating coherence in the transport) can be enhanced or decreased by an order of magnitude. Both these complimentary effects are important for practical applications Main ; Geier . The regime where we observe the optimal transport is accompanied by decrease of DC with temperature, the aspect which is absent in the corresponding homogeneous case. We have focussed on a restricted parameter regime to highlight the most interesting results arising due system inhomogeneities, in systems with simple potential. It is known that in the present model depending on system parameters current decreases with temperature, the effect akin to noise induced stability Nis . However, in this regime we have not observed any dramatic effect on DC as well as Péclet number. It is not clear whether noise induced stability (NIS) can enhance the coherence in the motion. The phenomenon of stochastic resonance (SR) in absence external ac field Sr is seen in this model. SR is characterized by the observation of peak in the particle mobility as function of system parameters such as $`T`$ and $`F`$ in certain parameter space. The analysis of Péclet number in this parameter space does not show any surprising features, so as to correlate with SR. This is due to the fact that SR occurs in the hight $`T`$ or high $`F`$ regime, where barriers to motion are absent. To observe the peak in the Péclet number the existence of barrier seems to be essential. To clarify the relation between SR, NIS and Péclet number one requires further detailed analysis.


format: text

text:
# Magnetic instability in a sheared azimuthal flow

## 1 Motivation

It is thought that turbulence in accretion discs arises from a magneto-rotational instability (MRI), where a magnetic field destabilises a rotating velocity field which decreases outwardly. This instability was discovered by Velikhov (1959) and Chandrasekhar (1961) when studying the motion of a fluid between two concentric cylinders (Taylor-Couette flow). It was only years later, when it was realized that Couette flow can be interpreted as model of Keplerian flow, that the implications for astrophysics were fully appreciated (Balbus & Hawley 1991). There have been many studies of the MRI recently, from numerical simulations of accretion discs (Brandenburg et al. 1995) to nonlinear calculations in spherical geometry (Drecker et al. 2000). In particular, in a recent paper Rüdiger & Zhang (2001) analysed the linear stability of hydromagnetic Couette flow and showed that, if the magnetic Prandtl number is less than unity, azimuthal Couette flow is more easily destabilised with a magnetic field than without. They found that the instability extends into the region of parameter space which is Rayleigh-stable without a magnetic field, which is important since Keplerian rotation is Rayleigh-stable. Rüdiger & Zhang were also able to study the MRI instability at magnetic Prandtl numbers Pm as small as $`0.001`$, towards the limit relevant to liquid sodium and gallium ($`\text{Pm}10^5`$) which are used in current MHD dynamo experiments (Tilgner 2000).

The aim of this paper to is to extend the investigation of Rüdiger & Zhang (2001). We explore the instability as a function of the speed of the outer cylinder in the Rayleigh-stable region, which is the parameter space of astrophysical interest (Rüdiger & Zhang considered only one nonzero ratio of outer to inner cylinder’s rotation). We also determine the effect of changing the magnetic Prandtl number, extending the range studied by Rüdiger & Zhang.

## 2 Model

We consider an incompressible fluid contained between two coaxial cylinders of inner radius $`R_1`$ and outer radius $`R_2`$ which rotate at prescribed angular velocities $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$. We use cylindrical coordinates $`(r,\varphi ,z)`$ and assume that a uniform magnetic field of strength $`B_0`$ is applied in the axial $`z`$ direction. At small angular velocities the flow is purely azimuthal (circular-Couette flow) and has magnitude $`V_0(r)=ar+b/r`$ which corresponds to the rotation law,

$$\mathrm{\Omega }(r)=\frac{V_0}{r}=a+\frac{b}{r^2},$$
(1)
The constants $`a`$ and $`b`$ are determined by the no-slip boundary conditions at the cylinder’s walls, so we have,

$$a=\mathrm{\Omega }_1\frac{(R_1^2/R_2^2\mathrm{\Omega }_2/\mathrm{\Omega }_1)}{(1R_1^2/R_2^2)},b=\mathrm{\Omega }_1R_1^2\frac{(1\mathrm{\Omega }_2/\mathrm{\Omega }_1)}{(1R_1^2/R_2^2)}.$$
(2)
At sufficiently high angular velocities Couette flow becomes unstable. The resulting velocity and magnetic fields $`V(r,\varphi ,z,t)`$ and $`B(r,\varphi ,z,t)`$ are determined by the MHD equations which we write in dimensionless form as,

$$\frac{V}{t}+(V\mathbf{})V=\mathbf{}p+^2V+\frac{Q}{\text{Pm}}(\mathbf{}B)B,$$
(3)
$$\frac{B}{t}=\frac{1}{\text{Pm}}^2B+\mathbf{}(VB),$$
(4)
$$\mathbf{}B=0,\mathbf{}V=0,$$
(5)
where $`p`$ is the pressure. In writing (3) and (4) we used $`\delta =R_2R_1`$ as unit of length, $`\delta ^2/\nu `$ as unit of time and $`B_0`$ as unit of magnetic field, where $`\nu `$ is the kinematic viscosity. The governing dimensionless parameters of the problem are the Reynolds numbers $`\text{Re}_1`$, $`\text{Re}_2`$, radius ratio $`\eta `$, and rotation ratio $`\mu `$,

$$\text{Re}_1=\frac{R_1\mathrm{\Omega }_1\delta }{\nu },\text{Re}_2=\frac{R_2\mathrm{\Omega }_2\delta }{\nu },\eta =\frac{R_1}{R_2},\mu =\frac{\mathrm{\Omega }_2}{\mathrm{\Omega }_1},$$
(6)
together with Chandrasekhar’s number $`Q`$ and the magnetic Prandtl number Pm,

$$Q=\frac{B_0^2\delta ^2\sigma }{\rho \nu },\text{Pm}=\frac{\nu }{\lambda },$$
(7)
where $`\rho `$ is the density, $`\sigma `$ is the electrical conductivity, $`\lambda =1/(\sigma \mu _0)`$ is the magnetic diffusivity and $`\mu _0`$ is the permeability.

We solve the MHD equations by direct time stepping from a small seeding initial condition and determine the stability of the Couette solution in different regions of parameter space. Our numerical method for 3D nonlinear flow is detailed in Willis & Barenghi (2002). For this linear study it suffices to say that the formulation is based on representing $`V`$ and $`B`$ with suitable potentials which are spectrally expanded over Fourier modes $`\mathrm{exp}[\mathrm{i}m\varphi +\mathrm{i}\alpha z]`$ in the azimuthal and axial directions and over Chebyshev polynomials in the radial direction. The axial wavelength of the disturbance is $`2\pi /\alpha `$. We assume no-slip boundary conditions for $`V`$, and electrically insulating boundaries conditions for $`B`$,

$$\frac{B_r}{z}=\frac{B_z}{_m}\frac{_m}{r},\frac{1}{r}\frac{B_z}{\varphi }=\frac{B_\varphi }{z},$$
(8)
where $`_m(r)`$ denotes the modified Bessel functions, $`I_m(\alpha r)`$ on $`R_1`$, and $`K_m(\alpha r)`$ on $`R_2`$. The time stepping is based on a combination of second order accurate Crank-Nicolson and Adams-Bashforth methods. The resulting computer code was tested against published results with and without magnetic field (Chandrasekhar 1961; Roberts 1964; Marcus 1984; Jones 1985; Barenghi 1991).

## 3 Results

It is well known that, in the absence of magnetic field, inviscid Couette flow is linearly stable provided that the celebrated Rayleigh criterion ($`\mu >\eta ^2`$) is satisfied.

Figure 1 shows the result of our calculations for the stability of dissipative Couette flow to axisymmetric perturbations, for radius ratio $`\eta =0.5`$ and magnetic Prandtl number $`\text{Pm}=1`$, as a function of the applied magnetic field. We plot the result in terms of $`\sqrt{Q}`$ rather than $`Q`$ in order to make direct comparison with the work of Rüdiger & Zhang (2001). The first curve refers to the case in which the outer cylinder is fixed, $`\mu =0`$, which is the most studied case in the fluid dynamics literature. It is apparent that the presence of a magnetic field makes the flow more unstable. The critical Reynolds number, which is $`\text{Re}_{1c}=68.2`$ for $`Q=0`$, decreases with increasing $`Q`$ and has a minimum at $`Q=39`$. The most unstable mode is $`m=0`$ over the range for $`Q`$ in Fig. 1. From here we consider axisymmetric disturbances only. The critical axial wavenumber $`\alpha `$ decreases significantly from $`3.1`$ to $`1.7`$ over the range, and varies like $`1/\sqrt{Q}`$ thereafter. This stiffening eventually restabilises the flow and $`\text{Re}_{1c}`$ increases like $`\sqrt{Q}`$ for strong fields.

The initial destabilisation is consistent with the finding of Rüdiger & Zhang (2001); the small difference between their $`\text{Re}_{1c}`$ and ours is certainly due to the different boundary conditions for $`B`$ (they assumed pseudo-vacuum conditions and we assume insulating conditions). The second curve of Fig. 1 refers to the case $`\mu =\eta ^2`$ (the Rayleigh line), which separates stable and unstable regions in the absence of a magnetic field. The curve well illustrates the striking destabilising effect of the magnetic field.

Figure 2 shows stability boundaries in the $`\text{Re}_1`$ vs $`\text{Re}_2`$ plane for different values of the imposed magnetic field. The axial wavenumber $`\alpha `$ does not vary a great deal along the boundaries, its dependence being principally determined by the strength of the field $`Q`$.

It is apparent that even a small value of $`Q`$ is enough to make the boundary cross the Rayleigh line (the upper dotted curve in the figure). The destabilising effect of the magnetic field is so large that the stability boundary drops toward the region of solid body rotation ($`\mathrm{\Omega }_1=\mathrm{\Omega }_2`$ or $`\mu =1`$), which is the lower dotted line. However, if the applied field is strong enough the flow can be restabilised, in accordance with the results of Fig. 1.

The destabilisation becomes more important the larger the magnetic Prandtl number Pm. Figure 3 shows results for a moderate value of applied magnetic field, $`Q=10`$, at increasing values of Pm. It is apparent that the stability boundary drops much below the Rayleigh line $`\mu =\eta ^2`$ and, for large enough Pm, becomes asymptotic to the solid body rotation line $`\mu =1`$.

Finally, in Fig. 4 we find that on the Rayleigh line the stability follows a power law,

$$\text{Re}_{1c}\text{Pm}^\beta ,\beta =0.500\pm 0.002,$$
(9)
over a surprising range of Pm. The error in $`\beta `$ is based on a fit through all points calculated. For $`\eta =0.5`$ the minimum $`\text{Re}_{1c}`$ occurs at $`Q=86`$ and $`\alpha =2.1`$. To verify this result we use the WKB analysis of Ji et al. (2001). The local dispersion relation they find for the Taylor-Couette flow is identical to the one derived for accretion discs in the incompressible limit (Sano & Miyama 1999). From the dispersion relation the necessary and sufficient condition for stability to axisymmetric perturbations is,

$$(\text{Pm}+S^2)^2(1+ϵ^2)+2\zeta \text{Rm}^22(2\zeta )\text{Rm}^2S^20,$$
(10)
where $`S`$ is the ratio of diffusive to Alfvén timescales, Rm is the magnetic Reynolds number, and $`\zeta =1/(r\mathrm{\Omega })_r(r^2\mathrm{\Omega })`$. On the Rayleigh line, from (1) and (2), $`\mathrm{\Omega }(r)`$ takes the simple form, $`\mathrm{\Omega }=b/r^2`$, and so $`\zeta =0`$ at all $`r`$. We said earlier that in Fig. 2 the wavenumber changes little along the boundary, as its dependence is principally on the strength of the imposed field. The wavenumber parameter $`ϵ`$ of Ji et al. (2001) then depends only on $`Q`$, so $`ϵ=ϵ(Q)`$. Ignoring factors that depend only on $`ϵ`$ in the substitutions,

$$S^2=\text{Pm}Q,\text{Rm}=\text{Pm}\text{Re},$$
(11)
(10) leads to the condition for stability,

$$\text{Re}_1\frac{A(Q)}{\sqrt{\text{Pm}}},$$
(12)
which confirms our result (9). This relation (12) should hold for all Pm as we have not been required to take any limits. The function $`A(Q)`$ has a minimum for some $`Q`$ and $`\alpha (Q)`$. For $`\eta =0.5`$ we found the minimum occurs at $`Q=86`$, $`\alpha =2.1`$.

Whilst the linear result (12) should hold for very small Pm, from the work of Wendt (1933) and Taylor (1936) a turbulent (nonlinear) instability is observed at large Reynolds numbers ($`\text{Re}_12.5\times 10^5`$) with this radius ratio. Extrapolating down one magnitude to $`\text{Pm}=10^5`$ for laboratory fluids gives $`\text{Re}_{1c}2\times 10^5`$, a similar value, for the magneto-rotational (linear) instability.

A radial truncation of 12 Chebyshev modes was found sufficient for all calculations. In particular, for the calculation at $`\text{Pm}=10^4`$ in Fig. 4, convergence was tested by increasing the truncation and also decreasing the timestep. The fractional numerical error in $`\text{Re}_{1c}`$ is estimated at approximately $`10^6`$.

## 4 Discussion

Our calculations show that many rotation laws of the form $`\mathrm{\Omega }(r)=a+b/r^2`$ which are hydrodynamically stable (that is to say, they satisfy the Rayleigh criterion) become linearly unstable when a magnetic field is applied. Our results confirm the finding of Rüdiger & Zhang (2001) and extend them in the Rayleigh stable region.

We have determined the instability at magnetic Prandtl numbers Pm one order of magnitude smaller than Rüdiger & Zhang’s, towards the small magnetic Prandtl number limit, which is relevant to possible MHD dynamo experiments with liquid sodium and gallium. Although the power law $`\text{Re}_{1c}\text{Pm}^{0.5}`$ that we find on the Rayleigh line ($`\mu =\eta ^2`$) is slightly different from theirs ($`\text{Re}_{1c}\text{Pm}^{0.65}`$ on $`\mu =\frac{1}{3}`$), it confirms their conjecture that the nonlinear instability found by Richards & Zahn (1999) and the MRI are likely to occur at Reynolds numbers of the same order of magnitude.

We also find that the flow becomes particularly unstable if the magnetic Prandtl number is greater than unity. The instability boundary in the $`\text{Re}_1`$ vs $`\text{Re}_2`$ plane rapidly tends towards the solid body rotation line. This enhanced instability for large Pm is consistent with earlier results of Kurzweg (1963). His boundary conditions were selected such as to avoid mathematical difficulties but for small Pm agreed well with the results of Chandrasekhar (1961). The significance of the instability in this case is linked to the possibility (Kulsrud & Anderson 1992; Brandenburg 2001) that large values of Pm exist in central regions of galaxies.


text:
# Stellar Parameters of Wolf-Rayet Stars from Far-UV to Mid-IR Observations

## 1. Recent Observational and Theoretical Progress

This article will focus on recent determinations of physical parameters for Galactic and Magellanic Cloud WR stars from UV to IR diagnostics. Observationally, the Far-Ultraviolet Spectroscopic Explorer (FUSE) has provided an impressive database of $`\lambda `$912–1187Å spectroscopy for Galactic and Magellanic Cloud WR stars (see Willis et al. these proc.) to supplement previous UV datasets obtained with the International Ultraviolet Explorer (IUE) and Hubble Space Telescope (HST). At longer wavelengths, Infrared Space Observatory (ISO) observations of WR stars have now been analysed. Much of the recent observational progress with Wolf-Rayet stars has involved the acquisition of high quality optical spectroscopy for individual stars beyond the Magellanic Clouds (Drissen, these proc.), plus X-ray spectroscopy of single and binary WR stars (e.g. Skinner et al. 2002) neither of which topics will be discussed here.

Theoretical developments in the last few years have been more steady, with the (laborious) implementation of line blanketing into codes by elements other than CNO and Fe, which had already been discussed at the last hot star beach symposium, IAU Symp. 193. The major change has been the widespread use of such codes to analyse individual stars within a range of galaxies. At present, there are a variety of model atmosphere codes which consider sphericity and line blanketing and fall into two main types, outlined below.

CMFGEN (Hillier & Miller 1998) and the Gräfener et al. (2002) code make use of variants of the super-level approach to incorporate the effect of tens of thousands metal lines on the atmospheric structure within the radiative transfer code. CMFGEN can now simultaneously consider blanketing by individual ions of up to 30 elements, including C, N, O, Ne, Si, S, Ar, Ca, Fe and Ni (Hillier, these proc.), whilst Gräfener et al. consider CNO, Si plus Fe-group elements (Sc to Ni) grouped together in a single generic atom. This approach suffers the least number of approximations, but remains computationally demanding. Recent test calculations for early-type WC stars show (perhaps surprisingly!) good consistency between these two codes, including ionizing fluxes.

Alternatively, Schmutz (1997) and ISA-wind (de Koter et al. 1993, 1997) use separate codes to solve the radiative transfer problem and line blanketing, the latter making use of Monte Carlo techniques. The method had the great computational advantage that complete intensity-weighed effective opacity factors can be calculated separately from the transfer problem. On the negative side, the ionization and excitation equilibrium of metal species is approximate, dictating which lines are efficient at capturing photons for each point in the atmosphere. Test calculations for a late-type WN star between CMFGEN and ISA-wind show excellent agreement in derived stellar parameters, but rather poorer agreement for ionizing fluxes (Crowther et al. 1999).

The main effect of blanketing is to re-distribute extreme UV flux to longer wavelengths, reducing the ionization balance in the atmosphere, such that higher stellar temperatures (and luminosities) are required to match observed line profile diagnostics relative to unblanketed studies. This is illustrated in Table 1 for the prototypical Galactic early-type WC star HD 165763 (WR111) whose derived stellar luminosity has increased by a factor of 5–7 over the past decade. Differences in luminosities for HD 165763 between the recent studies of Hillier & Miller (1999) and Gräfener et al. (2002) most likely result from the inclusion of additional blanketing elements, which CMFGEN now routinely handles. Recent revisions to temperatures and luminosities of O supergiants have acted in the opposite sense, relative to previous plane-parallel unblanketed model analyses, such that common techniques are now employed throughout for O and WR stars (e.g. Crowther et al. 2002ab).

Clumping is now routinely, albeit approximately, handled in WR model atmospheric studies via radial dependent volume filling factors, $`f`$. Constraints on $`f`$ can be obtained from comparisons between red electron scattering wings and observations (e.g. Hillier 1991), although exact determinations generally prove elusive due to line blending, particularly in WC stars. Generally, $`f`$0.05–0.25 provide reasonable matches to observed line profiles (e.g. Hamann & Koesterke 1998), such that global mass-loss rates are reduced by a factor of $`1/\sqrt{f}24`$ relative to smooth models. The majority of line profiles behave rather insensitively provided $`\dot{M}/\sqrt{f}`$ remains constant with some exceptions (e.g. Herald et al. 2001). In WC stars, the line centre of UV resonance lines of C iii-iv reacts to changes in $`f`$, providing additional constraints on the filling factor (Crowther et al. 2002a).

## 2. WN properties

Although large samples of WN stars have not yet been thoroughly analysed using recent line blanketed codes, stellar temperatures of WN stars range from 30kK (at WN10), to $``$40kK (at WN8) and approaching 100kK for early-type WN (WNE) stars. Increased stellar temperatures, particularly for WNE stars, with correspondingly smaller radii, brings atmospheric models into much closer agreement with direct determinations from short period WN+O binaries, such as V444 Cyg (e.g. Moffat & Marchenko 1996). Clumped mass-loss rates from optical diagnostics generally lie in the range $`\dot{M}=10^{5.5\mathrm{}4.5}M_{}`$ yr-1. Mass-loss rates can also be obtained from radio determinations, although these too are subject to uncertainties in volume filling factors, and are limited to stars within a few kpc. Typical wind velocities span a wide range, well correlated with spectral type: 300 km s-1 at WN10 to $``$2000 km s-1 at WN3-4.

Relative to recent unblanketed results, Herald et al. (2001) obtained a smaller core radius and slightly enhanced luminosity for HD 96548 (WR40, WN8) allowing for blanketing including Ne, Ar, Ca and Fe, such that the heavily blended Fe iv-v forest in the UV is very well reproduced. This naturally leads to the potential of late-type WN (WNL) stars as providing diagnostics of heavy elemental abundances from UV and far-UV (FUSE) spectroscopy. Morris et al. (2000) carried out a combined optical–infrared analysis of another WN8 star AS 431 (WR147) revealing overall excellent consistency between visual, near-IR and mid-IR ISO line diagnostics, such that reliable stellar parameters may be obtained from solely IR observations, providing all necessary diagnostics are covered. Meanwhile, Crowther et al. (1999) used the H ii ejecta nebula M1–67 as a sensitive diagnostic of the ionizing flux distributions predicted by analyses of the central star 209 BAC (WR124, WN8) using CMFGEN and ISA-wind, favouring the former in this case.

In general, WNL stars are relatively H-rich (H/He$``$0.5-2 by number), whilst WNE stars do not contain hydrogen, although exceptions do occur. The question of whether weak-lined WN5–6 stars in R136 are moderately enriched in helium (H/He$``$3, Crowther & Dessart 1998) or unprocessed (H/He$``$10, de Koter et al. 1997) remains open. Where studies of WN stars have been carried out, N is extremely enriched and C depleted, consistent with CN-cycle processed material, with O generally more difficult to constrain (e.g. Herald et al. 2001). In principal, other elemental abundances may be obtained from UV synthesis, but the most direct method is via analysis of mid-IR fine structure lines. Morris et al. (2000) obtained lower limits of $``$0.5 solar for the abundance of Ne, S and Ca in AS 431 from ISO observations, with individual determinations sensitive to the adoption of smooth or clumped winds (see e.g. J.D. Smith, these proc.).

WN spectral types mimic underlying stellar temperatures rather well, at least within a particular (metallicity) environment. However, since WN subtypes depend on the relative strength of nitrogen emission lines, samples within different galaxies do not necessarily behave in the same manner. This is because N iii-v classification lines have differing sensitivities to abundance which will be a factor of five times lower for a SMC WN star than its counterpart in the Solar neighbourhood. Crowther (2000) demonstrated that a WN5 star in the SMC would have a later spectral type in the Milky Way if its stellar parameters were kept constant, and an earlier spectral type at lower metallicity. This tendency broadly explains the observed trend earlier WN spectral types when one compares the statistics of WN stars in the Galaxy and Magellanic Clouds.

## 3. WC properties

Due to the overwhelming effect of metal line blanketing on the atmospheric structure of WC stars, considerable effort has recently gone into their analysis. Stellar temperatures range from $``$45kK (WC9) to $``$60kK (WC8) to $``$100kK (WCE) and likely well in excess of this value for some WO stars. Mid-IR ISO observations (e.g. De Marco et al. 2000) and far-UV FUSE spectroscopy (e.g. Crowther et al. 2000) have been combined with the usual UV and optical diagnostics to derive stellar parameters. The use of similar techniques for Galactic (e.g. Dessart et al. 2000) and Magellanic Cloud (Crowther et al. 2002a) WC stars permits their relative properties to be investigated in detail. Clumped WC mass-loss rates lie in the range $`\dot{M}=10^{5.0\mathrm{}4.4}M_{}`$ yr-1. As for WN stars, wind velocities are well correlated with spectral type: 1100 km s-1 at WC9, to 3000 km s-1 at WC4.

Carbon abundance determinations of WC stars have long used the He ii $`\lambda `$5411 and C iv $`\lambda `$5471 recombination lines. However, conflicting results may result if limited model atoms are used, or if trace elements are neglected, such that Gräfener et al. (1998) obtained C/He=0.32 by number (40% by mass) for HD 32125 (WC4), in contrast with the recent determination of C/He=0.13 by number (25% by mass) from Crowther et al. (2002a). Carbon abundances in Galactic and LMC WC stars span a similar range, 0.1 $``$ C/He $``$ 0.4, such that spectral types are not determined by carbon abundances, in conflict with predictions by Smith & Maeder (1991). Hamann et al. (these proc.) suggest a relatively uniform carbon abundance of C/He$``$0.25 for most Galactic WCE stars. Sand 2 (LMC, WO) does appear to have a higher C/He ratio than WC stars, but determinations for such stars are hindered by severe blending caused by it’s broad lines.

Oxygen abundances are rather more difficult to tightly constrain, since most diagnostics lie around $`\lambda `$3000Å, and span a much wider range of ions – O iii-vi in most WC stars – than carbon. Indeed, the classification line at $`\lambda `$5590 is formed from a blend of O iii and O v lines. Consequently, O/He determinations can generally not be obtained better than a factor of $``$2. Oxygen in Galactic and LMC WC stars ranges from 0.02 $``$ O/He $``$0.1, again with the probable exception of Sand 2 (Dessart et al. 2000; Crowther et al. 2000, 2002a). As for WN stars, other elemental abundances are most readily determined from mid-IR fine structure lines, provided clumping factors can be determined. For $`\gamma `$ Vel (WC8+O, WR11), Dessart et al. (2000) find that S is within 20% of the cosmic value, whilst Ne is enhanced by a factor of $``$8, in reasonable agreement with recent evolutionary predictions.

If WC spectral types do not result from chemical abundance variations, why are LMC WC stars systematically of earlier spectral type than Galactic WC stars? Temperature certainly plays a role, but for an answer we must first consider the WC classification lines – C iii $`\lambda `$5696 and C iv $`\lambda \lambda `$5801-12 – in greater detail. The upper level of $`\lambda `$5696 has an alternative decay via $`\lambda `$574 with a branching ratio of 147:1. Consequently, $`\lambda `$5696 only becomes strong when $`\lambda `$574 is optically thick, i.e. if the stellar temperature is low or the wind density is high. A comparison between Galactic WC5–7 stars and LMC WC4 stars reveals that their temperatures are comparable, such that the wind densities of the Galactic sample must be higher than the LMC stars. Crowther et al. (2002a) demonstrate that this is indeed the case, such that the principal difference amongst WCE stars is differing wind densities – high for WC6-7 stars and low for WC4 stars, with

$$\mathrm{log}(\dot{M})=1.38\mathrm{log}(L/L_{})12.35$$
(1)
for the LMC stars, and Galactic WC stars $``$0.25 dex higher, in good agreement with the generic WR mass-loss luminosity calibration of Nugis & Lamers (2000).

## 4. A Metallicity dependence for Wolf-Rayet stars?

A wind density origin for the subtype distribution amongst WC stars in Local Group galaxies can be understood if the mass-loss rates of Wolf-Rayet stars are metallicity dependent. The difference between Galactic and LMC WCE stars amounts to only a factor of two, which is consistent with a metallicity dependence of $`Z^{0.50.7}`$ as predicted by radiation driven winds for O stars (e.g. Vink et al. 2001). Reduced wind densities for WC stars at lower metallicity, preferentially affecting C iii $`\lambda `$5696, naturally explains the trend towards WC4 (and WO) stars, as is observed in the Magellanic Clouds and IC10 (Crowther et al. these proc). Although the atmospheres of WC stars are composed of mostly He, C and O, the heavier elements (Ne, Ar, Fe, Ni…) initiate their winds, as in O stars.

A WR metallicity dependence is not currently adopted in evolutionary calculations, which generally typically adopt $`\dot{M}=\beta \times 10^7(M/M_{})^\alpha `$ for all H-free WR stars (Langer 1989). For WC stars Langer adopted $`\alpha `$=2.5 and $`\beta `$=1.0. This relationship is compared to recent mass-loss rates and masses for Galactic and LMC WC stars in Fig. 1, revealing very poor agreement. LMC WC stars are well matched using $`\alpha `$=2.0 and $`\beta `$=0.9, similar to that discussed by Cherepashchuk (2001) from independent methods.

If WC stars are metallicity dependent, one would expect the same of WN stars. However, results to date for WN stars has been less clear (Crowther & Smith 1997; Hamann & Koesterke 2000), most likely due to greater scatter in wind properties at a given spectral type. Nevertheless, the trend is also towards weaker winds, as is seen most readily via the single WN stars in the SMC. The question of whether WR mass-loss rates are metallicity dependent has a major impact on the hardness of their ionizing flux distributions, as discussed by L.J. Smith et al. (these proc.). Denser winds cause softer Lyman continuum fluxes, which would imply that few WR stars produce nebular He ii ionizing photons at high metallicity, with the opposite true at low metallicities.

## 5. Pre-Supernova WR Masses

Recent theoretical progress in our analysis of Wolf-Rayet stars has led us to the claim that derived stellar parameters are ‘robust’. Since there is a well established theoretical mass-luminosity relation for hydrogen-free WR stars, those which are members of close period binaries offer us the possibility to check derived spectroscopic luminosities. De Marco et al. (2000) obtained a stellar luminosity of the WC8 component of $`\gamma `$ Vel which implies a mass of 9$`M_{}`$ (via the Heger & Langer 1996 relationship), in remarkable agreement with 9.5$`M_{}`$ derived from the binary orbit.

We can employ evolutionary models to obtain remaining lifetimes, which together with an appropriate mass-loss luminosity relation (e.g. Eqn. 1), permits an estimate of the final, pre-SN mass. This approach can be followed for all WR stars. However, uncertainties in evolutionary models and lack of recent results for WN stars suggest that WC stars represent our best candidates for the determination of pre-SN masses. Taking the LMC WC4 star HD 32402 as an example, Crowther et al. (2002a) derive a current mass of 21$`M_{}`$, a remaining lifetime in the range 1–4$`\times 10^5`$yr, and so a pre-SN mass of 14$`M_{}`$ (non-rotating model) or 19$`M_{}`$ (rotating model). From a sample of Galactic WC stars at known distance, final pre-SN masses lie in the range 7–14$`M_{}`$, versus 11–19$`M_{}`$ using LMC stars.

A number of Type-Ic supernovae that have been studied in the past few years have probable WR precursors, including SN 1998bw, for which a CO (ejected) core mass of 13.8$`M_{}`$ was derived by Iwamoto et al. (1998). This agrees rather well with the expected pre-SN masses of LMC WC stars. There has been considerable interest in SN 1998bw since it was also GRB 980425. SN 2002ap is another nearby Type-Ic for which a WR star is one of the few likely progenitors (Smartt et al. these proc.). Unfortunately, most hypernovae or ‘collapsar’ models require rapid rotation, whilst rotation rates are lacking for almost all WR stars. The general consensus is that most WR stars have spun-down – indeed practically all WC stars are thought to be spherical, as inferred from spectropolarimetry of their winds (Harries et al. 1998).

### Acknowledgments.

Many thanks to the Royal Society for continued financial support, and to my collaborators without whom the work presented here would not have been possible, John Hillier and Luc Dessart in particular.

## References

Cherepashchuk A.M., 2001, Astron. Rep. 45, 120

Crowther P.A., 2000, A&A 356, 191

Crowther P.A., Dessart L. 1998, MNRAS 296, 622

Crowther P.A., Smith L.J., 1997, A&A 320, 500

Crowther P.A., Pasquali A., De Marco O., et al. 1999, A&A 350, 1007

Crowther P.A., Fullerton A.W., Hillier D.J. et al. 2000, ApJ 538, L51

Crowther P.A,, Dessart, L., Hillier D.J. et al., 2002a, A&A in press (astro-ph/0206233)

Crowther P.A., Hillier D.J., Evans C.J. et al., 2002b, ApJ in press (astro-ph/0206257)

de Koter, A., Schmutz W., Lamers H.J.G.L.M., 1993, A&A 277, 561

de Koter, A., Heap, S.R., Hubeny, I. 1997, ApJ 477, 792

De Marco, O., Crowther P.A., Schmutz W., et al., 2000, A&A 358, 187

Dessart L., Crowther P.A., Hillier, D.J. et al., 2000, MNRAS 315, 407

Gräfener, G., Hamann, W-.R., Hillier, D.J., Koesterke, L., 1998, A&A, 329, 190

Gräfener, G., Koesterke, L., Hamann, W-R., 2002, A&A 387, 244

Hamann W.-R., Koesterke L., 1998, A&A 335, 1003

Hamann W.-R., Koesterke L., 2000, A&A 360, 647

Harries T.J., Hillier D.J., Howarth I.D., 1998 MNRAS 296, 1072

Heger A., Langer N., 1996 A&A 318, 421

Herald J.E., Hillier D.J., Schulte-Ladbeck, R.E., 2001, ApJ 548, 932

Hillier D.J., 1989, ApJ 347, 392

Hillier D.J., 1991, A&A 247, 455

Hillier D.J., Miller D.L., 1998, ApJ, 496, 407

Hillier D.J., Miller D.L., 1999, ApJ, 519, 354

Iwamoto K., Mazzali, P.A., Nomoto K., et al. 1998, Nat 395, 672

Langer N., 1989, A&A 220, 135

Moffat A.F.J., Marchenko S.V., 1996, A&A 305, L29

Morris P.W., van der Hucht K.A., Crowther P.A. et al. 2000, A&A 353, 624

Nugis T., Lamers H.J.G.L.M. 2000, A&A 360, 227

Schmutz W., 1997, A&A 321, 268

Schmutz W., Hamann, W.-R., Wessolowski, U. 1989, A&A 210, 236

Skinner S.L., Zhekov, S.A., Güdel, M., Schmutz W., 2002, ApJ 572, 477

Smith L.F., Maeder A. 1991, A&A 241, 77

Vink J.S., de Koter A., Lamers H.J.G.L.M, 2001, A&A 369, 574


text:
# Local Permutations of Products of Bell States and Entanglement Distillation

## I Introduction

We study mixed-state multi-copy entanglement distillation protocols for pairs of qubits. We start from $`n`$ identical copies of a Bell diagonal state of $`2`$ qubits and end up, after local operations and classical communication, with $`m<n`$ Bell diagonals (possibly statistically dependent) with higher joint fidelity than $`m`$ copies of the original Bell diagonal state. For non-Bell-diagonal initial states one can first perform $`n`$ separate optimal single-copy distillation protocols to make them Bell-diagonal Verstraete et al. (2001). Our protocol can be used in a recurrence scheme followed by the hashing protocol as in Bennett et al. (1996a, b). We propose a protocol with $`n=4`$ and $`m=1`$ that does significantly better than existing protocols. Our results can be used to find even better protocols for other values of $`n`$ and $`m`$ that are tuned for specific initial states.

We see three main reasons for studying entanglement distillation protocols. The first and most obvious reason is that entanglement distillation protocols are a means of obtaining states that are closer to maximally entangled pure states, as needed in typical applications like teleportation, from mixed states that can be reached by sending one qubit of an entangled pair through a realistic channel. A second reason to study distillation protocols is that asymptotic protocols yield a lower bound for “entanglement of distillation”, an important measure of entanglement, that is in itself a lower bound for any sensible measure of entanglement Horodecki (2001). In this context we also mention the upper-bounds on entanglement of distillation obtained in Rains (2001). A third reason is that multi-copy entanglement protocols can be considered as applications of entanglement, where more can be done in the presence of entangled pairs than without. We hope that studying these mechanisms will reveal some information on the important problem of how exactly the presence of entanglement enables one to do things that are impossible without.

Multi-copy mixed state entanglement distillation for qubit pairs was first studied in Bennett et al. (1996a, b). An improved variant of the two-copy protocol in that paper was described in Deutsch et al. (1996) under the title of quantum privacy amplification. These protocols as well as ours start from $`n`$ identical qubit pairs in a Bell-diagonal state, shared by Alice and Bob. A crucial ingredient of these protocols is a local unitary operation, performed by Bob and Alice on their $`n`$ qubits, that results globally in a permutation of the $`4^n`$ possible tensor products of Bell states. The key ingredient of this paper is a characterization, by means of a binary matrix group, of all possible local permutations of the products of Bell states. This enables a search for the best protocol within this setting. In section II we study local permutations of products of Bell states. In section III we discuss the new protocols. In section IV we discuss the combination of our protocols with a recurrence scheme and the hashing protocol and show the strength of our protocols by computer simulations.

## II Local Permutations of products of Bell states

In this section we study the class of local unitary operations that can be performed by Alice and Bob locally and result in a permutation of the $`4^n`$ (tensor) products of Bell states, where $`n`$ is the number of qubit pairs. These local permutations are the key ingredient of the new distillation protocols described in the next section.

We will code products of Bell states by binary vectors by assigning two-bit vectors to the Bell states as follows

$$\begin{array}{ccccc}|\mathrm{\Phi }^+\hfill & =& \frac{1}{\sqrt{2}}(|00+|11)\hfill & =& |B_{00}\hfill \\ |\mathrm{\Psi }^+\hfill & =& \frac{1}{\sqrt{2}}(|01+|10)\hfill & =& |B_{01}\hfill \\ |\mathrm{\Phi }^{}\hfill & =& \frac{1}{\sqrt{2}}(|00|11)\hfill & =& |B_{10}\hfill \\ |\mathrm{\Psi }^{}\hfill & =& \frac{1}{\sqrt{2}}(|01|10)\hfill & =& |B_{11}.\hfill \end{array}$$
(1)
A product of $`n`$ Bell states is then described by a $`2n`$-bit vector, e.g. $`|B_{001101}=|B_{00}|B_{11}|B_{01}=|\mathrm{\Phi }^+|\mathrm{\Psi }^{}|\mathrm{\Psi }^+`$.

We will also exploit a correspondence between Bell states and Pauli matrices

$$\begin{array}{ccc}|\mathrm{\Phi }^+\hfill & & \frac{1}{\sqrt{2}}\sigma _{00}=\frac{1}{\sqrt{2}}\sigma _0=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 1\end{array}\right]\hfill \\ |\mathrm{\Psi }^+\hfill & & \frac{1}{\sqrt{2}}\sigma _{01}=\frac{1}{\sqrt{2}}\sigma _x=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}\hfill 0& \hfill 1\\ \hfill 1& \hfill 0\end{array}\right]\hfill \\ |\mathrm{\Phi }^{}\hfill & & \frac{1}{\sqrt{2}}\sigma _{10}=\frac{1}{\sqrt{2}}\sigma _z=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 1\end{array}\right]\hfill \\ |\mathrm{\Psi }^{}\hfill & & \frac{1}{\sqrt{2}}\sigma _{11}=\frac{1}{\sqrt{2}}\sigma _y=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}\hfill 0& \hfill i\\ \hfill i& \hfill 0\end{array}\right]\hfill \end{array}$$

A tensor product of $`n`$ Bell states is then described by a tensor product of Pauli matrices, e.g. $`|\mathrm{\Phi }^+|\mathrm{\Psi }^{}|\mathrm{\Psi }^+\frac{1}{\sqrt{8}}\sigma _0\sigma _y\sigma _x=\frac{1}{\sqrt{8}}\sigma _{00}\sigma _{11}\sigma _{01}`$. We will also use longer vector subscripts to denote such tensor products, e.g.

$$\sigma _{001101}=\sigma _{00}\sigma _{11}\sigma _{01}.$$
(2)
In this representation of pure states of $`2n`$ qubits as $`2^n\times 2^n`$-matrices $`\stackrel{~}{\mathrm{\Psi }}`$, local unitary operations $`|\psi (U_AU_B)|\psi `$, in which Alice acts on her $`n`$ qubits (jointly) with an operation $`U_A`$ and Bob on his with an operation $`U_B`$, are represented by

$$\stackrel{~}{\mathrm{\Psi }}U_A\stackrel{~}{\mathrm{\Psi }}U_B^T.$$
(3)
We are now in a position to state the main result of this section:

###### Theorem 1

(i) If a local unitary operation (3) results in a permutation of the $`4^n`$ tensor products of $`n`$ Bell states, this permutation can be represented on the binary vector representations (1) as an affine operation

$$\begin{array}{cc}& \varphi :_2^{2n}_2^{2n}:xAx+b\hfill \\ \text{with }\hfill & A_2^{2n\times 2n},b_2^{2n}\hfill \\ \text{and }\hfill & A^TPA=P,\hfill \\ \text{where }\hfill & P=\text{diag}[\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],\mathrm{},\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]]\hfill \end{array}$$
(4)
(ii) Conversely, any such permutation can be realized by local unitary operations.

Note that all multiplication and addition should be done modulo 2. We call a matrix $`A`$ satisfying $`A^TPA=P`$ $`P`$-orthogonal. The affine and linear transformations considered are invertible and therefore amount to a permutation of $`_2^{2n}`$. In the sequel we sometimes directly refer to the linear transformations as permutations.

Proof: We first prove part (i). One can easily check that $`\sigma _v\sigma _w\sigma _{v+w}`$, where $`v`$ and $`w`$ are binary vector indices as in (2), and the $``$-sign means equal up to a complex phase (in this case $`1`$,$`i`$,$`1`$ or $`i`$). Such a phase is irrelevant as the Pauli matrices here are matrix representations of pure state vectors.

Assume now that $`U_A`$ and $`U_B`$ indeed result in a permutation $`\pi :_2^{2n}_2^{2n}`$, then the null vector is mapped to some vector $`v=\pi (0)`$. Accordingly $`U_A\sigma _0U_B^T\sigma _v`$. Since $`\sigma _0`$ is the identity matrix, we have $`U_B\sigma _vU_A^{}`$ where  denotes complex conjugation. If we want to represent $`\pi `$ by $`xAx+b`$ we clearly have to choose $`b=v`$. Note that (3) now reads $`\stackrel{~}{\mathrm{\Psi }}U_A\stackrel{~}{\mathrm{\Psi }}U_A^{}\sigma _b.`$

We now have to show that the permutation $`\pi ^{}:x\pi (x)+b`$, which maps $`\stackrel{~}{\mathrm{\Psi }}`$ to $`U_A\stackrel{~}{\mathrm{\Psi }}U_A^{}`$, is a linear $`P`$-orthogonal map $`\pi ^{}:xAx`$. Linearity of binary maps means that sums are mapped to sums $`\pi ^{}(v+w)=\pi ^{}(v)+\pi ^{}(w)`$. This is clearly true since $`U_A\sigma _{v+w}U_A^{}U_A\sigma _v\sigma _wU_A^{}=U_A\sigma _vU_A^{}U_A\sigma _wU_A^{}`$. Furthermore, it can be verified using the commutation and anticommutation laws for Pauli matrices, that $`\sigma _v`$ and $`\sigma _w`$ are commutable operators if and only if $`v^TPw=0`$. Since $`\sigma _v`$ and $`\sigma _w`$ are commutable if and only if $`U_A\sigma _vU_A^{}`$ and $`U_A\sigma _wU_A^{}`$ are commutable, it must hold that $`v^TA^TPAw=v^TPw`$ for all $`v`$ and $`w`$, which proves $`A^TPA=P`$.

To prove part (ii), we will first consider $`n=2`$ and show that all permutations $`\pi :xAx`$ with $`A^TPA=P`$ can be generated with the operations $`\varphi _u:\stackrel{~}{\mathrm{\Psi }}U_u\stackrel{~}{\mathrm{\Psi }}U_u^{}`$ with $`U_u=e^{i\frac{\pi }{4}\sigma _u}=\frac{1}{\sqrt{2}}(I+i\sigma _u)`$ (with $`u_2^4`$). (This is also true for $`n>2`$, but generators affecting more than $`2`$ qubits at a time will not be needed). Using $`\sigma _v\sigma _w=(1)^{v^TPw}\sigma _w\sigma _v`$, it can be shown that $`\varphi _u`$ translated into the binary language results in a permutation $`\pi _u:xx+u(u^TPx)=(I+uu^TP)x`$.

We will now show that the group of permutations generated by the permutations $`\pi _u`$ is isomorphic to $`S_6`$, the group of all permutations of $`6`$ elements. Next we will show that the group of $`P`$-orthogonal $`4\times 4`$ matrices contains $`6!=720`$ elements, which proves that all $`P`$-orthogonal permutations are generated. Since no permutation (except the identity) is commutable with all permutations, $`S_6`$ is isomorphic to the group of transformations $`\chi _q:S_6S_6:pqpq^1`$ where $`p`$ and $`q`$ are permutations of $`6`$ elements. Such a transformation $`\chi _q`$ is completely determined by specifying the images of the 15 commutations $`p_{i,j}`$, permutations on $`\{1,2,3,4,5,6\}`$ that permute $`i`$ and $`j`$. This holds because any permutation is a composition of such commutations and $`\chi _q(p_1p_2)=\chi _q(p_1)\chi _q(p_2)`$. Note that the image under $`\chi _q`$ of a commutation is again a commutation. As a result $`S_6`$ is isomorphic to the group of permutations of $`15`$ elements obtained by restricting $`\chi _q`$ to the commutations. We will show that this is exactly the group of permutations generated by the generators $`\pi _u`$ (which can be considered as permutations of $`15`$ elements as $`0000`$ can be left out, being always mapped to itself). To this end we establish the following correspondence between nonzero $`4`$-bit vectors and commutations $`\gamma :up_{i,j}`$:

$$\begin{array}{cccc}& 0001p_{5,6},\hfill & 0010p_{4,6},\hfill & 0011p_{4,5},\hfill \\ 0100p_{2,3},\hfill & 0101p_{1,4},\hfill & 0110p_{1,5},\hfill & 0111p_{1,6},\hfill \\ 1000p_{1,3},\hfill & 1001p_{2,4},\hfill & 1010p_{2,5},\hfill & 1011p_{2,6},\hfill \\ 1100p_{1,2},\hfill & 1101p_{3,4},\hfill & 1110p_{3,5},\hfill & 1111p_{3,6}.\hfill \end{array}$$

It can be verified that $`\gamma (\pi _u(x))=\chi _{\gamma (u)}(\gamma (x))`$ for all $`u`$ and $`x`$. So $`\pi _u`$ and $`\chi _{\gamma (u)}`$ realize the same permutation of $`15`$ elements. As a consequence, also products $`\pi _{u_1}\mathrm{}\pi _{u_k}`$ realize the same permutations as products $`\chi _{\gamma (u_1)}\mathrm{}\chi _{\gamma (u_k)}`$. This finally establishes the isomorphism between $`S_6`$ and the permutations generated by the $`\pi _u`$.

It remains to be shown that there are $`6!`$ $`P`$-orthogonal $`4\times 4`$ matrices. It follows from $`A^TPA=P`$ that $`A`$ is $`P`$-orthogonal if and only if all pairs of columns of $`A`$ represent commutable $`\sigma _{a_i}`$ except for the first and second or the third and fourth column. Therefore to make an arbitrary $`P`$-orthogonal matrix, the first column $`a_1`$ can be chosen to be any nonzero 4-bit vector ($`15`$ choices), the second column should satisfy $`a_1^TPa_2=1`$ (one linear condition yielding $`8`$ possible $`a_2`$), the third column should be commutable with $`a_1`$ and $`a_2`$ (two linear conditions yielding $`3`$ choices after excluding $`0000`$) and finally the fourth column should be commutable with $`a_1`$ and $`a_2`$ and noncommutable with $`a_3`$ (three linear conditions, yielding $`2`$ possibilities). This results in $`15832=720=6!`$ possibilities. This ends the proof for $`n=2`$.

For $`n>2`$ we turn to the matrix picture and show that every $`P`$-orthogonal matrix $`A`$ can be reduced to the identity matrix by two-qubit operations, i.e. $`4\times 4`$ $`P`$-orthogonal matrices embedded in an identity matrix on rows and columns $`2k+1,2k+2,2l+1,2l+2`$ for some $`k,l\{0,\mathrm{},n1\}`$. We concentrate on two columns of $`A`$ at a time, first $`1`$ and $`2`$, then $`3`$ and $`4`$ and so on and transform them to the corresponding columns of the identity matrix with two-qubit operations. Assume, without loss of generality, that we are working on columns $`1`$ and $`2`$, then we name $`K^{(k,l)}=A_{\{2k+1,2k+2,2l+1,2l+2\},\{1,2\}}`$. If the two columns of $`K^{(k,l)}`$ are commutable they can be thought of as the first and third column of a $`4\times 4`$ $`P`$-orthogonal matrix and can be reduced by a two-qubit operation to the first and third column of an identity matrix. If the two columns of $`K^{(k,l)}`$ are noncommutable they can be reduced to the first and second column of an identity matrix. One can see that by combining such two qubit operations the first two columns of $`A`$ can be reduced to the first two columns of an identity matrix. Due to the commutability relations between the columns of $`A`$, as a result, also the first two rows become the first rows of an identity matrix. One can now proceed in a similar way with the next pairs of columns until the whole matrix is reduced to the identity matrix. The composition of the inverses of all two-qubit operations that were applied yields a decomposition of $`A`$ into two-qubit operations that can be realized by local unitary operations as shown above. This ends the proof. $`\mathrm{}`$

In the proof we saw that linear transformations ($`b=0`$) correspond to operations with $`U_B=U_A^{}`$, i.e. $`\stackrel{~}{\mathrm{\Psi }}U_A\stackrel{~}{\mathrm{\Psi }}U_A^{}`$. The matrices $`U_A`$ that under this action map tensor products of Pauli matrices to tensor products of Pauli matrices possibly with a minus-sign are known to form the Clifford group, studied in Gottesman ; Gottesman (1997) in the context of quantum error correction and quantum computation. The P-orthogonal matrices form a group that is isomorphic to a quotient group of the Clifford group. The Clifford group is known to be generated by CNOT operations and one-qubit operations that map Pauli matrices to Pauli matrices. It is possible that this knowledge may be used to give other proofs for the theorem above. However, we think that our set of generators and the isomorphism between $`P`$-orthogonal $`4\times 4`$ matrices and permutations of $`6`$ elements are worthwhile results in their own. It also follows that the CNOT operation should be decomposable in terms of our generators (at least up to phase factors, but we can do better). One can easily verify that $`\text{CNOT}=\frac{1+i}{\sqrt{2}}e^{i\frac{\pi }{4}\sigma _{1000}}e^{i\frac{\pi }{4}\sigma _{1001}}e^{i\frac{\pi }{4}\sigma _{0001}}`$. Note that the first and last operation are actually 1-qubit operations.

## III Mixed state multi-copy entanglement distillation <br>from pairs of qubits

The distillation protocols presented in this paper can be summarized as follows.

1. Start from $`n`$ identical independent Bell diagonal states with entanglement. This yields a mixture of $`4^n`$ tensor products of Bell states.
2. Apply a local permutation of these $`4^n`$ products of Bell states as described in the previous section. As a result the $`n`$ qubit pairs get statistically dependent.
3. Check whether the last $`nm`$ qubit pairs are $`|\mathrm{\Phi }`$-states ($`|\mathrm{\Phi }^+`$ or $`|\mathrm{\Phi }^{}`$). This can be accomplished locally by measuring both qubits of each pair in the $`|0,|1`$ basis, and checking whether both measurements yield the same result.
4. If all measured pairs were $`|\mathrm{\Phi }`$-states, keep the first $`m`$ pairs. This is a new mixture of $`4^m`$ products of Bell states.

This is a generalization of a protocol with $`n=2`$ and $`m=1`$, presented in Bennett et al. (1996a, b). In that protocol the applied local permutation consisted of a bilateral CNOT operation by Alice and Bob. In our protocol, we will only consider linear permutations ($`b=0`$) as one can easily see nothing can be gained by considering affine permutations. In the next section we discuss how to choose the local permutation so as to obtain a good protocol. The main result of this section is a formula for the resulting state of $`m`$ pairs as a function of the permutation of Bell states performed in step 2 of the protocol:

###### Theorem 2

If Alice and Bob apply the above protocol, starting from $`n`$ independent identical copies of a Bell diagonal state $`p_{00}|\mathrm{\Phi }^+\mathrm{\Phi }^+|+p_{01}|\mathrm{\Psi }^+\mathrm{\Psi }^+|+p_{10}|\mathrm{\Phi }^{}\mathrm{\Phi }^{}|+p_{11}|\mathrm{\Psi }^{}\mathrm{\Psi }^{}|`$ with $`p_{00}p_{01}p_{10}p_{11}`$, and with entanglement, i.e. $`p_{00}>\frac{1}{2}`$, with in step $`2`$ a local permutation, $`\pi :xAx`$ with $`A^TPA=P`$, the resulting state of the remaining $`m`$ qubit pairs is given by

$$2^{nm}\underset{y_2^{2m}}{}\left(\frac{_{x𝒮+PA^TP\overline{y}}p_x}{_{x𝒮}s_x}\right)|B_yB_y|$$
(5)
where $`𝒮`$ is the subspace spanned by the rows of $`AP`$ with indices $`2m+2,2m+4,\mathrm{},2n`$,

$$\left[\begin{array}{c}s_{00}\\ s_{01}\\ s_{10}\\ s_{11}\end{array}\right]=\left[\begin{array}{cccc}\hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 1& \hfill 1& \hfill 1& \hfill 1\\ \hfill 1& \hfill 1& \hfill 1& \hfill 1\end{array}\right]\left[\begin{array}{c}p_{00}\\ p_{01}\\ p_{10}\\ p_{11}\end{array}\right],$$

$`\overline{y}`$ is $`y`$ extended with $`2(nm)`$ zeros, and the long vector indices of $`p`$ and $`s`$ and $`B`$ behave like the indices of $`\sigma `$ in the previous section, e.g. $`p_{001101}=p_{00}p_{11}p_{01}`$.

Proof: After the permutation, and before the measurement, the state of the $`n`$ qubit pairs is given by $`_{x_2^{2n}}p_x|B_{Ax}B_{Ax}|`$. The states $`|B_{Ax}`$ with $`(Ax)_{2m+2},(Ax)_{2m+4},\mathrm{},(Ax)_{2n}=0`$ yield $`|\mathrm{\Phi }`$-states and will be kept. These are the states $`|B_{Ax}`$ for which $`x`$ is commutable with the rows $`a_{2m+2},\mathrm{},a_{2n}`$ of $`AP`$. If we call the subspace of these vectors $`x`$, $``$, the success rate (probability of keeping the first $`m`$ pairs) is $`_xp_x`$. Among the states that are kept, the ones with $`(Ax)_j=y_j,j=1,\mathrm{},2m`$ yield $`|B_y`$-states. Together with the conditions for being kept, these are $`2m+(nm)`$ independent linear conditions, yielding a coset of an $`(nm)`$-dimensional subspace of $`_2^{2n}`$. This subspace must be $`𝒮`$ since the latter is $`(nm)`$-dimensional and satisfies all homogeneous conditions (with $`y=0`$) by the $`P`$-orthogonality of $`A`$. The right coset is obtained by adding $`PA^TP\overline{y}`$ (a combination of the first $`2m`$ rows of $`AP`$, determined by $`y`$). As a result the state of the first $`m`$ pairs after the measurement is

$$\underset{y_2^{2m}}{}\left(\frac{_{x𝒮+PA^TP\overline{y}}p_x}{_xp_x}\right)|B_yB_y|$$

If all coefficients (for all $`y`$) are calculated, the denominator $`_xp_x`$ can be calculated as the sum of the $`2^{2m}`$ numerators. If only one coefficient is needed (for instance if only the fidelity of the end state is needed) the denominator can be calculated in a more efficient way as $`_xp_x=2^{(nm)}_{x𝒮}s_x`$. One can easily verify that $`s_v=_{x_2^{2n}}(1)^{v^TPx}p_x`$ (First verify for 2 bits and then extend). Therefore $`_{v𝒮}s_v=_{x_2^{2n}}(_{v𝒮}(1)^{v^TPx})p_x`$. If $`x`$ commutes with all $`v𝒮`$, $`(1)^{v^TPx}=1`$ for all $`v𝒮`$ and $`_{v𝒮}(1)^{v^TPx}=2^{nm}`$. If $`v𝒮`$, one can easily show that half the coefficients $`(1)^{v^TPx}`$ are $`1`$ and half are $`1`$. Now the states $`x`$ are exactly the ones for which $`x`$ is commutable with all elements of $`𝒮`$. Therefore $`_{x𝒮}s_x=2^{nm}_xp_x`$. This concludes the proof. $`\mathrm{}`$

## IV Recurrence schemes

With the above formula for the end state of the protocol (Theorem 2), it is possible to derive good protocols by searching over all possible values for the relevant rows of the $`P`$-orthogonal matrix $`A`$ and optimizing some quality measure. Typically this measure will depend on the fidelity of the end state and the success rate of the protocol (the probability of having $`|\mathrm{\Phi }`$-states in the measured pairs). In that case, one only needs the first coefficient (the fidelity) and the denominator (the success rate) in (5), which both only depend on $`𝒮`$, a space spanned by only $`nm`$ rows of $`𝒮`$. Although this drastically limits the search space, it still grows exponentially with growing $`n`$.

Therefore, to come up with schemes for large $`n`$, one needs to use the recurrence scheme, as was proposed for $`n=2`$ and $`m=1`$ in Bennett et al. (1996a, b). If $`m=1`$, this scheme means that the above protocol is performed $`n`$ times (with the same local permutation) and the $`n`$ identical end states are taken as the input for a new step. Of course more than two steps are possible too. One could also envision recurrence schemes with $`m1`$, for instance, combining two end states of an $`n=4`$,$`m=2`$-protocol to yield the input for a second step with $`n=4`$. In that case however the input for the second step would no longer consist of $`n`$ independent pairs. Although this only requires a minor modification of the above results ($`p_x`$ and $`s_x`$ can no longer be interpreted as products of $`p_{00},\mathrm{},p_{11}`$), we will not consider this case in this paper.

To end up with almost pure Bell states, the recurrence scheme can best be combined with the hashing protocol as in Bennett et al. (1996a, b). The hashing protocol is the best known asymptotic protocol (for $`n\mathrm{}`$) but can only be applied to Bell diagonal states with high enough fidelity. The combined protocol then consists of first applying a few recurrence steps and then switching to the hashing protocol.

The best known $`n=2`$,$`m=1`$-recurrence scheme is the one of Deutsch et al. (1996). In our language it amounts to a scheme with a $`4\times 4`$ $`P`$-orthogonal matrix whose last line is $`\mathrm{11\hspace{0.17em}11}`$. It can be proven that this scheme yields the best achievable fidelity after one step (though not achieved with the best success rate) for initial probabilities that are ordered $`p_{00}>p_{01}p_{10}p_{11}`$. For this reason it is also best to apply a pair per pair transformation after each recurrence step, which reorders the probabilities of the end state if they are not ordered. (One can easily find such one-pair transformations using the theory of section II or equivalently using the local operations of Bennett et al. (1996a, b).)

Although it is probably best to search for a new protocol for every given initial state, we propose below a protocol which we think is good if one does not have the time for such a search. We show by computer simulations that it performs better than the $`n=2`$ scheme.

Our scheme is an $`n=4`$,$`m=1`$ recurrence scheme combined with hashing, and with as the last step possibly an $`n=2`$,$`m=1`$-step if this can lead to better performance. For the local permutation (determined by the $`P`$-orthogonal matrix $`A`$) we choose a permutation that is found experimentally to often lead to the best fidelity after one step, when starting with ordered probabilities. For this reason we also apply a reordering in between recurrence steps as discussed for the $`n=2`$-protocol above. The chosen local permutation corresponds to an $`8\times 8`$ $`P`$-orthogonal matrix $`A`$ whose fourth, sixth and eighth row span the space spanned by $`\{\mathrm{10\hspace{0.17em}11\hspace{0.17em}11\hspace{0.17em}10},\mathrm{01\hspace{0.17em}10\hspace{0.17em}11\hspace{0.17em}00},\mathrm{11\hspace{0.17em}10\hspace{0.17em}10\hspace{0.17em}11}\}`$. This can be achieved by the operations

$$\begin{array}{cc}U_A=U_B^{}=\hfill & e^{i\pi /4\sigma _{\mathrm{10\hspace{0.17em}01\hspace{0.17em}00\hspace{0.17em}00}}}e^{i\pi /4\sigma _{\mathrm{01\hspace{0.17em}00\hspace{0.17em}00\hspace{0.17em}01}}}\hfill \\ & e^{i\pi /4\sigma _{\mathrm{10\hspace{0.17em}00\hspace{0.17em}11\hspace{0.17em}00}}}e^{i\pi /4\sigma _{\mathrm{00\hspace{0.17em}01\hspace{0.17em}10\hspace{0.17em}00}}}.\hfill \end{array}$$
(6)
In this realization the first and second row of the $`P`$-orthogonal matrix $`A`$ are $`\mathrm{01\hspace{0.17em}10\hspace{0.17em}00\hspace{0.17em}10}`$ and $`\mathrm{10\hspace{0.17em}10\hspace{0.17em}10\hspace{0.17em}10}`$. These rows are needed to compute the reordering operations between the steps, for although the three values of $`p_{01}^{}`$, $`p_{10}^{}`$ and $`p_{11}^{}`$ after one step of the protocol are fixed, their order is not. (The three cosets of $`𝒮`$ in $``$ in equation (5) are fixed but not their order).

This realization was found by exhaustive search over all operations that can be realized by $`4`$ consecutive elementary two-qubit operations. If, for protocols with larger $`n`$ for instance, no such simple realization can be found in a reasonable amount of time, one can always find a realization using the theory of section II but this can increase the total amount of work for the distillation protocol. This was also one of the reasons for choosing $`n=4`$ in the proposed protocol.

As a performance measure we have chosen the expected average number of input pairs needed per output Bell state in an asymptotic protocol (the inverse of the asymptotic yield). The number of recurrence steps was also chosen as to optimize this measure. Fig. 1 shows the performance for our method ($`n=4`$,$`m=1`$-recurrence with the local permutation realized by $`U_A`$ as in (6), with reordering between the steps, possibly one last $`n=2`$,$`m=1`$-step, and optimal switching to hashing protocol) and the method of Deutsch et al. (1996) with reordering between the steps and optimal switching to the Hashing protocol. For the sake of simplicity the figure only shows the results for Werner states (with $`p_{00}=F>\frac{1}{2}`$ and $`p_{01}=p_{10}=p_{11}=\frac{1F}{3}`$), but the method also performs better for non-Werner states.

To do better than this protocol for a specific initial Bell diagonal state, one can do several things depending on the amount of computing time available. One can try recurrence schemes with higher $`n`$ and even higher $`m`$, but the amount of time needed increases fast with increasing $`n`$. There is of course no obligation to take the same local permutation in consecutive recurrence steps. One can also consider distilling more than one end state at once. Making two states with two $`n=4`$,$`m=1`$-protocols is just a special case of a non-optimal $`n=8,m=2`$-protocol. One can of course search for better ones if one has the time. In this case, the two obtained Bell states will not be independent but as the fidelity goes to $`1`$, their dependence will vanish. Also two consecutive recurrence steps, say two $`n=2`$,$`m=1`$-steps, can be considered as one bigger non-optimal step, in this case with $`n=4`$,$`m=1`$. So if one has the time one can in theory always go for a one shot protocol (no recurrence), but if one combines with the recurrence scheme one can always afford lower initial entanglement with the same amount of computing time.

## V Conclusion

We have derived new protocols for distillation of entanglement from mixed states of two qubits. The protocols were based on a characterization of the group of all locally realizable permutations of the $`4^n`$ possible tensor products of $`n`$ Bell states. Our protocols perform significantly better than existing protocols as was shown by computer simulation. We also indicated how to derive even better protocols for specific initial states.

###### Acknowledgements.

Our research is supported by grants from several funding agencies and sources: Research Council KULeuven: Concerted Research Action GOA-Mefisto 666 (Mathematical Engineering), several PhD/postdoc & fellow grants; Flemish Government: Fund for Scientific Research Flanders (several PhD/postdoc grants, projects G.0256.97 (subspace), G.0240.99 (multilinear algebra), research communities ICCoS, ANMMM)); Belgian Federal Government: DWTC (IUAP IV-02 (1996-2001) and IUAP V-22 (2002-2006): Dynamical Systems and Control: Computation, Identification & Modelling); The European Commission: Esprit project: DICTAM