text:
# Orientational field-dependence of low-lying excitations in mixed state of unconventional superconductors

## Abstract

Orientational field-dependence of the zero energy density of states (ZEDOS) is calculated for superconductors with the polar state (line node), axial state (point node) and 3D d-wave state. Depending on the gap topology and relative field direction the field dependencies of ZEDOS sensitively differ, providing us a useful and practical method to identify the gap topology. It is also demonstrated that for d-wave state the field rotation in the basal plane shows a sizable oscillation ($``$3%) of ZEDOS. This is directly measurable in low-T specific heat experiment in the mixed state.

Much attention has been focused on various unconventional superconductors, ranging from high Tc cuprates, heavy Fermion materials, Boro-carbides, MgB2, to Skutterudite PrOs4Sb12. The unconventionality is associated with the gap anisotropy of the orbital function in addition to the spin structure of the Cooper pair. It is quite important to determine the detailed nodal topology of the gap function; either point or line node and their location on the Fermi surface. The determination of these characteristics is expected to lead to an understanding to the pairing mechanism of exotic superconductors sigrist .

There are several experimental methods to probe the gap anisotropy. One can basically distinguish the line and point nodes because these give rise to different and distinct power law temperature (T) dependence in various physical quantities. As for the orientation of these nodes, there is only a few ways to probe it. The field-dependent thermal conductivity $`\kappa (H)`$ and the polarization dependent sound attenuation are typical ones. By measuring $`\kappa (H,\alpha )`$ for different field direction $`\alpha `$ one can detect the location of the node in principle because the nodal quasiparticles (QP) under $`H`$ with zero-energy transport heat current. In fact, a series of experiments by Izawa et al. izawa  have determined the location of nodes in several systems. These transport measurements are, however, inevitably involved by the scattering time effect and localization effect of nodal QP, which hamper the determination of the nodal direction in some cases d-wave .

Here we propose another method based on thermodynamics: The Sommerfeld coefficient $`\gamma `$ of the $`T`$-linear specific heat at lower $`T`$ is most fundamental physical quantity in Fermionic systems of interest. Since the nodal QP created around a vortex core in the mixed state sensitively reflect its gap structure, the angle-dependent $`\gamma (H,\alpha )`$ can yield characteristic oscillation pattern relative to the nodal position under a fixed $`H`$. Recent angle-resolved $`\gamma (H,\alpha )`$ measurement on YNi2B2C by Park et al park  demonstrates a fourfold oscillation in the basal plane whose amplitude $`4\%`$, nicely coinciding with $`\kappa (H,\alpha )`$ experiment by Izawa et al. izawa . They agree with the nodal direction (), but disagree with the topology (point (line) in the latter (former)). The observed oscillation amplitude ($`4\%`$) in $`\gamma (H,\alpha )`$ is far off the theoretical prediction ($`30\%`$) based on the so-called Doppler shift argument vekhter  which is qualitative in nature (see discussion in Ref. dahm ). In this paper we calculate the zero-energy density of states (ZEDOS) in the mixed state for various situations; its direction-dependence for the polar and axial state and the angular dependence for d-wave state, in order to examine the experimental feasibility. We develop a full three dimensional (3D) computation based on quasi-classical framework which is valid for superconductors with $`k_F\xi 1`$ ($`k_F`$ Fermi wave number and $`\xi `$ the coherence length). This kind of calculations gives quantitatively reliable results.

Anisotropic pairing is routinely analyzed within the separable model of pairing potential $`V(𝒌,𝒌^{})=V_0\mathrm{\Omega }(𝒌)\mathrm{\Omega }(𝒌^{})`$. Then order parameter takes the following form: $`\mathrm{\Delta }(𝒓,𝒌)=\mathrm{\Psi }(𝒓)\mathrm{\Omega }(𝒌)`$. In the clean limit quasiclassical equations read as

$$\left[2\mathrm{}\omega _n+\mathrm{}𝒗\left(+\frac{2\pi i}{\mathrm{\Phi }_0}𝑨\right)\right]f=2\mathrm{\Psi }(𝒓)\mathrm{\Omega }(\varphi ,\theta )g,$$
(1)
$$\left[2\mathrm{}\omega _n\mathrm{}𝒗\left(\frac{2\pi i}{\mathrm{\Phi }_0}𝑨\right)\right]f^{}=2\mathrm{\Psi }^{}(𝒓)\mathrm{\Omega }(\varphi ,\theta )g.$$
(2)
Here $`\mathrm{}\omega _n=\pi T(2n+1)`$ with integer $`n`$ are Matsubara frequencies, $`𝒗`$ is Fermi velocity, $`\mathrm{\Phi }_0`$ is flux quantum, and $`f`$, $`f^{}`$, $`g`$ are Green’s functions integrated over energy normalized so that $`ff^{}+g^2=1`$. Fermi surface is assumed to be sphere. Order parameter $`\mathrm{\Psi }(𝒓)`$ and vector-potential $`𝑨(𝒓)`$ are obtained selfconsistently from the following equations

$$\mathrm{\Psi }(𝒓)\mathrm{ln}\frac{T_c}{T}=2\pi T\underset{\omega _n>0}{}\left[\frac{\mathrm{\Psi }(𝒓)}{\mathrm{}\omega _n}\mathrm{\Omega }(\varphi ,\theta )f\right],$$
(3)
$$\times \times 𝑨(𝒓)=\frac{4\pi ^2\mathrm{}N_0T}{\mathrm{\Phi }_0}\mathrm{Im}\underset{\omega _n>0}{}g𝒗.$$
(4)
Average over Fermi surface is denoted as $`\mathrm{}`$. First, polar state with a line node $`\mathrm{\Omega }(\varphi ,\theta )=\sqrt{3}\mathrm{cos}\theta `$ and axial state with point nodes $`\mathrm{\Omega }(\varphi ,\theta )=\sqrt{3/2}\mathrm{sin}\theta `$ are analyzed. Polar and azimuthal angle refer to the coordinate system with $`z`$-axis that coincide with $`c`$ crystal direction. Factors $`\sqrt{3}`$ and $`\sqrt{3/2}`$ assure that average of $`|\mathrm{\Omega }(\varphi ,\theta )|^2`$ over spherical Fermi surface is unity. We are interested in Green’s function $`g(𝒓,𝒗,\omega )`$ that describe QP excitations associated with vortices. The QP density of states $`N(E)`$ with energy $`E`$ relative to Fermi level is defined as

$$\frac{N(E)}{N_0}=\frac{N(E,𝒗)}{N_0}=\overline{\mathrm{Re}g(𝒓,𝒗,\omega 0^+iE)},$$
(5)
where $`N_0`$ is ZEDOS in the normal state. Green’s function $`g(𝒓)`$ spatially averaged over vortex lattice unit cell is denoted as $`\overline{g(𝒓)}`$. We focus on ZEDOS at low temperatures. This is because low temperature specific heat is $`C_s=\gamma (B)T=2\pi ^2\mathrm{}^2N(E=0,B)T/3`$. Therefore the equations are solved for $`T=0.1T_c`$. It is sufficient to know the Green’s function only in the vortex lattice unit cell which is divided in mesh $`41\times 41`$. Once the order parameter and vector potential are obtained selfconsistently Eqs. (1) and (2) are solved again for $`\omega 0^+`$. Typically we choose $`\omega =0.001\pi T_c/\mathrm{}`$. Method of solution is extended from ichioka  and details will be described elsewhere pedjanew . Here we present the results.

In Fig. 1 field dependence of ZEDOS for polar state with a line node is shown for $`Hc`$ (full circles) and $`Hc`$ (empty circles). Ratio $`N(E=0,B)/N_0`$ at low $`T`$ is equal to $`\gamma (B)/\gamma _N`$, where $`\gamma (B)T`$ ($`\gamma _NT`$) is low $`T`$ specific heat of superconducting (normal) phase. We discuss the power-law exponent of $`B`$ dependence of ZEDOS. It is difficult to fit the data with a single power-law function $`(B/H_{c2})^\beta `$. At least in low field we can estimate $`N(E=0,B)/N_0(B/H_{c2})^{0.35}`$ for $`Hc`$, i.e. very steep increase with field. Here is the explanation. The most important contribution to ZEDOS is coming from QP that flow in the plane perpendicular to the applied field. For $`Hc`$ geometry, those QP experience zero energy gap. They are easily excited and extended outside of the vortex core even in low field in comparison to $`s`$-wave superconductors. The outcome is steep increase of ZEDOS with field. Experimentally similar small exponent is observed in MgB2. The physics is analogous to the case of polar state. Small exponent is coming from the small gap at the $`\pi `$-band nakai  in MgB2 while coming from the line node in the polar state. For perpendicular orientation $`Hc`$ (empty circles) the problem is analogous to that of two dimensional (2D) d-wave case for fields along $`c`$-axis. Power law with exponent $`\beta 0.45`$ calculated here should be compared with self-consistent calculation on cylindrical Fermi surface and 2D $`d`$-wave gap function which reveals $`N(E=0,B)/N_0(B/H_{c2})^{0.43}`$ power law ichioka . For this geometry $`Hc`$, QP in plane $`H`$ experience zero gap only if their momentum is in basal plane. Therefore, they are more difficult to excite compared to parallel geometry $`Hc`$, hence the exponent $`\beta `$ is bigger.

It is important to emphasize that ZEDOS in Fig. 1a) is plotted against $`B/H_{c2}`$, where $`H_{c2}`$ is different for each geometry. For polar state there is a large anisotropy of upper critical field $`H_{c2}^{}2H_{c2}^{}`$. When plotted versus magnetic field, Fig. 1b), ZEDOS lines crosses at some critical field $`B_{cr}`$. Therefore by rotating magnetic field from $`c`$-axis toward the basal plane ZEDOS may increase or decrease depending on field value.

For axial state with two point nodes and $`Hc`$, energy gap is small only for small fraction of QP that flow along the $`c`$-axis, which makes small contribution to total ZEDOS. In this geometry ZEDOS resembles that in $`s`$-wave superconductors. Most of the low energy QP are trapped at vortex cores, at least in low field, thus $`N(E=0,B)/N_0B/H_{c2}`$. This is confirmed by numerical calculation shown in Fig. 2a) (full circles). For field $`Hc`$ (empty circles) power-law exponent is smaller than that in line node polar state for both field geometries, Fig 1a). Roughly speaking, the larger is the angular area of suppressed gap, the faster ZEDOS is increasing with induction $`B`$. Note that the field dependence $`N(E=0,B)/N_0(B/H_{c2})\mathrm{ln}(B/H_{c2})`$ predicted by the Doppler shift calculation for point node case mineev  is far off the present result in Fig. 2, warning us its validity. Similar to the polar state, ZEDOS curves for two field directions cross at some critical field $`B_{cr}`$ as shown in Fig. 2b).

We show the importance of the mutual arrangement of line node and magnetic field on ZEDOS. In this sense it is interesting to examine 3D version of $`k_x^2k_y^2`$ symmetry of the gap function which is given by: $`\mathrm{\Omega }(\varphi ,\theta )=\sqrt{15/4}\mathrm{sin}^2\theta \mathrm{cos}2\varphi `$. The form of the gap function is a natural choice for spherical Fermi surface. The question is, can we guess for which field direction in the basal plane ZEDOS is maximum based on the calculation for polar state? Namely, if the magnetic field $`Hc`$ is along the node, then all QP flowing perpendicular to the field experience zero energy gap (see inset in Fig. 3) analogous to polar state with $`Hc`$. If the field is along the antinode direction then QP that flow perpendicular to the field experience zero gap only if their momenta are along the c-axis, analogous to polar state and $`Hc`$. These simple qualitative arguments, suggesting $`N(E=0,antinode)<N(E=0,node)`$, are misleading since our calculation gives opposite result. It is because one must take into account the contribution from QP that are flowing at some angle with respect to the field in this 3D problem. It is instructive to see how angle resolved ZEDOS $`N(E=0,𝒗)`$, averaged over angle $`\varphi ^{}`$ in plane perpendicular to field, changes with angle $`\theta ^{}`$ between QP velocity $`𝒗`$ and vortex axis. We plot this quantity, multiplied with weighting factor $`\mathrm{sin}\theta ^{}`$ in Fig. 3 at some very low field. Total ZEDOS is area under the curve. Two field directions perpendicular to the c-axis are considered, node and antinode. For QP that flow perpendicular to the magnetic field, $`\theta ^{}=90^{}`$, ZEDOS for node field direction is about 3 times larger than for antinode field direction. However for antinode direction QP experience gap node as long as the angle between QP velocity and field direction is $`45^{}<\theta ^{}<90^{}`$. In spite of the $`\mathrm{sin}\theta ^{}`$ weighting factor there is a significant contribution to the total ZEDOS that are coming from these QP. As a result at low field $`N(E=0,antinode)`$ is larger than $`N(E=0,node)`$ by a few percents.

In Fig. 4a) field dependence of ZEDOS for node and antinode field direction is shown. Exponent $`\beta 0.32`$ for antinode direction differs from $`\beta 0.45`$ of polar state ($`Hc`$) in Fig. 1. This difference comes from the different power expansion of the gap function in the node vicinity. In the antinode case, in plane perpendicular to the field, gap function can be approximated as $`|\mathrm{\Omega }(\varphi ,\theta )|\mathrm{\Omega }_0\theta ^2`$ in the vicinity of gap node $`\theta =0`$. On the other hand, $`|\mathrm{\Omega }(\varphi ,\theta )||\theta \pm \pi /2|`$ near node $`\theta =\pm \pi /2`$ for polar state and $`Hc`$. The latter case is analogous to 2D d-wave function and $`Hc`$ giving $`\beta 0.45`$. It was shown by Barash and Svidzinsky barash  that temperature dependence of specific heat is closely related to the exponent $`n`$ of gap function power expansion near node. The larger is exponent $`n`$ the faster specific heat $`C_s`$ is increasing with $`T/T_c`$. Similar qualitative arguments can be applied to field dependence of specific heat.

At $`B=B_{cr}`$ two ZEDOS curves cross, the same as for polar and axial state, see Fig. 4b). In the inset of Fig. 4b), ratio $`R=N(antinode)/N(node)`$ is plotted against induction $`B`$. Physics of the crossing is very simple and will be explained on the 3D d-wave case. In the present model Fermi surface is assumed to be sphere and upper critical field anisotropy is determined by the gap function. For our simple 3D d-wave case $`H_{c2}^{node}<H_{c2}^{antinode}`$. Therefore, for a fixed high field $`N(E=0,node)>N(E=0,antinode)`$ because along the node field direction the superconductor is closer to the normal state and ZEDOS is closer to the normal state value $`N_0`$. On the other hand, for $`BH_{c2}`$ the QP excitations probe the gap structure since the biggest contribution to the ZEDOS is coming from the delocalized QP. It was calculated $`N(E=0,node)<N(E=0,antinode)`$.

The value of crossing field $`B_{cr}`$ depends on the Fermi surface model. Upper critical field is also affected by the Fermi surface anisotropy, and can reverse the sign of four-fold $`H_{c2}`$ oscillations in the basal plane. For example, in YNi2B2C the gap node is along $`[100]`$, $`[010]`$ directions, implying that those are also the directions of $`H_{c2}`$ minima. But in borocarbides Fermi surface is highly anisotropic. If we accept that LuNi2B2C has similar electronic structure as Y-borocarbide, then minimum of upper critical field along $`[110]`$ direction metlushko  implies the decisive role of Fermi surface on $`H_{c2}`$ anisotropy. Thus, in this case the high field inequality should be $`N(E=0,node)<N(E=0,antinode)`$. Since we expect that Fermi surface anisotropy has no role in low-field ZEDOS, then the sign of four-fold ZEDOS oscillation should be the same for all $`B`$, i.e. there is no crossing of two ZEDOS lines.

In Fig. 5 low field angular dependence of ZEDOS for field rotating in the basal plane is shown. Four-fold oscillations is what one expects from the symmetry of the gap function. Angular variation is $`3\%`$ at low fields, which is measurable with present experimental techniques park . Note that in the 2D d-wave case Doppler-shift calculation vekhter  estimates angular variation as large as $`30\%`$. Parabolic-like minimum in angular dependence of ZEDOS is in contrast with cusp-like minimum in 2D d-wave case (and cylindrical Fermi surface) vekhter ; hirschfeld ; schachinger . Cusp-like minimum in thermal conductivity angular dependence is predicted in the $`s+g`$ model (point node) and observed experimentally izawa . While shape of ZEDOS minimum can rule out some forms of the gap function it can not provide the unique answer on the question of node topology (line or point). It is necessary to tilt the field out of the basal plane and study ZEDOS to gain additional information. This was done by measuring thermal conductivity izawa  and specific heat park  with fields rotating around $`c`$-axis in YNi2B2C.

We have studied the orientational field-dependence of the nodal QP with zero-energy in the mixed state for the three representative gap functions, namely the axial, polar and d$`_{x^2y^2}`$ states. Our computation is based on quasi-classical approach for 3D Fermi sphere. We have demonstrated that the orientational dependent and angle-resolved specific heat measurements are an ideal tool to distinguish line and point nodes and to locate the nodal direction free from scattering time or localization effects associated with transport experiments and also that this can be feasible in the present-day technical limitations. When conducting field-rotation experiment, it is important to keep the field low ($`B<B_{cr}`$) to probe the intrinsic gap structure.

We thank Y. Matsuda, T. Sakakibara, K. Izawa, I. Vekhter, M. B. Salamon and T. Park for useful discussions.


format: text

text:
# The Synoptic Swift Synergy – Catching Gamma-Ray Bursts Before They Fly

## 1 Introduction

The gamma-ray burst GRB 021004 was detected by HETE II at 12:06 UT on the 4th of October 2002 (Shirasaki et al. 2002). Observations after about 9 minutes from the trigger revealed a fading optical transient (Fox 2002), which was densely sampled in several bands, especially at early times. The afterglow of GRB 021004 has shown several unusual features (Mirabal et al. 2002; Salamanca et al. 2002; Moller et al. 2002; Bersier et al. 2003). Perhaps its most unique feature was that the field had been observed shortly before the gamma-ray burst itself was detected. Wood-Vasey et al. (2002) give a limiting unfiltered magnitude of 21.4 for observations on the day before the burst and 22.3 integrated over the year before the burst.

Astronomy is on the threshold of a new era where large portions of the sky are surveyed deeply and regularly. The question arises what is likelihood of getting photometry of a gamma-ray burst precursor, specifically if supernovae precede gamma-ray bursts as in the supranova model (Vietri & Stella 1998)? Although the flux upper limits for GRB 021004 are not strigent enough to constain theoretical models of gamma-ray bursts, the high-burst localization rate of Swift combined with the fast sky coverage of the Sloan Digital Sky Survey (SDSS) and later Pan-STARRS, LSST and the Supernova-Acceleration Probe (SNAP) could provide important constraints on gamma-ray burst precursors.

During the first and second years of operation of the Swift mission, SDSS will scan approximately 3000 square degrees (or 7% of the sky) each year (SDSS Collaboration 2001). Over this area it will detect point sources down to $`R23.2`$. If Swift or subsequent missions are operational in 2006, the Pan-STARRS program will observe 20,000 square degrees every four days (or 50% of the sky) to a limiting magnitude of $`R24.2`$ (Kaiser et al. 2002). Finally, potentially beginning in 2010, SNAP will cover 15 square degrees every four days with each observation reaching a limiting magnitude of $`R28`$ (Kim et al. 2002), coadding observations over a month would go one magnitude deeper. The SNAP lensing survey will cover 300 square degrees over five months to a similar limiting magnitude.

Long gamma-ray bursts are thought to be associated with the collapse of a massive star, a supernova. Specifically, in the collapsar model, the formation of a black hole in the center of the star results in relativistic jets which pierce the envelope of the star (MacFadyen & Woosley 1999). Along the axis of the jets, the collapsing star appears as a gamma-ray burst, and the supernova reaches its peak a few weeks after the GRB. Vietri & Stella (1998) proposed an alternative model in which the gamma-ray burst accompanies the delayed collapse of a quickly spinning neutron star which is more massive than the maximum mass of a non-rotating neutron star. The neutron star may take several months or years after the supernova to spin down to the critical frequency and collapse.

In this Letter, I will estimate the number of gamma-ray burst events with photometry which overlaps on the sky but shortly precedes in epoch from SDSS and other surveys and compare the flux limits with the expected flux from a supernova which may precede the gamma-ray burst.

## 2 Gamma-Ray Burst Overlap with Future Surveys

To calculate how often sufficiently deep photometry will precede the observation of a gamma-ray burst on the sky, several ingredients are required: a model for the spectral-energy distribution as a function of time of a supernova associated with a GRB, an estimate of the luminosity-rate function of GRBs as a function of redshift ($`\dot{\varphi }(z,L)`$), a model for the field of view of the gamma-ray burst detector ($`\mathrm{\Omega }_{\text{GRB}}=2`$ for Swift) and its detection threshold ($`P_1`$) and the rate of sky coverage of the photometric program ($`_{\text{photo}}`$) and its detection threshold ($`R_{\text{lim}}`$). Porciani & Madau (2001) provide models for $`\dot{\varphi }(z,L)R_{\text{GRB}}(z)\psi (L)`$. The rate of GRBs, $`R_{\text{GRB}}(z)`$, is taken to be proportional to the star-formation rate, and the luminosity function of GRBs, $`\psi (L)`$, is constrained by the BATSE GRB number counts. The rate of overlapping photometry is given by the product of the rate of sky coverage with a integral over the assumed cosmological distribution of GRBs,

$`{\displaystyle \frac{\mathrm{d}N_{\text{overlap}}}{\mathrm{d}t}}(P>P_1,R<R_{\text{lim}})`$ $`=`$ $`{\displaystyle \frac{_{\text{photo}}\mathrm{\Delta }t}{4\pi }}{\displaystyle \frac{\mathrm{d}N_{\text{total}}}{\mathrm{d}t}}(P>P_1,R<R_{\text{lim}})`$ (1)
$`{\displaystyle \frac{\mathrm{d}N_{\text{total}}}{\mathrm{d}t}}(P>P_1,R<R_{\text{lim}})`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_{\text{GRB}}}{4\pi }}{\displaystyle _0^{z:R(z)=R_{\text{lim}}}}dz{\displaystyle _{L(P_1,z)}^{\mathrm{}}}dL{\displaystyle \frac{\mathrm{d}V(z)}{\mathrm{d}z}}{\displaystyle \frac{\dot{\varphi }(z,L)}{1+z}}.`$ (2)
For lack of a better model for the evolution of a supernova associated with a GRB, I will assume that SN2001ke (Garnavich et al. 2003) is a prototype for this class, and furthermore that a supernovae associated with a GRB maintains its peak brightness for a period $`\mathrm{\Delta }t=14(1+z)`$ days in the observer’s frame and otherwise it is undetectable (see Reichart 1999; Bloom et al. 1999, for other GRB-associated supernovae). It is reasonable to use the median value of $`z`$ for GRBs whose associated supernova are brighter than the magnitude limit of the particular photometric survey. However, to be highly conservative, I will take $`\mathrm{\Delta }t=14`$ days to calculate the rate of overlap.

If a survey covers the same area of sky more often than once per interval $`\mathrm{\Delta }t`$ such as the SNAP supernovae search and Pan-STARRS, the rate of sky coverage $`_{\text{photo}}`$ should only account for the first visit in each period $`\mathrm{\Delta }t`$; for example $`_{\text{photo}}`$ for Pan-STARRS is $`2\pi `$ per fourteen days. The additional visits during each fortnight do not increase $`_{\text{photo}}`$ but they do allow the survey to probe deeper by coadding the successive images.

According to the original supranova model (Vietri & Stella 1998), the supernova may reach its peak at any time up to several years before the GRB, so this calculation implicitly assumes that both the GRB survey and the photometric survey will be operating at the appropriate times. $`L(P_1,z)`$ is the luminosity of a GRB at a redshift $`z`$ which is detected at a count-rate of $`P_1`$, and $`R(z)`$ is the $`R`$-band apparent magnitude of a GRB-associated supernovae at a redshift $`z`$. Both of these functions include the $`k`$correction (Hogg 1999) and assume the cosmographic parameters, $`\mathrm{\Omega }_M=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. Fig. 1 plots the number of GRBs detected by Swift per year whose associated supernova would be brighter at its peak than a particular $`R`$-band magnitude.

The results shown in Fig. 1 assume the SF1 model of Porciani & Madau (2001). This model provides a conservative lower limit for the overlap. It predicts that Swift will localize about 110 bursts per year — the more generous estimates range up to 300 bursts per year (Myers 2002). Furthermore, this model predicts that the bursts detected will be a higher redshifts than other models, so the accompanying supernovae will be fainter and more difficult to detect.

Table 1 gives the overlap rate between various photometric surveys and the Swift GRB localization mission. What is striking is that the shallow but wide Pan-STARRS and LSST surveys will perform much better than any of the other surveys. Furthermore, if supernovae precede GRBs, Pan-STARRS and LSST each will detect nearly ten GRB-associated supernovae per year. If it finds none, it would place severe constraints on the supranova model for GRBs. It must be emphasized that this rate of overlap is extremely conservative. It assumes a low Swift burst localization rate and a distribution of GRBs skewed to high redshift (therefore, faint assocated supernovae). The actual rate of overlap will probably be higher if both programs operate simultaneously. Furthermore, Swift will generate a catalog of burst positions and redshifts. One will be able to cross-correlate a posteriori this catalog with earlier Pan-STARRS or LSST observations and exclude the appearance of transients to $`R25`$ over a wide range of epochs preceding the burst yielding definitive constraints on GRB-progenitors independent of assumptions about the GRB luminosity function and its evolution.

This calculation of the overlap rate assumes that either the GRB localization program or the photometric survey studies random portions of the sky. In fact both the Swift mission and all of the photometric surveys avoid studying the region of the sky near the sun. Although the average rate of overlap over a year in given by the formulae above and the values in the tables, the chance of detecting the supernova associated with GRB is somewhat higher than average if the supernova precedes the GRB by less than three months or between nine and fifteen months. If the supernova precedes the GRB by six to nine months, the chance of detecting it is somewhat lower than average. However, this seasonal variation is smaller than the uncertainities in the GRB luminosity function.

## 3 Discussion

The philosophy employed for finding gamma-ray burst precursors is somewhat different that what is necessary for finding supernovae or microlensing events. Because the precursors will be sought after the gamma-ray burst is detected and localized, it is not necessary to have more than one epoch of data from the particular region of sky before the burst. Even a single epoch would yield important constraints. Furthermore, unless the cadence of the observations is sufficiently low (no more than biweekly), the repeated observations of the same patch of sky do not improve the chances of catching a precursor (unless one coadds the data to probe deeper), because supernova typically evolve over the course of weeks. Consequently, although supernova and microlensing surveys have a large data rate of high quality photometry, because of their relative lack of sky coverage and depth they do not contribute much to the detection rate of precursors. From another point of view, only a small fraction ($`<10^4`$) of supernovae result in GRBs directed toward us, so one would typically have to find at least $`10^4(4\pi /\mathrm{\Omega }_{\text{GRB}})`$ supernovae in a blind search before finding a single GRB-associated supernova.

The best bets are the large deep wide surveys of the sky. SDSS is the prototype and Pan-STARRS and LSST should deliver results. There is a small possibility that SDSS will catch a supernova before a gamma-ray burst providing important evidence for the supranova model for gamma-ray bursts (it may have done so already). Pan-STARRS or LSST if it overlaps with a high-locatization-rate GRB mission such as Swift will be able to provide important constraints on gamma-ray-burst models. Specifically, it will be able to exclude the possibility that GRBs follow supernovae within a year.

I would like to acknowledge useful discussions with Kris Stanek, Robert Lupton and Bob Kirshner and helpful suggestions from the anonymous referee. I was supported by the Chandra Postdoctoral Fellowship Award # PF0-10015 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073.


text:
# Instabilities and Bifurcations of Nonlinear Impurity Modes

## I Introduction

Nonlinear localized modes in discrete systems, also called intrinsic localized modes or discrete breathers, are self-localized states that exist due to the interplay between lattice coupling and nonlinear effects. Spatially localized modes of the discrete nonlinear Schrödinger (DNLS) equation, known as discrete solitons, have appeared in many diverse areas of physics, such as biophysics, nonlinear optics, solid state physics  and, more recently, in the studies of the Bose-Einstein condensates in optical lattices  and photonic-crystal waveguides and circuits .

In application to nonlinear guided-wave optics, discrete solitons have been suggested to exist in nonlinear waveguide arrays , and they have been observed a decade later in AlGaAs waveguide structures  (see also a review paper ). Several other issues, such as the discreteness-induced soliton dynamics and soliton interaction, have been investigated. Regarding applications, it was recently shown that discrete solitons in two-dimensional networks of nonlinear waveguide arrays can be used to realize various functional operations, such as blocking, routing, logic functions, time gating, etc. .

On the other hand, localized impurities are known to play a crucial role in numerous physical systems, not only introducing interesting wave scattering phenomena , but also creating the possibility for the excitation of impurity modes, which are spatially localized oscillatory states at the impurity sites . The relevant phenomenology has been recognized as important in a variety of physical settings ranging from defect modes in superconductors  to the dynamics of the electron-phonon interactions  and from the propagation of light in dielectric super-lattices with embedded defect layers  to defect modes in photonic crystals .

It is crucially important to examine the interplay between disorder and nonlinearity. For the continuum problems, this interplay is known to lead to the existence of symmetric and asymmetric impurity modes and their interesting stability properties (see Refs.  and references therein). Here we concentrate on the study of the discrete systems such as the DNLS model. Some earlier results for the discrete lattices focused on the scattering of localized modes by impurities in the DNLS equation , as well as in the Klein-Gordon model . The effect of the impurity on localization properties of the nonlinear lattice and stability of localized modes was studied, to the best of our knowledge, only for lattices with anharmonic coupling .

However, the earlier theoretical predictions  and recent experimental results  for the interaction of the discrete solitons with structural defects in arrays of AlGaAs optical waveguides suggest that discreteness should play a crucial role in the properties and stability of nonlinear localized modes. In this paper, we study, in the framework of the DNLS model, how the interaction of a discrete localized mode with a repulsive nonlinear impurity leads to different localized states, as well as examine, both analytically and numerically, the stability criteria for the bifurcations between different localized modes near the impurity site.

The existence of multiple localized states near the impurity site can be understood by means of simple physics. Indeed, in the continuum approximation, a two-hump localized mode centered at the defect  is known to be unstable with respect to an exponential growth of antisymmetric linear perturbations . The above-mentioned instability leads to the motion of the localized mode away from the defect , and for small distances $`\xi `$ between the mode center and the defect, the effective interaction energy can be presented in the form $`H_{\mathrm{int}}H_0A\xi ^2`$, where $`A`$ is defined by the defect parameters. However, this situation becomes quite different in discrete systems because of the so-called “Peierls relief”, an effective periodic potential due to the lattice discreteness. Such potential was first calculated for small-amplitude breathers in a nonlinear elastic chain described by the discrete nonlinear Klein-Gordon equation  applying an approach suggested for kinks of the Frenkel-Kontorova model . The potential has the following structure, $`H_pH_1B\mathrm{cos}(2\pi \xi /h)`$, where $`h`$ is the lattice spacing. Thus, a competition between these two potentials defines a sequence of stable positions of a localized mode near the defect. This qualitative picture is valid in the limit of small $`h`$ and, generally speaking, the problem should be solved numerically. Below, we consider this problem in the framework of the DNLS model with a single on-site nonlinear defect.

The paper is organized as follows. In Sec. II we introduce our model and discuss the analytical results obtained for the family of the stationary states localized near the defect site. Then, in Sec. III we study numerically symmetric localized modes and their symmetry-breaking instability, including the analysis of the linear eigenvalue problem. Section IV includes the results of the instability-induced dynamics of localized modes and the study of the families of asymmetric localized modes which appear as stable states due to a balance of the defect repulsion and the effective trapping potential of the lattice. We also find the families of asymmetric modes and demonstrate their stability. Finally, Sec. V concludes the paper.

## II Model and Analytical Results

We consider a DNLS model with a nonlinear impurity on one site of the lattice. To examine the interplay between discreteness, nonlinearity and disorder, we study the case of a focusing nonlinearity in the waveguide array with a repulsive nonlinear defect. The DNLS equation describing our model can be written in the following dimensionless form,

$`i\dot{u}_n+{\displaystyle \frac{C}{2}}\mathrm{\Delta }_2u_n+|u_n|^2u_n=\alpha \delta _{n,n_0}|u_n|^4u_n,`$ (1)
where $`\mathrm{\Delta }_2u_n(u_{n+1}+u_{n1}2u_n)`$, $`u_n`$ describes the complex envelope of the electric field in the $`n`$-th waveguide, while the dot stands for the spatial derivative along the array propagation direction (which we treat here as time $`t`$).

In this study, we select the quintic nonlinearity for the defect, in order to differentiate its effect from the lattice as well as to model non-Kerr nonlinear response often observed in experiment (see, e.g. Ref. ). However, we expect that the results will be qualitatively valid for any type of the power-law nonlinearity of the order $`n`$. The generalized nonlinearity of the defect can also appear in the theory of photonic circuits .

The coupling coefficient $`C`$ can be expressed through the lattice spacing $`h`$ as follows $`C=1/h^2`$, and it characterizes the effective diffraction in the array. We should notice that the lattice spacing $`h`$ may not coincide with the relative distance between the waveguide in an array, and these values are connected in a more complicated way .

Equation (1) stems from the Hamiltonian

$`H={\displaystyle \underset{n}{}}{\displaystyle \frac{1}{2}}\left(C|u_{n+1}u_n|^2|u_n|^4\right)+{\displaystyle \frac{\alpha }{3}}|u_{n_0}|^6,`$ (2)
upon setting its derivative with respect to $`u_n^{}`$ equal to $`i\dot{u}_n`$.

The linear waves of the form $`u_n=u_0\mathrm{exp}(i\mathrm{\Omega }tiknh)`$ have the (phonon) spectrum $`\mathrm{\Omega }=(2/h^2)\mathrm{sin}^2(kh/2)`$ consisting of a band of negative frequencies ($`\mathrm{\Omega }<0`$).

In the homogeneous (i.e., no impurity or equivalently $`\alpha =0`$), continuum (i.e., $`h0`$) limit, the model (1) has an exact localized solution of the form (in the long-wave approximation, i.e. $`h\eta 1`$)

$`u(x,t)={\displaystyle \frac{\eta e^{i\mathrm{\Lambda }t}}{\mathrm{cosh}[\eta (x\xi )]}},`$ (3)
where $`\mathrm{\Lambda }=\eta ^2/2`$, $`\mathrm{\Lambda }`$ is the frequency of the breather (and should be positive to avoid resonances with the linear spectrum of plane waves), and $`\eta `$ essentially defines an effective mass of the breather excitation in the long wave approximation according to $`M_{eff}=2\eta `$.

A natural approach to the study of a discrete model is to use the expression (3) as an ansatz solution in the Hamiltonian (2) in order to obtain information about the static solution of Eq. (1) and its stability. The last term of Eq. (2) is then easy to evaluate in the framework of the ansatz (3). On the other hand, the first two terms (in the parenthesis) of Eq. (2) can be evaluated using the Poisson summation formula  which states that

$`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}f(\beta n)={\displaystyle \frac{\sqrt{2\pi }}{\beta }}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}F\left({\displaystyle \frac{2m\pi }{\beta }}\right),`$ (4)
where $`F`$ is the Fourier transform of $`f`$,

$`F(k)={\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}f(x)e^{ikx}𝑑x.`$ (5)
Using the above remarks and formulas for the Hamiltonian (2), we can evaluate $`H=H(\xi ;\{\eta ,h\})`$ (up to constant, $`\xi `$-independent terms) as the following

$$H=H_p+H_{\mathrm{int}},$$
(6)
where

$`H_\mathrm{p}={\displaystyle \frac{4\pi ^2}{h^2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{m\mathrm{cos}(2m\pi \xi /h)}{\mathrm{sinh}(m\pi ^2/\eta h)}}\left[{\displaystyle \frac{1}{h^2}}+{\displaystyle \frac{\eta ^2}{3}}\left(1+{\displaystyle \frac{m^2\pi ^2}{\eta ^2h^2}}\right)\right],`$

$`H_{\mathrm{int}}={\displaystyle \frac{\alpha \eta ^6}{3\mathrm{cosh}^6(\eta \xi )}}.`$

Notice that the first part of the expression pertaining to the summation assumes the correct asymptotic (i.e., small $`h`$) form used e.g., in Eq. (3.5) in , as well as in references therein. Finally, it is worth mentioning that the series in Eq. (6) has terms which are exponentially smaller with respect to the leading order terms and hence even keeping a small number of terms in the series should yield reasonably accurate results. In the leading order approximation, the Hamiltonian can be approximated as follows

$`H{\displaystyle \frac{16\pi ^2}{h^4}}e^{\pi ^2/h\eta }\mathrm{cos}(2\pi \xi /h)+{\displaystyle \frac{\alpha \eta ^6}{3\mathrm{cosh}^6\eta \xi }}`$ (7)
The criterion for obtaining a steady state in the reduced variable formalism of $`H=H(\xi ;\{\eta ,h\})`$ requires that

$`{\displaystyle \frac{H}{\xi }}|_{\xi =0}=0.`$ (8)
It is worth noting that the fact that the solutions of Eq. (8) will be static solutions of Eq. (1) has been rigorously proved recently by Kapitula , through the use of a Lyapunov-Schmidt reduction. Given that the ansatz of Eq. (3) is not exact, this result should be interpreted in the orbital sense (i.e., there are orbits in the configuration space which are close to the approximate solution of Eq. (8) in conjunction with Eq. (3), in the appropriate norm). In particular, in this case one can directly observe that $`\xi =0`$ will always be a solution of Eq. (8). Additional solutions are present close to $`\xi =nh`$ ($`n𝒵`$), and unstable ones close to $`\xi =(n+1/2)h`$. In the homogeneous limit, the pulse-like steady states are exactly at these points, as is well-known, but now the presence of the impurity modifies their exact location, albeit slightly due to the exponential decay of the impurity potential away from its location at $`\xi =0`$.

Of particular interest is the stability of the solution at $`\xi =0`$. One naturally expects that when the impurity is sufficiently strong and repulsive, it will reverse the local stability picture close to $`\xi =0`$, converting the relevant local energy minimum into a local energy maximum. We should note here that as per the positive definite nature of the Hessian matrix in the focusing case (for a single pulse), the local minima of the energy landscape $`H=H(\xi ;\{\eta ,h\})`$ will correspond to stable solutions. The opposite will be true for local maxima. The Hessian matrix here is given by $`P/\mathrm{\Lambda }`$, where $`P`$ is the power or (squared) $`L^2`$ norm of the solution $`P=u_2^2`$ . Recall that the $`L^p`$ norm is defined as $`u_p=(_n|u_n|^p)^{1/p}`$. The $`L^{\mathrm{}}`$ norm (used also below) is $`u_{\mathrm{}}=\mathrm{max}_n|u_n|`$. $`P`$ can be thought of as the number of elementary excitations present in the field. In the optics literature, this is often called the “energy” of the solution.

Hence, the instability criterion in the reduced, one degree of freedom dynamical system for the relevant steady state will be

$`{\displaystyle \frac{^2H}{\xi ^2}}|_{\xi =0}<0,`$ (9)
which can be restated in the form of the asymptotic approximation, through the leading order term of the series in Eq. (6), as follows $`\alpha >\alpha _c`$, where

$`\alpha _c={\displaystyle \frac{32\pi ^4}{h^6\eta ^8}}e^{\frac{\pi ^2}{\eta h}}.`$ (10)
In terms of the full series, the critical value reads

$`\alpha _c={\displaystyle \frac{8\pi ^4}{\eta ^8h^6}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{m^3}{\mathrm{sinh}(m\pi ^2/\eta h)}}\left[{\displaystyle \frac{1}{h^2}}+{\displaystyle \frac{\eta ^2}{3}}\left(1+{\displaystyle \frac{m^2\pi ^2}{\eta ^2h^2}}\right)\right],`$ (11)
which can be easily confirmed to have as leading order behavior Eq. (10).

The result (11) describes the stability threshold of the discrete mode localized at the impurity site $`n=0`$, so that when the impurity is sufficiently strong and repulsive, it will reverse the local stability of the mode close to $`\xi =0`$, converting the relevant local energy minimum into a local energy maximum. However, the results (8), (9) are more general, and they suggest the existence of many different stable stationary states located at certain distances from the impurity according to the local minima of the effective energy, the minima are produced by the effective periodic potential due to the model discreteness. For all such states, the physical picture based on the reduced Hamiltonian is valid only approximately, and more detailed numerical studies are required.

## III Localized modes and their stability

We now turn to numerical results. We numerically constructed solutions of Eq.(1), for the chain with $`N=200`$ sites, localized at the impurity site $`n_0=100`$, using a Newton method with an initial condition given by Eq. (3), assuming $`\mathrm{\Lambda }=0.5`$ ($`\eta =1`$). Upon convergence, linear stability analysis of the resulting pulses was performed. This was done by linearizing using the ansatz: $`u_n=u_{sol}+ϵ\mathrm{exp}(i\mathrm{\Lambda }t)[a_n\mathrm{exp}(i\omega t)+b_n\mathrm{exp}(i\omega ^{}t)]`$; $`u_{sol}`$ is the exact solution (around which we are linearizing) and the $`\omega `$’s denote the linearization eigenfrequencies. Then to O($`ϵ`$) we obtain the relevant (linearization) eigenvalue problem

$`\omega \left(\begin{array}{c}a_k\\ b_k^{}\end{array}\right)=𝐉\left(\begin{array}{c}a_k\\ b_k^{}\end{array}\right),`$

where $`𝐉`$ is the linear stability (Jacobian) matrix of the form

$`𝐉=\left(\begin{array}{cc}\frac{F_i}{u_j}& \frac{F_i}{u_j^{}}\\ \frac{F_i^{}}{u_j}& \frac{F_i^{}}{u_j^{}}\end{array}\right),`$

where $`F_i=\frac{1}{2}C(u_{i+1}+u_{i1}2u_i)u_i^2u_i^{}+\alpha u_i^3u_{i}^{}{}_{}{}^{2}\delta _{i,n_0}`$. When the corresponding eigenfrequencies $`\omega `$ are real, then the solution is linearly stable. On the contrary, the presence of a mode with a nonzero imaginary part in its eigenfrequency denotes the presence of an instability.

Figure 1 shows a typical mode profile obtained in different cases. These soliton profiles are similar to those discussed in the long-wave approximation in Refs. . For $`\alpha =0`$ (homogeneous case), the circles denote the discretization of the continuum profile in accordance with Eq. (3). The corresponding exact discrete homogeneous solution is shown by the stars. These results are shown for the lattice spacing parameter $`h=1`$. Notice, that discreteness induces a narrowing and increase of amplitude with respect to the continuum pulse. Finally, the pluses indicate the solution at the onset of the instability (that occurs as $`\alpha `$ is increased), for $`\alpha =0.22`$. An interesting observation is that at the onset of the instability, the pulse has essentially regained its continuum shape but for a significantly larger amplitude at the central site.

A more detailed study of the mode instability that occurs as $`\alpha `$ is increased is presented in Figs. 2 and 3. Figure 2 shows three different examples of the discrete localized mode and the corresponding linear stability results. The top panel is for $`\alpha =0.1`$, in the subcritical case, where the stability of the mode can be determined by the absence of imaginary eigenfrequencies (see the top right panel in Fig. 2). The middle row shows a slightly supercritical value of the impurity strength ($`\alpha =0.3`$, while $`\alpha _c0.22`$), while a strongly unstable, supercritical mode for $`\alpha =0.5`$ is shown in the bottom panel. It is important to remark here that a perhaps counter-intuitive result of our numerical investigations is that the instability point is not coincident with the point (in parameter space) where the symmetric mode becomes two-humped. We have generically observed that the mode becomes two-humped for considerably larger (than the instability threshold) values of the defect strength $`\alpha `$.

Similar calculations were performed for different values of the discreteness parameter $`h`$ and a two-parameter stability diagram is shown in Fig. 3. Above the different curves the (centered at $`\xi =0`$) solutions are unstable, while the opposite is true below the curves. The solid line denotes the exact numerical result, the dashed line is the prediction of Eq. (10), while the dash-dotted line is the prediction of Eq. (11) including $`m=50`$ terms in the summation. One can observe that even though the asymptotic result gradually fails (as it should) for larger $`h`$, the full series prediction of Eq. (11) remains very close to the fully numerical result. The difference can be well accounted for by the approximations involved in the continuum ansatz.

Additional insight on the appearance of the instability is given in Fig. 4, where we examine the mode stability for different values of $`\alpha `$ (at $`h=1`$). The bottom panel of the figure shows the trajectory of the main eigenmodes of the point spectrum of the linearization around the solution. The solid line shows the translational or pinning (antisymmetric) mode, while the dashed line shows the edge or breathing (spatially symmetric) mode . Notice that for $`h=1`$ (cf. Fig. 1 of Ref. ), the pinning mode has a larger frequency than the breathing mode. The band edge of the continuous spectrum which lies at $`\mathrm{\Lambda }=0.5`$  is shown by a dash-dotted line. As $`\alpha `$ is increased, it is clear that the breathing mode is “repelled” by the presence of the impurity towards the phonon band edge. On the contrary, the antisymmetric mode gradually approaches the origin and becomes unstable (dotted part of the relevant curve) for $`\alpha >0.22`$ (thereafter, for $`\alpha >\alpha _c`$ there appears an imaginary pair of eigenfrequencies). Perhaps it is even more interesting to examine the occurrence of this instability in the light of the top panels of the figure. The latter show the $`L^2`$ and $`L^{\mathrm{}}`$ norms of the solution as a function of $`\alpha `$. It is noteworthy that the occurrence of the instability coincides with the change of concavity of the $`L^2`$ norm, while it also coincides with a maximum of the $`L^{\mathrm{}}`$ norm. The above results suggest two additional numerical criteria that can be used to identify the appearance of the instability, namely

$`{\displaystyle \frac{^2u_2}{\alpha ^2}}=0,`$ (12)
and

$`{\displaystyle \frac{u_{\mathrm{}}}{\alpha }}=0;{\displaystyle \frac{^2u_{\mathrm{}}}{\alpha ^2}}<0.`$ (13)
## IV Dynamical effects and asymmetric modes

To examine the dynamical development of the instability, once the solution becomes unstable for the supercritical values of $`\alpha >\alpha _c`$, we have performed direct numerical simulations of Eq. (1) with an initial condition consisting of an exact, unstable localized mode for $`\alpha =0.3`$, perturbed by small noise (of uniform distribution, and amplitude ranging from 0.0001 to 0.1). The scenario that is described was found to be generic for such perturbations of the unstable solution. It was thus found that as is shown in the top left panel of Fig. 5 for $`t=200`$, the unstable mode sustains a symmetry breaking, which cleaves it into another stable localized mode centered essentially at the neighboring (previous) lattice site ( see also the bottom right panel), and to a small-amplitude mobile pulse that propagates along the chain and carries an energy excess. It is worth noting that the original (unstable) eigenmode had 3 “main” sites oscillating at an amplitude of $`0.9`$, while the resulting breather is very strongly localized at a single site (of amplitude $`1.7`$).

The trapping of the unstable mode by the neighboring site is consistent with the physical picture described by the effective Hamiltonian (6), when the mode overcomes the neighboring Peierls-Nabarro barrier and becomes localized by the next minimum. However, for smaller values of the lattice spacing $`h`$, when the Peierls-Nabarro barrier becomes negligible, the unstable localized mode can start its motion through the lattice under the action of the initial perturbation and as a result of the development of the symmetry-breaking instability. An example of the latter type is shown in Fig. 6 for the defect strength $`\alpha =1`$ and $`h=0.5`$.

One of the most interesting findings of our study of the time evolution as a result of the mode instability is the effective “relaxation” of the unstable mode towards the mode centered at the neighboring (to the defect) site. The latter result prompted us to explore the modes in the vicinity of the defect. Due to symmetry, our study was restricted to modes on one side of the impurity site. The bifurcation diagram of the relevant branches of solutions is shown in Fig. 7 obtained for $`h=0.5`$. We have observed the relevant phenomenology to be generic, but report it for this value of the lattice spacing that necessitates narrower ranges of parameter sweeps (i.e., the same observations can be obtained for larger $`h`$, but due to the weaker coupling, they occur for considerably larger values of the defect strength $`\alpha `$).

The natural starting point for our analysis is the continuation from the case of an homogeneous lattice ($`\alpha =0`$), where the stationary modes are well-known to exist on a lattice site and between two consecutive lattice sites. Hence, in the homogeneous case we consider the branches centered at $`n_0=100`$ (branch a), $`n_0=99.5`$ (branch b), $`n_0=99`$ (c), $`n_0=98.5`$ (d), $`n_0=97.5`$ (e) and $`n_0=97`$ (f). We observe that at the instability point $`a_c0.00252`$, the stable branch a and the unstable branch b merge. The resulting branch is always unstable thereafter. For any other branch apart from the central one, we have found that the relevant stable (node) and unstable (saddle) solutions exist for a parameter interval that depends on the distance of the central site of the branch from the defect site; naturally, the branches whose central site is more remote from the defect are less affected by it (and exist for larger intervals of defect strengths). These branches eventually disappear in saddle-node bifurcations, as is observed for branches c and d and e and f in Fig. 7. A set of spatial profiles corresponding to the different branches are shown in Fig. 8 for the different branches of Fig. 7 and for $`\alpha =0.0005`$ (corresponding to the vertical dotted line in Fig. 7).

Disappearance of the relevant branches of static solutions is consonant with the mobility of the unstable modes for smaller values of $`h`$ observed in Fig. 6 above. For smaller values of $`h`$, the neighboring branches disappear for smaller defect strengths, hence allowing for no static solutions in the vicinity (in configuration space) of the unstable solution, thereby resulting in the mobility of the latter (i.e., the solution cannot be trapped in neighboring potential wells as these have disappeared). However, in the opposite case of large $`h`$, the lattice is strongly discrete and an adiabatic growth of the defect strength $`\alpha `$ leads to the mode motion from one stable state to its neighboring state (via an intermediate unstable mode) through a cascade of the bifurcation points.

## V Conclusions

We have studied the existence and stability of nonlinear localized modes in the framework of the DNLS model. We have used a variational formalism to obtain the analytical result for the occurrence of an instability due to the interplay between attractive nonlinearity and repulsive impurity. We have examined the variational prediction by means of numerical bifurcation theory, linear stability analysis, as well as direct integration of the discrete lattice model. We have found the variational prediction (and its refinements) to be in reasonable agreement with the numerical results, attributing the discrepancies to the variation of the original profile with respect to the pre-selected ansatz. We have observed that the instability develops much before the solution becomes two-humped, and we have developed numerically motivated criteria for the identification of the instability. The study of the dynamical evolution of the unstable modes revealed the potential for mobile localized modes for smaller values of the lattice spacing, while it results in relaxation to the neighboring potential wells and the corresponding stable, static solutions. Finally, the bifurcation diagram has been explored as a function of the defect strength, using continuation from the homogeneous lattice case. We have found that for the branch corresponding to the defect site, the stable and nearest unstable solution branches merge at the instability point (producing a thereafter always unstable branch), while the neighboring sites corresponding stable and unstable branches disappear through collision via saddle-node bifurcations. Such a bifurcation scenario is in a perfect agreement with the concept of the Peierls potential that affects the mode motion in discrete systems.

PGK would like to thank J. Cuevas for stimulating discussions and for bringing to his attention the results of  prior to publication. He would also like to acknowledge partial support through a Faculty Research Grant of the University of Massachusetts and the Clay Foundation for support through a Special Project Prize Fellowship and the National Science Foundation for support through DMS-0204585. YSK acknowledges a support from the US Air Force-Far East Office and the Australian Research Council. ASK would like to acknowledge partial support from Royal Swedish Academy of Science.


text:
# Strategy updating rules and strategy distributions in dynamical multiagent systems

\[

## Abstract

In the evolutionary version of the minority game, agents update their strategies (gene-value $`p`$) in order to improve their performance. Motivated by recent intriguing results obtained for prize-to-fine ratios which are smaller than unity, we explore the system’s dynamics with a strategy updating rule of the form $`pp\pm \delta p`$ ($`0p1`$). We find that the strategy distribution depends strongly on the values of the prize-to-fine ratio $`R`$, the length scale $`\delta p`$, and the type of boundary condition used. We show that these parameters determine the amplitude and frequency of the the temporal oscillations observed in the gene space. These regular oscillations are shown to be the main factor which determines the strategy distribution of the population. In addition, we find that agents characterized by $`p=\frac{1}{2}`$ (a coin-tossing strategy) have the best chances of survival at asymptotically long times, regardless of the value of $`\delta p`$ and the boundary conditions used.

\]

The Minority Game (MG) is a successful model describing a population of competing and evolving individuals. This complex system has been explored extensively in the last few years, see e.g.,  and references therein. The present work is mainly motivated by the recent results of .

In this toy model, a population of $`N`$ agents with limited information and capabilities repeatedly compete for a limited global resource, or to be in the minority. The desire to be in a minority group is found in many real life situations, such as: financial markets, traffic jams, or among a group of predators (who wish to hunt in areas with fewer competitors).

At each round of the game, every individual has to choose whether to be in room ‘0’ (e.g., choosing to sell an asset) or in room ‘1’ (e.g., choosing to buy an asset). At the end of each turn, agents belonging to the smaller group (the minority) are the winners, each of them gains $`R`$ points (the “prize”), whereas the others lose a point (the “fine”). The agents share a common look-up table, containing the outcomes of recent occurrences. This allows the determination of a “predicted trend” in the system, which is followed by each agent with probability $`p`$, known as the agent’s “gene” value.

In the evolutionary formulation of the model (EMG) agents are allowed to evolve (“mutate”) their strategies based on past experience. If an agent score falls below some value $`d`$, he mutates – its gene value is modified. In this sense, each agent tries to learn from his past mistakes, and to adjust his strategy in order to perform better.

A remarkable conclusion deduced from the EMG  is that a population of competing agents tends to self-segregate into opposing groups characterized by extreme behavior. It was realized that in order to flourish in such situations, an agent should behave in an extreme way ($`p=0`$ or $`p=1`$) . On the other hand, in many real life situations the prize-to-fine ratio may take a variety of different values . A different kind of strategy may be more favorite in such situations. In recent studies it was found  that an intriguing phase transition exist in the model: “confusion” and “indecisiveness” take over when the prize-to-fine ratio falls below some critical value, in which case agents characterized by a coin-tossing strategy ($`p=\frac{1}{2}`$) perform better than extreme ones. In such circumstances agents tend to cluster around $`p=\frac{1}{2}`$ (see Fig. 1 of Ref. ) rather than self-segregate into two opposing groups.

In  we have considered a uniform strategy updating rule in which the new strategy (of a mutating agent) is chosen uniformly within the range $`0p1`$. Burgus, Ceva and Perazzo  have recently considered the same model problem, with an updating rule of the form $`pp\pm \delta p`$, where $`\delta p<\frac{1}{2}`$, and found that the population tends to form an M-shaped strategy distribution in the $`R<1`$ case. In the present work we further explore this system, and provide some new insights that extend and link the results of  to those of .

First, we would like to stress the importance of the chosen boundary conditions in the case of an updating rule of the form $`pp\pm \delta p`$ . Figure 1 displays the long-time averaged gene distribution $`P(p)`$ of the agents for two different types of boundary conditions: periodic and reflective. One finds that for periodic boundary conditions the population tends to cluster at intermediate gene values. The curve between the two peaks, located at $`p=\delta p`$ and $`p=1\delta p`$, is almost flat, while agents with extreme gene values ($`p0`$ and $`p1`$) perform much worse (we shall shortly demonstrate that the gene distribution may also have an inverse-U shape, depending on the precise values of $`R`$ and $`\delta p`$). On the other hand, the gene distribution is almost flat for reflective boundary conditions.

The underlying mechanism which is responsible for this important difference are the temporal oscillations observed in the winning probabilities of the agents . Figure 2 displays the time dependence of the winning probability of a $`p=0`$ agent (the winning probability of a central agent, with $`p=\frac{1}{2}`$, is practically constant in time). We consider three distinct cases, characterized by: (i) $`\delta p=0.1`$ with periodic boundary conditions, (ii) $`\delta p=0.1`$ with reflective boundary conditions, and (iii) uniform updating rule. One finds smaller oscillation amplitudes and longer periods for reflective boundary conditions, as compared to the case of periodic boundary conditions. This implies that, for reflective boundary conditions the performance of extreme agents ($`p=0`$ and $`p=1`$) becomes quite similar to the performance of central agents (characterized by $`p=\frac{1}{2}`$), implying a flatter gene distribution for these boundary conditions. On the other hand, for periodic boundary conditions one finds that the temporal oscillations are much more similar to the uniform case studied in  (as compared to the case of reflective boundary conditions). Indeed, the ratio $`P(\frac{1}{2}):P(0)`$ for periodic boundary conditions is very similar to the corresponding ratio in the uniform case (compare Fig. 1 with Fig. 1 of ).

Figure 3 shows the strategy distribution of the population for different prize-to-fine ratios, and with $`\delta p1`$. The results demonstrate the existence of a stable phase characterized by an inverse-U shaped gene distribution. However, unlike the uniform case , the critical value of $`R`$ which separates the inverse-U distribution from the M-shaped one does not equal $`1`$ (in the $`N\mathrm{}`$ limit).

In Fig. 4 we display $`P(p)`$ for various different $`\delta p`$ values with periodic boundary conditions. We find that the peaks of the strategy distribution (for prize-to-fine ratios which are large enough to allow an M-shaped gene distribution) occurs at $`p=\delta p`$ and its symmetric counterpart $`1\delta p`$. Regardless of the value of $`\delta p`$, the agents do not self segregate – the extreme strategies ($`p=0`$ and $`p=1`$) perform worst. The strategy distribution moves smoothly into an inverse-U shape in the limit of $`\delta p=\frac{1}{2}`$ . Figure 5 displays the same results for reflective boundary conditions, where $`\delta p=1`$ is equivalent to the uniform updating rule .

Figure 6 displays the average lifespan $`<`$$`L(p)`$$`>`$ of the agents. In order to get a better picture of the lifespan distribution, we also plot $`<`$$`L(p)`$$`>`$+$`\sigma _L(p)`$ as a function of the gene value $`p`$. Here $`\sigma _L(p)`$ is the root mean square separation of the lifespans. In this case, one finds an inverse-U shaped distribution (with the peak occurs at $`p=\frac{1}{2}`$). This implies that agents characterized by $`p=\frac{1}{2}`$ (a coin-tossing strategy) have the best chances of survival at asymptotically long times, as predicted analytically in . This important feature is explained by the global currents in the gene-space, which reduce the value of $`\sigma _L(p=0)`$, and have a negligible effect on $`\sigma _L(p=\frac{1}{2})`$ . We emphasize that these results hold true for both periodic and reflective boundary conditions.

The efficiency of the system is defined as the number of agents in the minority room, divided by the maximal possible size of the minority group, $`(N1)/2`$. Figure 7 displays the efficiency as a function of the length scale $`\delta p`$. The system’s efficiency is a monotonic decreasing function of $`\delta p`$. This is caused by the fact that larger $`\delta p`$ values imply larger temporal oscillations in the occupation numbers of the rooms, thus decreasing the number of agents in the winning group (and increasing the number of agents in the losing room).

Finally, we would like to address the last point raised in . It is claimed that the fluctuations in the average gene value $`<`$$`p`$$`>`$ have been considered in . However, the oscillatory behavior of $`<`$$`p`$$`>`$, which is an highly important feature of the system’s dynamics was not observed in . Rather, Burgos et. al.  find a non-oscillatory value for $`<`$$`p`$$`>`$-$`\frac{1}{2}`$, see Eq. $`(15)`$ of . We have shown, on the other hand, that the quantity $`<`$$`p`$$`>`$-$`\frac{1}{2}`$ displays temporal oscillations, with well defined frequency and amplitude . It is important to distinguish between regular temporal oscillations of the physical quantities (such as $`<`$$`p`$$`>`$) discussed in , as opposed to thermal fluctuations discussed in . Thermal fluctuations of a thermodynamic system are essentially random in nature, whereas we have found regular oscillations, which are characterized by well defined frequency and amplitude. The oscillatory nature of $`<`$$`p`$$`>`$  has been proven to be an essential feature which is responsible for the dynamical phase transition (from self-segregation to clustering) observed in the EMG . We would like to emphasize that these oscillations exist also for complex systems with a strategy updating rule of the form $`pp\pm \delta p`$, regardless of the value of $`\delta p`$ and the type of boundary conditions used (see. Fig. 2).

ACKNOWLEDGMENTS

SH thanks a support by the Dr. Robert G. Picard fund in physics. This research was supported by grant 159/99-3 from the Israel Science Foundation. EN thanks a support by the Horwitz foundation.


text:
# Cloning the entanglement of a pair of quantum bits

## Abstract

It is shown that any quantum operation that perfectly clones the entanglement of all maximally-entangled qubit pairs cannot preserve separability. This “entanglement no-cloning” principle naturally suggests that some approximate cloning of entanglement is nevertheless allowed by quantum mechanics. We investigate a separability-preserving optimal cloning machine that duplicates all maximally-entangled states of two qubits, resulting in 0.285 bits of entanglement per clone, while a local cloning machine only yields 0.060 bits of entanglement per clone.

Quantum entanglement is known to be a resource that is central to many quantum information processes such as quantum teleportation, quantum cryptography, or quantum computing nielsen . In view of this, much work has been devoted to defining measures of entanglement or to investigating the best information-theoretical use of entanglement. Despite the fact that entanglement is a very fragile resource, extremely sensitive to decoherence, several techniques have been developed in order to overcome decoherence, namely quantum error correction or entanglement purification (see BDSW ). Out of these many studies of entanglement, none has so far addressed the issue of whether (and how well) entanglement can be cloned.

There has been a vast literature about the question of cloning quantum states. First of all, the no-cloning theorem has been stated WZ , which precludes the perfect copying of an arbitrary quantum state. Then, imperfect quantum cloning machines (QCM) have been introduced, which duplicate an arbitrary qubit state with the highest possible fidelity BH . Since then and more recently, a large variety of QCMs have been devised, with the purpose of cloning equally well a given set of states in a space of arbitrary dimension (see Werner ).

In this paper, we raise the question of whether quantum entanglement itself can be cloned or not. In order to simplify our analysis, we restrict ourselves to qubit pairs (dimension $`2\times 2`$). We show that the requirement of perfectly cloning the entanglement carried by a qubit pair in an arbitrary maximally-entangled (ME) state is incompatible with the requirement that separable qubit pairs remain unentangled via cloning. Of course, if we restrict ourselves to four orthogonal ME states such as the Bell states

$$|\mathrm{\Phi }^\pm =(|00\pm |11)/\sqrt{2},|\mathrm{\Psi }^\pm =(|01\pm |10)/\sqrt{2},$$
(1)
then we can very well make a Bell measurement of the original pair and subsequently prepare an arbitrary number of clones in the measured ME state. However, this procedure does not work properly on the linear combinations of Bell states that remain ME states. Even if the clones are not required to be close to the original ME 2-qubit state but solely to be entangled, it remains impossible to fully preserve entanglement.

Here we show that it is nevertheless possible to clone part of the original entanglement, much in the same way quantum states can be cloned imperfectly. We construct a QCM that is universal over the set of ME 2-qubit states and argue that it effects an optimal cloning of entanglement. For this purpose, we exploit the property that the set of ME 2-qubit states is isomorphic to the set of real states in 4 dimensions BDSW , from which we construct an optimal symmetric cloner that is covariant under local unitaries. This cloner maximizes the amount of entanglement left in the clones (about 0.285 ebit per clone), while product states remain unentangled. Finally, we consider asymmetric QCMs and investigate how entanglement is distributed among the clones by these transformations.

Entanglement no-cloning principle: Entanglement cannot be cloned perfectly, that is, if a quantum operation can be found that perfectly duplicates the entanglement of all ME states, then it necessarily does not preserve separability (some separable states become entangled after cloning).

Proof: We restrict ourselves to 2 qubits and consider two orthogonal ME states, e.g., $`|\mathrm{\Phi }^\pm `$. Assume that the entanglement of these states is perfectly cloned, i.e., the output states of the clones remain ME even though they may differ from the input state. The most general cloning transformation $`U`$ preserving the entanglement of these two states can be written as

$`|\mathrm{\Phi }^\pm |A\stackrel{U}{}|e_a^\pm |e_b^\pm |A^\pm ,`$ (2)
where $`|A`$ denotes the initial state of the ancilla and the blank copy, while $`|A^\pm `$ is the ancilla state after cloning. Thus, the states of the two clones $`|e_a^\pm `$ and $`|e_b^\pm `$ are some ME states. Now, the linear combination $`|\stackrel{~}{\mathrm{\Phi }}=(|\mathrm{\Phi }^++i|\mathrm{\Phi }^{})/\sqrt{2}=(e^{i\pi /4}|00+e^{i\pi /4}|11)/\sqrt{2}`$ is still a ME state. By linearity, the above transformation yields the following output state

$`|\stackrel{~}{\mathrm{\Phi }}|A\stackrel{U}{}(|e_a^+|e_b^+|A^++i|e_a^{}|e_b^{}|A^{})/\sqrt{2}.`$ (3)
In order to preserve the full entanglement within each clone, a necessary condition is that either $`|e_a^+=|e_a^{}`$ or $`|e_b^+=|e_b^{}`$. However, in each of these two cases, at least one of the clones is left in a ME state that is independent of the input state (within the space spanned by $`|\mathrm{\Phi }^\pm `$) regardless of it being separable or not. For example, in the first case, if the input is the separable state obtained as the linear combination $`|s=(|\mathrm{\Phi }^++|\mathrm{\Phi }^{})/\sqrt{2}=|00`$, then the transformation gives

$`|s|A\stackrel{U}{}|e_a(|e_b^+|A^++|e_b^{}|A^{})/\sqrt{2}.`$ (4)
Clearly, the separability is not preserved here since the first clone is maximally entangled. We therefore conclude that no perfect cloning of entanglement is possible. $`\mathrm{}`$

As a consequence, only imperfect QCMs that approximately reproduce the entanglement while preserving separability can be implemented. In the rest of this paper, we will be interested in separability-preserving QCMs that yield clones with the highest achievable entanglement for all ME input states. As shown later on, finding these QCMs is strongly related to finding transformations that clone optimally and equally well the set of 2-qubit ME states. Consider an arbitrary 2-qubit pure state

$$|\mathrm{\Phi }=\underset{i=0}{\overset{3}{}}n_i|e_i$$
(5)
written in the orthonormal basis made of the Bell states with particular phases (sometimes referred to as the magic basis BDSW ):

$`|e_0=|\mathrm{\Phi }^+,|e_1=i|\mathrm{\Phi }^{},|e_2=i|\mathrm{\Psi }^+,|e_3=|\mathrm{\Psi }^{},`$ (6)
where the amplitudes $`n_i`$ are normalized as $`_{i=0}^3|n_i|^2=1`$. In this basis, the entanglement of formation $`E`$ of the state $`|\mathrm{\Phi }`$ can be expressed in a very simple way as

$$E(𝒞(\mathrm{\Phi }))=H\left(\frac{1}{2}+\frac{1}{2}\sqrt{1𝒞(\mathrm{\Phi })^2}\right),$$
(7)
where $`H`$ is the binary entropy function and

$$𝒞(\mathrm{\Phi })=|\underset{i}{}n_i^2|$$
(8)
is called the concurrence BDSW ; HW . Clearly, any real linear combination (up to an irrelevant global phase) of the magic basis elements is a ME state since $`𝒞`$ (and therefore $`E`$) is then equal to one. Furthermore, every ME state can be expressed as a real linear combination of the magic basis elements. For this reason, the problem of cloning the set of ME 2-qubit states boils down to constructing a transformation that optimally clones all real 4-dimensional states in the magic basis.

This particular transformation can be found by following a method inspired from Cerf . We assume that the qubit pair to be cloned is itself maximally entangled with a reference system (another qubit pair). We then explicitely construct the most general joint state describing the reference $``$, the output clones $`a`$ and $`b`$, and an ancilla $`𝒜`$ after the cloning transformation

$$|𝒮_{,a,b,𝒜}=\underset{i,j,k,l=0}{\overset{3}{}}s_{ijkl}|l_{}|i_a|j_b|k_𝒜.$$
(9)
(The reference, the two clones, and the ancilla are all 4-dimensional systems here.) This state serves to completely define the cloning transformation: the result of cloning the ME state $`|\mathrm{\Phi }=_in_i|e_i`$ (with real $`n_i`$’s) is obtained by projecting $``$ onto the complex conjugate $`|\mathrm{\Phi }^{}`$ (which is equal to $`|\mathrm{\Phi }`$ here) Cerf . Thus, the most general cloning transformation is defined as

$$|\mathrm{\Phi }\underset{i,j,k,l=0}{\overset{3}{}}s_{ijkl}n_l|i_a|j_b|k_𝒜.$$
(10)
At this point, we impose the additional condition that the QCM is covariant under $`SU(2)\times SU(2)`$ in the computational basis (or, equivalently, under $`SO(4)`$ in the magic basis). This means that the QCM acts similarly in all bases connected by local unitaries to the computational basis. This restriction is natural since it guarantees that all states equivalent up to local unitaries (thereby having the same entanglement) result in equally entangled clones. A sufficient condition for covariance is patrick

$$|𝒮_{,a,b,𝒜}=R^4|𝒮_{,a,b,𝒜},$$
(11)
where $`R`$ is any real rotation matrix in $`SO(4)`$. This requirement implies that $`s_{ijkl}`$ is an invariant tensor of rank four, that is, $`s_{ijkl}=R_{ii^{}}R_{jj^{}}R_{kk^{}}R_{ll^{}}s_{i^{}j^{}k^{}l^{}}`$. A main simplification here results from the fact that the generic form of such a tensor is

$$s_{ijkl}=A\delta _{il}\delta _{jk}+B\delta _{jl}\delta _{ik}+C\delta _{kl}\delta _{ij},$$
(12)
with the normalization condition on (10) imposing that

$$4(|A|^2+|B|^2+|C|^2)+2\mathrm{Re}(AB^{}+AC^{}+BC^{})=1.$$
(13)
For a symmetric cloner, the permutation symmetry between the two clones imposes furthermore that $`A=B`$, so that we are left with a transformation depending on two parameters, $`A`$ and $`C`$. If we use the cloning fidelity as a figure of merit, Eqs. (10) and (12) result in

$`F=\mathrm{\Phi }|\rho _{a,b}|\mathrm{\Phi }=7|A|^2+|C|^2+4\mathrm{Re}(AC^{}),`$ (14)
where $`\rho _{a(b)}`$ denotes the reduced density matrix of clone $`a(b)`$. This expression can be maximized under the normalization constraint Eq. (13), giving

$`A={\displaystyle \frac{1}{3}}\left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{\sqrt{13}}}\right)^{1/2},C={\displaystyle \frac{A}{2}}(\sqrt{13}3),`$ (15)
with the corresponding fidelity

$$F=\frac{5+\sqrt{13}}{12}0.7171.$$
(16)
Interestingly, this maximization procedure yields a fidelity that saturates the upper bound derived from the no-signaling condition in patrick . We therefore conclude that the optimal covariant cloner of ME 2-qubit states is characterized by Eq. (15) and yields the fidelity (16). This is slightly higher than the fidelity of the universal 4-dimensional cloner, namely $`F=7/10`$ Werner ; Cerf , as expected since the ME states form a subset of the 2-qubit states.

We now generalize Eq. (14) to asymmetric cloning transformations ($`AB`$). The cloning fidelities become

$`F_a=4|A|^2+|B|^2+|C|^2+2\mathrm{Re}(AB^{}+AC^{}+BC^{}),`$
$`F_b=4|B|^2+|A|^2+|C|^2+2\mathrm{Re}(AB^{}+AC^{}+BC^{}).`$
Assuming $`A`$, $`B`$, and $`C`$ to be real, we eliminate $`A`$ and $`C`$ from Eqs. (13) and (Cloning the entanglement of a pair of quantum bits), to get

$`F_a(B,F_b)=3B^2+{\displaystyle \frac{F_b+1}{2}}+{\displaystyle \frac{\sqrt{3B^2+F_b}B}{2}}`$
$`\times \sqrt{18B^2+18B\sqrt{3B^2+F_b}15F_b+6}.`$ (18)
We then maximize $`F_a`$ over $`B`$ as a function of $`F_b`$. Figure 1 displays the resulting $`F_a`$ as a function of $`F_b`$. The mid-point value of this curve lies at $`F_a=F_b=F`$, as expected. We also confirm that the fidelity of a clone is equal to one when the other one is in a completely mixed state, i.e., when its fidelity equals 1/4.

We have numerically confirmed the optimality of this class of asymmetric cloners with the use of a technique based on semidefinite programming Fiurasek01 ; Fiurasek03 . The cloning transformation is a linear trace-preserving completely positive (CP) map that can be represented by a positive semidefinite operator $`S`$ on the tensor-product space of the input and output states. The cloning fidelities can be expressed as $`F_{a(b)}=\mathrm{Tr}[SR_{a(b)}]`$, with appropriately defined $`R_{a(b)}0`$ Fiurasek03 . The optimal asymmetric cloner can be obtained by maximizing $`F=pF_a+(1p)F_b`$, where $`p[0,1]`$ is the asymmetry parameter. The resulting fidelities coincide with those obtained from Eq. (Cloning the entanglement of a pair of quantum bits) up to the machine precision.

Let us now investigate the entanglement properties of this cloning transformation and show that it is also optimal with respect to our original goal, namely cloning the amount of entanglement. Let us start by checking that it preserves separability. For the ansatz (9), we have $`S=\mathrm{Tr}_𝒜(|𝒮𝒮|)`$, and the CP map that describes the relationship between the input and the clone $`a`$ ($`b`$) can be characterized by $`S_{a(b)}=\mathrm{Tr}_{b(a)}[S]0`$. Since the positive partial transpose (PPT) criterion is a necessary and sufficient separability condition for a qubit pair, a sufficient condition for these two maps to preserve separability is that $`S_a`$ and $`S_b`$ represent PPT operations Rains01 . The map $`S_{a(b)}`$ is PPT if $`S_{a(b)}^{T_{1,1^{}}}0`$, where $`T_{1,1^{}}`$ denotes partial transposition with respect to the first qubit of the original and the clone $`a(b)`$. An explicit analytical calculation shows that if $`s_{ijkl}`$ is an invariant rank-4 tensor (12), then $`S_{a(b)}^{T_{1,1^{}}}=S_{a(b)}0`$, hence a covariant cloner necessarily preserves separability.

It is therefore natural to maximize the output entanglement for the other extreme case, namely when the original qubit pair is maximally entangled. The amount of entanglement left in the clones will be measured here by the entanglement of formation $`E`$, which can be evaluated by using the extended definition of the concurrence $`𝒞`$ for mixtures HW . The entanglement of formation of an arbitrary 2-qubit state $`\rho `$ is given by $`E(\rho )=E(𝒞(\rho ))`$, where $`𝒞(\rho )=\mathrm{max}(0,\lambda _1\lambda _2\lambda _3\lambda _4)`$ and the $`\lambda _i`$’s are the eigenvalues, in decreasing order, of the Hermitian matrix $`\stackrel{~}{\rho }\sqrt{\sqrt{\rho }\rho ^{}\sqrt{\rho }}`$. Here $`\rho ^{}`$ denotes the complex conjugate of $`\rho `$ when expressed in the magic basis HW . For a generic covariant cloner, Eq. (12), the reduced density matrices of the clones of an input ME state $`|e_i`$ can always be written as

$`\rho _{a,b}=F_{a,b}|e_ie_i|+{\displaystyle \frac{1F_{a,b}}{3}}{\displaystyle \underset{ji}{}}|e_je_j|,`$ (19)
so the clones are left in a mixture of (generalized) Bell states, or more precisely in a Werner state. This is consistent with the fact that $`\rho _{a,b}`$ are real density matrices in the magic basis because the input state is real HW . Hence, $`\stackrel{~}{\rho }_{a,b}=\rho _{a,b}`$ so the concurrence of the clones simply reduces to $`𝒞_{a,b}=\mathrm{max}(0,2F_{a,b}1)`$. Therefore, for a covariant cloner, maximizing $`E`$ reduces to maximizing $`F`$. Consequently, the cloner characterized by Eq. (15) is the entanglement cloner we were looking for, provided that covariance is taken for granted.

The corresponding entanglement of formation of the clones, $`E_a`$ and $`E_b`$, is shown in Fig. 1 for different values of $`F_b`$. As expected, the entanglement of formation of clone $`b`$ vanishes for $`F_b1/2`$. Conversely, the entanglement of formation of clone $`a`$ vanishes when $`F_b0.8984`$ (that is, when $`F_a1/2`$). Note that the entanglement of formation of a clone is equal to one only when its fidelity is one, thus confirming that a fully asymmetric cloner (a trivial cloner which outputs the original and a random clone) is the only solution if we want to fully conserve the original entanglement of 1 ebit. Finally, note that the sum of the entanglement of formation of the two clones (also shown in Fig. 1) never exceeds one, meaning that the entanglement cloner does not create more entanglement than that contained in the original ME state.

For the symmetric cloner, the entanglement of formation of both clones is equal to $`E_a=E_b=0.2847`$ bits. The optimality of this result has been verified using numerical optimization where the structure of the cloning transformation was based on the parametrization proposed in Cerf  and the maximized quantity was the concurrence instead of the fidelity. Up to irrelevant local unitaries (which decrease $`F`$ while keeping $`𝒞`$ and $`E`$ constant), we recovered the same cloning transformation. This strongly suggests that restricting ourselves to covariant QCMs is justified, so the cloner of ME states also optimally clones the amount of entanglement.

Finally, the properties of our entanglement cloner can be analyzed in the intermediate case of non-maximally entangled input states. In Figure 2, we plot the entanglement of formation of the clones as a function of that of the original state $`\alpha |00+\sqrt{1\alpha ^2}|11`$ ($`0\alpha 1`$) for the symmetric cloner. No entanglement is cloned below a critical value $`E_{\mathrm{in}}=0.161`$ bits. Then, the output entanglement increases monotonically to reach its maximum $`E_{\mathrm{out}}=0.2847`$ bits at $`E_{\mathrm{in}}=1`$ bit.

In conclusion, we have shown that the quantum entanglement of an unknown ME qubit pair cannot be perfectly cloned if, at the same time, product states are required to be cloned into unentangled qubit pairs. In other words, a separability-preserving QCM cannot perfectly duplicate the entanglement of the set of ME states. Only imperfect QCMs do exist. As a first step, we have constructed an optimal symmetric entanglement cloner which is universal over the set of ME states, and whose fidelity saturates the no-signaling upper bound patrick . This cloner yields imperfect clones with 0.285 ebits if the original qubit pair contains 1 ebit, while unentangled pairs are cloned into separable states. The distribution of entanglement among the clones has also been investigated using an asymmetric cloner. Similarly to the situation that prevails when cloning quantum states, this no-go theorem for entanglement cloning might be exploited in order to imagine new quantum key distribution schemes. For example, one could imagine a protocol where the eavesdropper is only able to apply a local cloning on a ME state instead of the above global cloning. One can check that applying the optimal universal qubit cloner BH  locally on each qubit of a ME state reduces the fidelity of the clones to $`7/12`$. Since the reduced density matrix of the clones is again a Werner state, we obtain $`𝒞=1/6`$ leading to only about 0.060 ebits per clone. Thus, the fact that independent local operations on each qubit are less efficient than joint operations could be used to give an advantage to the authorized parties.

We thank T. Durt and S. Iblisdir for helpful discussions. We acknowledge funding from the Communauté Française de Belgique under grant ARC 00/05-251, from the IUAP programme of the Belgian government under grant V-18, and from the EU under project CHIC (IST-2001-32150). JF also acknowledges funding from the grant LN00A015 of the Czech Ministry of Education.