text:
# Production of the Large Scale Superluminal Ejections of the Microquasar GRS 1915+105 by Violent Magnetic Reconnection

## 1 Introduction

Galactic and extragalactic accretion-powered sources often exhibit quasi-periodic variability and ejection phenomena that may offer important clues to the general physical processes in the inner regions of all classes of accreting sources, from young stars to X-ray binaries and AGNs.

The X-ray source GRS 1915+105 was the first galactic object to show evidence of superluminal radio ejection (Mirabel & Rodriguez 1994, 1999) which is normally interpreted as due to relativistic jets. At a distance of $`12.5`$ kpc and probably with a 10 solar mass black hole in the center of a binary system, this microquasar offers an excellent laboratory for the investigation of black hole accretion and associated jet phenomena revealing, in time scales of minutes to months, a profusion of variable emission structures in a large range of wavelengths (Dhawan, Mirabel, & Rodríguez 2000).

From a detailed analysis of the observed RXTE X-ray variability of this source, Belloni et al. (2000) have distinguished 12 different modes of variability that they interpreted as transitions between three basic states: a hard spectral-index state (denoted by C) with low X-ray luminosity, and two states (denoted by A and B) with softer spectra that were modeled as thermalized emission from an accretion disk, with state B being more luminous and hotter than state A. They also found that after a long steady hard class, the most frequent class of variability was one characterized by a 30-minute cycle.

Variability also has been observed in the radio and IR wavebands. Mirabel et al. (1998) reported simultaneous X-ray, IR and radio observations of these 30-minute cycles during which the X-ray emission is found to oscillate between states A/B and C. The IR/radio emission of the plasma (which is consistent with synchrotron emission from an expanding nuclear jet; see below) seems to be associated with the times of recovery from the X-ray dips (see Figs. 2 and 3 from Mirabel et al. 1998), i.e., with the times the source is changing from a period in state C to a period in which it rapidly oscillates between states A/B and C.

Large scale superluminal radio ejections have been also detected with MERLIN emerging from the source when this was at the end of a 3-week extended phase of low/hard X-ray state (also designated the $`plateau`$ state; Fender et al. 1999; see also Fender et al. 2002). Similarly, in a compilation of (VLBA) radio and (RXTE and BATSE) X-ray data, taken during several weeks, Dhawan, Mirabel & Rodriguez (2000; hereafter DMR) have distinguished two main states of the system, a $`plateau`$ and a $`flare`$ state. The plateau state is characterized by a flat radio spectrum ($`S_\nu \nu ^0`$) coming from a compact region of size of a few AU, and flux density 1-100 mJy. During this phase, the associated RXTE (2-12 keV) soft X-ray emission is weak, while the BATSE (20-100 keV) hard X-rays are strong. On the other hand, during the flare phase, optically thin ejecta are superluminally expelled up to thousands of AU with fluxes up to 1 Jy at $`\lambda 13`$ cm ($`S_\nu \nu ^{0.6}`$) with rise time $`1`$ day and which fade after several days. The soft X-rays also flare during this phase and exhibit high variability, while the hard X-rays fade for few days before recovering (see Figs. 1 and 4 of DMR). The radio imaging over a range of wavelengths (13, 3.6, 2.0, and 0.7 cm) resolves the nucleus as a compact jet of length $`10\lambda _{cm}`$ AU. This nuclear jet varies in time scales of $`30`$ minutes during minor X-ray/radio bursts and reestablishes within $`18`$ hours after a major outburst. So far, GRS 1915 + 105 has been the only system among the galactic accreting black holes of stellar mass where both AU-scale steady jets and large scale superluminal ejections have been unambiguously observed.

It is generally believed that the X-ray emission of galactic black holes and active galactic nuclei arises from hot gas accretion flow or from an accretion-disk corona. In the first class of models, advection-dominated accretion flow (ADAF) of disk material into the black hole (e.g., Narayan, Mahadevan, & Quartaet 1997), as well as, simultaneous advection-dominated inflow-outflow (ADIOS) solutions (Blandford & Begelman 1999) have been proposed. In the second class, the presence of magnetic fields adds complexity to the solutions but plays an essential role in the disk-coronal coupling, in the energy and angular momentum transport, and in jet production and acceleration (e.g., Liu, Mineshige & Shibata 2002).

The AU-scale nuclear jet of GRS 1915+105 developed in the plateau phase has been interpreted as synchrotron emission from relativistic electrons (e.g., Falke & Biermann 1999), and its origin has been attributed to disk instabilities in an ADAF-type flow by Belloni et al. (1997a, b). According to this model, the 30-minute soft X-ray variability observed during the plateau state, which involves a transition between states C and A/B, is explained in terms of periodic evacuation and refilling of the inner disk region on time scales of seconds as a result of thermal viscosity. This process is believed to occur in the region where the radiation pressure dominates the gas pressure and the time required to refill the disk with X-ray emitting hot gas is set by the thermal viscous time scale. The observation of radio oscillations with similar period interval supports the idea that the synchrotron emission of the nuclear jet is fed by relativistic plasma injected from the inner disk during the dips in the soft X- ray flux (DMR).

In contrast, the large superluminal radio flares observed in the $``$ 500 AU scales cannot be explained by the same viscous disk instability model, because they eject an order of magnitude more mass than the AU-scale $`minijet`$ and the recovery time scale (of $``$ 18 hr) would require according to that model a much larger evacuation radius of disk material well beyond the radiation pressure-dominated region where the instability occurs.

Recently, different groups have proposed alternative scenarios to explain the variability of the GRS 1915+105 microquasar. Livio, Pringle & King (2003; see also King et al. 2004) have suggested that the inner region of the accretion disk would switch between two states. In one of them (state A/B), the accretion energy would be dissipated locally to produce the observed disk luminosity, and in the other (state C) the accretion energy would be converted into magnetic energy and emitted in the form of a relativistic jet (state C). They have attributed the transition between the two states to dynamo generation of a global poloidal magnetic field. Similarly, Tagger et al. (2004), also focusing mainly on the 30-minute cycles, have suggested that the alternations between the high-soft and the low-hard states and the production of a nuclear relativistic jet could be controlled by the successive accumulation and release of poloidal magnetic flux (e.g., by magnetic reconnection) in the innermost parts of the accretion disk. They have proposed a specific model to explain the quasi periodic variability in the disk emission based on the development of the accretion-ejection instability (Tagger & Pellat 1999). This is a global mode formed by spiral waves that travel back and forth in the inner regions of a disk threaded by a poloidal magnetic field with a ratio $`\beta `$ between the gas pressure and the magnetic pressure close to unity. The spiral waves become unstable by extracting energy and angular momentum from the inner region (thus increasing the accretion) and transferring it outward.

Also, in a contemporaneous work, de Gouveia Dal Pino & Lazarian (2004) have suggested the possibility that the superluminal large scale radio flares of GRS 1915+105 could be produced by relativistic plasma accelerated during violent magnetic reconnection events occurring in the magnetized corona just above the inner regions of the accretion disk. In the present work, we try to put together all the pieces of information above in a consistent picture to show that violent magnetic reconnection could indeed occur in the innermost parts of the disk and produce the large scale superluminal ejections when $`\beta `$ approaches one and the mass accretion energy becomes close to the critical value.

In ¶2 we draw the potential scenario for flare production and the superluminal ejecta. Then in ¶3, we try to quantify it by computing the disk and coronal parameters and evaluating the total amount of magnetic energy released during violent reconnection in the corona; and in ¶4, we evaluate both the electron and the synchrotron power-law spectra that is produced by the accelerated particles through first-order Fermi process in the reconnection site. We finally conclude (in ¶5) with a brief discussion of the implications of our results.

## 2 A Simple Possible Scenario

### 2.1 The BH magnetosphere

Magnetized accretion disks around rotating (Kerr) black holes (BHs) are frequently invoked to explain the high energy radiation and jet production and collimation in quasars and microquasars. The magnetic field lines originally frozen in the disk plasma will deposit along with accreting gas onto the BH horizon, therefore developing a magnetosphere around the later (see, e.g., Wang, Xiao, & Lei 2002). Near the horizon, no matter how chaotic may be the field threading the disk, the field through the hole will become very $`clean`$ and ordered. If the disk tries to deposit chaotic field on the hole, the field’s closed loops will destroy themselves on a time scale $`tR_H/c10^4`$ s $`M_{10}`$, leaving the field ordered. If on the other hand, magnetic pressure from field lines threading the disk temporarily push the hole’s field into a clumped configuration, it will spring back on this same time scale and make itself more uniform (MacDonald et al. 1986).

The magnetic coupling (MC) of a rotating black hole with its surrounding disk has been recently investigated in some detail (see e.g., Blandford 1999, Wang, Lei & Ma 2003, and references therein). The presence of some closed field lines connecting the BH with the inner disk edge (in a variant of the Blandford-Znajek process, 1977), can lead to energy and angular momentum transfer between the BH and the disk if they are rotating at different angular speeds. Though not a necessary condition, we will assume, for simplicity, in the present study that the BH and the inner disk edge are nearly co-rotating so that no significant angular momentum and energy transfer is occurring between them. The magnetic structure near the BH horizon will have a force-free configuration determined by the MC process and the strength of the field near the hole will be determined by the disk’s past history. The net flux today must be equal to the total flux that the disk has deposited in its entire past lifetime. These deposits will cause the net flux through the hole to fluctuate stochastically. The rms field will grow with time until it becomes so strong that it pushes its way back into the accretion disk and slows down accretion. The resulting maximum strength is likely to be such as to produce equilibrium between gas ram pressure and magnetic pressure in the innermost parts of the disk (Wang, Lei and Ma 2003), though Parker-Rayleigh-Taylor instabilities may reduce this strength (MacDonald et al. 1986).

### 2.2 The accretion disk and corona around the BH

In the standard model for an accretion disk around a black hole, gravitational energy is continually converted into kinetic, thermal, and magnetic energy by the joint action of viscosity, and magnetic field line stretching and reconnection. The viscosity can be mainly turbulent or magnetic, or both (Shakura-Sunyaev 1973).

In the process of angular momentum transport in the disks, shear is expected to convert an existing poloidal magnetic field into a toroidal component. Then, probably some of the toroidal flux will be used to generate new poloidal fields by raising magnetic loops in a dynamo process. Most of the energy will be either dissipated within the disk that will then emit thermalized radiation, or in the low density corona above and below the disk that will radiate thermally or by Compton scattering, or else the energy will be released in the form of a bulk flow, as a jet. In the last case, it is possible that essentially all the accretion energy of the inner disk is lost in the jet (see, Blandford & Payne 1982, Livio et al. 2003 and references therein).

Let us assume as in previous works (e.g., Merloni 2003, Livio et al. 2003, Tagger et al. 2004, de Gouveia Dal Pino & Lazarian 2004) that the inner parts of the accretion disk alternate between two states which are controlled by changes in the global magnetic field. Most of the magnetic field generated in a dynamo process in which the turbulence is driven by, e.g., the magneto-rotational instability (Balbus & Hawley 1991), is a tangled small-scale magnetic field within the disk. As remarked by Livio et al. 2003, as long as the magnetic field remains as such most of the dissipation will occur locally within the disk that will be radiative and very luminous. However, if a large scale poloidal magnetic field is generated then it is possible that an energetic outflow can be driven.

Numerical simulations of the dynamo process in accretion disks (e.g., Kudoh et al. 2002, Shibata 2003), are still unable to compute the generation of a global poloidal magnetic field, so that for the time being we must rely on physical estimations rather than on detailed modeling (Livio, Pringle and King 2003). Tagger et al. (2004), for example, have suggested that a gradual increase in the magnetic flux in the inner regions of the accretion disk could be achieved by the combined action of the turbulent dynamo associated, e.g., with the magneto-rotational instability (which would create a vertical flux of opposite polarities in the inner and outer disk regions), followed by advection of the magnetic flux with the gas to the inner disk regions. Though magnetic diffusivity could interrupt this magnetic flux accumulation, there is at least some observational evidence for it in the case of the magnetic field structure of our own galactic disk that shows a strong vertical field in the inner tens of pcs of the galactic center which is almost three orders of magnitude more intense than the horizontal magnetic fields measured elsewhere in the disk. Livio (1997; see also Livio et al. 2003, King et al. 2004), in addition, assuming a turbulent dynamo process with inverse cascade and reconnection of magnetic loops above and below the disk, have demonstrated that (provided that the distribution $`n(l)`$ of the lengths of the magnetic loops satisfies a relation $`n(l)l^{3/2}`$) large scale poloidal fields can be generated such that $`L_jL_{acc}`$, where $`L_j`$ is the jet power and $`L_{acc}`$ is the disk accretion luminosity.

Considering the aspects above, we assume here that during the $`plateau`$ state that precedes a radio flare, a large scale poloidal field is progressively built in the disk by a turbulent dynamo process. The action of buoyancy forces will also make the disk unstable against the Parker-Rayleigh-Taylor instability and horizontal magnetic field lines will raise from the disk forming large scale loops in the rarified hot corona. We further assume that once the dynamo process establishes a global poloidal field over a substantial region of the disk, it will be able to maintain that field for a period of time. The length of time the source spends in this state is hard to estimate without a complete physical model, however observations indicate that the system remains in the plateau state for prolonged periods of time.

The vertical field flux will give rise to a wind that will remove angular momentum from the disk, therefore significantly increasing the accretion rate (possibly at a rate greater than the rate due to the disk viscosity). Also, with the accumulation (trap) of vertical flux in the inner regions the ratio between the gas+radiation pressure to the magnetic field pressure ($`\beta `$) will soon decrease to one. Under these conditions, events of reconnection of magnetic field lines with opposite polarization will be inevitable, and in the innermost regions this process may become eventually very violent when enough magnetic energy is stored in the corona. We show below that this occurs when $`\beta 1`$ and the accretion rate approaches the critical value. We then argue that this could explain the large scale radio flares in GRS 1915+105 with the ejection of the observed relativistic blobs from the innermost regions.

### 2.3 Magnetic reconnection site in the inner disk-corona region

Considering the assumptions above, the resulting qualitative structure of the magnetosphere of the hole and the accretion disk must be like that shown in Figure 1 (which was adapted from Macdonald & Thorne 1982, and Macdonald et al. 1986). In the inner disk region (of radius $`R_X`$), there is a site which may be appropriate for violent magnetic reconnection. Surfaces of null poloidal field lines mediate the geometry of the open field lines anchored into the BH horizon with the opened lines of the disk wind and those connecting the disk with the BH horizon. Labeled as ”$`Y`$ $`netral`$ $`zone`$” in Figure 1, these magnetic null surfaces begin or end on $`Y`$ points. Across each null surface, the poloidal field suffers a sharp reversal of direction. According to the Ampère’s law, large electric currents must flow out of the plane shown in Figure 1, along the null surfaces, and in the presence of finite electric resistivity, dissipation of these currents will lead to reconnection of the oppositely directed field lines (e.g., Biskamp 1997, Lazarian & Vishniac 1999). We here investigate the possibility that the magnetic energy released by reconnection near the Y point region (see also Figure 2) is able to accelerate the plasma to relativistic velocities and produce the synchrotron radio flares (see below).

The release of magnetic energy by reconnection is believed to be the most promising mechanism for production of solar flares. Both direct and indirect astrophysical evidence show that the nature seems to know the way of producing fast reconnection in spite of the theoretical difficulties involved (see, e.g., Lazarian, Vishniac & Cho 2004). Y-type neutral zones are also present, for example, in eruptive magnetic arcades of the solar corona, like those associated to two-ribbon solar flares (see, e.g., Takasaki et al. 2004, Sakajiri et al. 2004). In these cases, a possible interpretation of the observations is that a current sheet forms between an emerging magnetic flux and the pre-existing magnetic flux of opposite polarization and magnetic reconnection in this zone will cause the formation of a filament/plasmoid. This will temporally inhibit further reconnection and establish a temporally stable magnetic configuration. During this stage, there will be considerable magnetic energy storage due to the continuing converging motion of the two fluxes. Eventually an eruption develops (for example, from sheared or twisted arcades lying underneath the magnetic null point), the filament/plasmoid is then ejected and a large quantity of magnetic energy is released (see, e.g., Takasaki et al. 2004).

Lazarian & Vishniac (1999) proposed a mechanism at which flares can naturally emerge as a result of an instability mediated by turbulence. In their model the reconnection rate increases due to magnetic field wandering caused by turbulence. Fluxes of non-parallel magnetic field lines reconnect very slowly if the turbulent velocities are much smaller than the Alfven velocity. However, outflows caused by reconnection increase the level of turbulence and therefore the reconnection rate. As a result of such a positive feedback a flare develops.

Whatever is the particular mechanism, the energy released by magnetic reconnection is expected to ultimately drive violent outward motions in the surrounding plasma that can produce copious amounts of solar cosmic rays via direct acceleration along the open magnetic field lines or via a first order Fermi process (e.g., Reames 1995). Similar mechanisms have been recently proposed to accelerate the plasma in magnetic reconnection zones of young born accretion-induced collapse pulsars to produce ultra-high energy cosmic rays (de Gouveia Dal Pino & Lazarian 2000, 2001).

The several stages described above for magnetic reconnection and solar flare production may be also operating in our accretion disk-corona-BH system, but since in the present study it is not crucial to follow the detailed evolution of the magnetic configuration and of the features that develop during the reconnection process, we will consider a simpler scenario. According to the discussion in the previous sections, we argue that whenever $`\beta `$ becomes small enough and the accretion rate attains values near the critical one in the inner disk region, then the lines of opposite polarization near the Y-zone will be pressed together by the accreted material sufficiently rapidly as to undergo violent reconnection events. Energy stored in the magnetic field will be released suddenly to heat the plasma and accelerate charged particles. Some of the high speed particles will travel downward to the BH or to the disk to the foot points of the magnetic-arches, yielding bright X-ray emission, others will spew outward giving rise to the relativistic radio ejecta (see below). <sup>1</sup><sup>1</sup>1In the scenario just described, the high accretion rate may also cause the compression and shear of magnetic loops that may exist in the lower parts of the inner disk corona, under the Y-zone dome. These could also undergo violent reconnection and even produce plasmoids that could eventually erupt and provoke particle acceleration as well, like in the Sun (see discussion in section 4).

We notice that the merging of the closed field lines of the loops and the open field lines in the disk further out from the inner radius can also lead to reconnection, however, we find that the amount of magnetic energy released in these regions is not sufficient to accelerate the particles to produce the observed radio luminosities in the large flares (see below).

For reconnection to occur in the inner disk region $`R_X`$, near the Y point, the poloidal flux from the disk and the flux anchored in the central magnetic structure surrounding the BH must have opposite signs, like in Figure 1. As pointed out by Tagger et al. (2004), this suggests that large ejections with violent flares will occur only when the system is in that configuration. Violent reconnection will not take place in the periods when both fluxes have the same sign. This dichotomy is observed for example in the interaction between the terrestrial and the solar wind magnetic field. Now, the magnetic flux trapped in the BH horizon is a sum of all magnetic flux accreted during the whole history of the source. If the magnetic flux in the disk is either produced by a turbulent dynamo in the disk itself or advected from a companion star, its sign is expected to change from time to time, although the characteristic time scale of this variation is difficult to compute in the absence of a detailed modelling of the dynamo in the disk or in the companion star. On the long term, a sufficient number of reconnection events might in turn cancel the flux trapped in the BH horizon and change its sign as well. In the case of the solar dynamo, it is not clear yet why the field reversals take place every 22 year. In the case of GRO 1915$`+`$105, Tagger et al. (2004) have suggested that this long term evolution of the field configuration could be of the order of one to few years so that part of the time the field fluxes in the inner disk edge and in the hole would be antiparallel and the rest of the time parallel.

As remarked above, observations of the Sun indicate that the amplification, dissipation and reconnection of the magnetic field lines is likely to occur not steadily, but sporadically, involving violent flares followed by more gradual recoveries to the pre-flare magnetic field configuration with no obvious periodicity. Similar behavior, i.e., violent flares followed by more slowly recovering of the pre-flare state is also detected in the observed superluminal radio ejecta of GRS 1915+105, thus suggesting that magnetic reconnection could be also associated to the production of relativistic blob ejections that rapidly move to the $`100`$ AU scales. In the next section, we try to quantify the ideas just described.

We notice that after reconnection, the partial destruction of the magnetic flux in the inner disk will make $`\beta `$ to increase and make the disk resumes the less magnetized condition of the plateau state with most of the energy being dissipated locally within the disk instead of in the outflow.

## 3 Building the Model for the Superluminal Ejections During the Flares

We adopt the field geometry of the magnetized accretion disk/corona as described in the previous section (Fig. 1).

We could first try to make an order-of-magnitude estimate of the amount of magnetic energy that could be extracted from the inner disk region by simply assuming that it is approximately equal to the disk accretion energy rate $`GM_{}\dot{M}/R_X`$. For a BH with mass $`M_{}=10M_{}=M_{10}`$ and Schwartzschild radius $`R_{}=2GM_{}/c^2=2.96\times 10^6M_{10}`$ cm, this results

$$\dot{W}_B\mathrm{\hspace{0.17em}10}^{39}\mathrm{erg}\mathrm{s}^1M_{10}\dot{M}_{19}R_{X,7}^1$$
(1)
where we have assumed that the inner radius of the accretion disk ($`R_X`$) corresponds approximately to the last stable orbit of the BH, $`R_X\mathrm{\hspace{0.17em}3}R_{}\mathrm{\hspace{0.17em}10}^7`$ cm, and that the accretion rate $`\dot{M}\dot{M}_{Edd}`$, where $`\dot{M}_{Edd}=1.9\times 10^{19}M_{10}`$ g s<sup>-1</sup> is the Eddington critical accretion rate (Shakura & Sunyaev 1973), $`\dot{M}_{19}`$ is $`\dot{M}`$ in units of $`10^{19}`$ g s<sup>-1</sup>, and $`R_{X,7}`$ is $`R_X`$ in units of $`10^7`$ cm. This estimate of $`\dot{W}_B`$ can be compared to the observed X-ray luminosity of GRS 1915+105 ($`10^{39}`$ erg s<sup>-1</sup>).

Now, in order to evaluate more precisely the amount of magnetic energy that can be extracted through violent reconnection from the Y neutral zone to produce superluminal ejecta, we need first to evaluate the physical conditions both in the magnetized disk and corona.

### 3.1 Disk and Coronal Parameters

To evaluate the disk quantities at $`R_X`$, we may employ the standard disk model (Shakura & Sunyaev 1973; Frank et al. 1992). In the inner disk region (named the X-region), the radiation pressure is larger than the thermal pressure and the disk will be $`radiation`$ pressure-dominated. In this case,

$$T_d\mathrm{\hspace{0.17em}1.2}\times 10^8\mathrm{K}\alpha _{0.1}^{1/4}M_{10}^{1/4}R_{X,7}^{3/4}$$
(2)
and

$$n_d\mathrm{\hspace{0.17em}2.8}\times 10^{17}\mathrm{cm}^3\alpha _{0.1}^1M_{10}^{1/2}\dot{M}_{19}^2R_{X,7}^{3/2}q_{0.82}^8$$
(3)
where, $`q=[1(R_{}/R_X)^{1/2}]^{1/4}=\mathrm{\hspace{0.17em}0.82}`$, $`T_d`$ and $`n_d`$ are the disk gas number density and temperature, respectively, $`\alpha _{0.1}`$ is the Shakura-Sunyaev viscous coefficient in units of 0.1 (e.g., Frank et al. 1992), and $`\dot{M}_{19}`$ is the disk accretion rate in units of $`10^{19}`$ g s<sup>-1</sup>.

The width $`H`$ of the disk in the radiation pressure-dominated region is given by

$$\frac{H}{R_X}\mathrm{\hspace{0.17em}0.7}\dot{M}_{19}R_{X,7}^1q_{0.82}^4$$
(4)
As remarked in the previous section, the violent magnetic reconnection events that may originate large scale relativistic ejections are expected to occur when the large scale magnetic field in the inner disk region, which is generated by dynamo process, attains an intensity as large as $`B_d^2/8\pi P_d`$, where $`P_d`$ is the disk pressure which is dominated by the radiation pressure, i.e., when $`\beta 1`$, giving

$$B_d\mathrm{\hspace{0.17em}6.9}\times 10^8\mathrm{G}\alpha _{0.1}^{1/2}M_{10}^{1/4}R_{X,7}^{3/4}$$
(5)
Recently, Liu, Mineshige, & Shibata (2002), have proposed a simple model to quantify the parameters of a magnetized disk-corona system. Assuming, like in the solar corona, that gas evaporation at the foot point of a magnetic flux tube quickly builds up the density of the corona to a certain value, and the tube radiates away the heating due to magnetic reconnection through Compton scattering, they derived the following coronal quantities (as a function of disk quantities) for a radiation pressure-dominated disk

$$n_c\mathrm{\hspace{0.17em}7.6}\times 10^{14}\mathrm{cm}^3\alpha _{0.1}^{15/16}\beta _1^{3/4}M_{10}^{25/32}\dot{M}_{19}^2R_{X,7}^{75/32}l_8^{3/4}$$
(6)
and

$$T_c\mathrm{\hspace{0.17em}4.9}\times 10^8\mathrm{K}\alpha _{0.1}^{15/32}\beta _1^{3/8}M_{10}^{25/64}\dot{M}_{19}^1R_{X,7}^{75/64}l_8^{1/8}$$
(7)
where $`n_c`$ and $`T_c`$ are the coronal gas number density and temperature, respectively, $`\beta _1`$ is the value of $`\beta `$ in the disk in units of 1, and $`l_8`$ is the scale height of the Y neutral zone in the corona in units of $`10^8`$ cm. Although the quantities above have been computed for the purposes of the present model only in the inner disk radius $`R_X`$, they are actually valid for any radius of the disk-corona within the radiation pressure-dominated region. We note that, according to the Liu et al.’s model, smaller values of $`\beta `$ would result more energy dissipation through the corona and larger coronal temperature and density, as indicated by Eqs. (6) and (7).

### 3.2 Rate of Magnetic Energy Release in the Y-type neutral zone

Let us assume that the poloidal magnetic field that raises in the corona just above the inner disk region ($`B_X`$) is of the order of the local disk magnetic field, $`B_d`$ (eq. 5), that has been produced and pushed to the inner disk region as described in the previous section (see also, e.g., Arons 1993, de Gouveia Dal Pino & Lazarian 2001). <sup>2</sup><sup>2</sup>2 We notice that this value is a little larger than the one estimated by assuming, e.g., that the ejected gas takes the bulk of the accretion power, i.e., $`L_jL_{acc}`$ (where $`L_jB_c^2R_X^2v_\varphi /2`$, $`L_{acc}\dot{M}v_\varphi ^2`$, and $`v_\varphi `$ is the Keplerian toroidal velocity in the disk; see also Livio et al. 2003).

The rate of magnetic energy that can be extracted from the Y-zone in the corona (above and below the disk) through reconnection is $`\dot{W}_B(B_X^2/8\pi )\xi v_A(4\pi R_X)Ł_X`$, where $`\xi =v_{rec}/v_A`$ is the reconnection efficiency factor, $`v_{rec}`$ is the reconnection velocity, $`v_A=B_X/(4\pi n_cm_p)^{1/2}`$ is the coronal Alfvén speed, $`m_p`$ is the hydrogen mass, and $`L_X`$ is the length of the reconnection region which according to Figure 3 is $`L_X/\mathrm{\Delta }R_Xv_A/\xi v_A\xi ^1`$. For the conditions here investigated we find that $`v_Ac`$. Substitution of these relations into the equation above gives

$$\dot{W}_B\frac{B_X^2}{8\pi }c(4\pi R_X)\mathrm{\Delta }R_X$$
(8)
In order to estimate $`\mathrm{\Delta }R_X`$ (or $`L_X`$), we need to model the reconnection process itself. Classical reconnection schemes use to provide too slow reconnection rates. However, there is ample observational evidence that reconnection in Astrophysics may take place at very large rates which are comparable with the Alfvén speed (see e.g., Dere 1996). In the framework of solar flares, considerable progress has been made in recent years towards the determination of more realistic reconnection rates during flare events based mostly on the separation speed of $`H\alpha `$ flare ribbons (e.g., Qiu et al. 2004; Jing et al. 2005; Sakajiri et al. 2004, Takasaki et al. 2004). These observations indicate that the reconnection speed can be as high as few tenths of the Alfvén velocity (Takasaki et al. 2004). Since in the conditions we deal with in the present analysis this speed approaches $`c`$, the corresponding reconnection rate can be rather large.

As mentioned earlier, in the turbulent reconnection model proposed by Lazarian & Vishniac (1999, henceforth LV99; 2000, see also Lazarian 2005 for an extensive review on this subject) wandering of magnetic field lines was suggested as the ultimate cause of fast reconnection. An alternative model (Bisckamp, Schwarz & Drake, 1997, Shay et al. 2003) attempts to explain fast reconnection appealing to plasma properties, such as anomalous resistivity (Parker 1979). These models do not exclude each other. For instance, while dealing with anomalous resistivity in most astrophysical situations one has also to deal with the effects of turbulence. At the same time the reconnection rates in the LV99 model may get an initial boost from the anomalous resistive effects. This boost may be important when the initial turbulent velocities are much smaller than the Alfven speed. In our present analysis, we focus our attention on regions of strong magnetic fields, like those around the BH in the X-region during nearly critical accretion rate. For those regions the resistivity should be anomalous and the fast reconnection should be guaranteed (with $`v_{rec}v_A`$).

Following LV99, we can estimate the width of the current sheet for which the resistivity should be anomalous:

$$\mathrm{\Delta }R_X=\frac{c\mathrm{\Delta }B}{4\pi n_cZev_{th,c}}$$
(9)
where $`\mathrm{\Delta }B2B_X`$ denotes the change of the magnetic field across the reconnection region, and $`v_{th,c}=(\frac{5/3kT_c}{m_p})^{1/2}`$ the thermal velocity of the ions of charge $`Ze`$ in the corona

$$v_{th,c}=\mathrm{\hspace{0.17em}2.65}\times 10^8\mathrm{cm}\mathrm{s}^1\alpha _{0.1}^{15/64}\beta _1^{3/16}M_{10}^{25/128}\dot{M}_{19}^{1/2}R_{X,7}^{75/128}l_8^{1/16}.$$
(10)
Substitution of Eqs. (5), (6), and (10) into Eq. (9) results

$$\left(\frac{\mathrm{\Delta }R_X}{R_X}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}0.0034}Z^1\alpha _{0.1}^{43/64}\beta _1^{15/16}M_{10}^{157/128}\dot{M}_{19}^{5/2}R_{X,7}^{599/128}l_8^{11/16}$$
(11)
Now, if we plug the result above and Eq. (5) into Eq. (8), we get

$$\dot{W}_B\mathrm{\hspace{0.17em}2.4}\times 10^{39}\mathrm{erg}\mathrm{s}^1\alpha _{0.1}^1M_{10}^{1/2}R_{X,7}^{1/2}\left(\frac{\mathrm{\Delta }R_X/R_X}{0.0034}\right)$$
(12)
and the corresponding reconnection time is

$$t_{rec}\frac{R_X}{\xi v_A}\mathrm{\hspace{0.17em}3.3}\times 10^4\mathrm{s}^1\xi ^1R_{X,7}$$
(13)
which indicates that the release of magnetic energy is very fast.

The analysis of Eqs. (11) and (12) shows that the adoption of a somewhat smaller viscosity parameter $`\alpha `$ (which would imply a more quiescent disk), by an order of magnitude, for instance, would require a disk accretion rate only few times larger in order to keep the values of $`\mathrm{\Delta }R_X/R_X`$ and $`\dot{W}_B`$ unaltered.

Eq. (12) gives the total expected amount of magnetic energy rate released by fast reconnection in the Y-zone during the flare state of GRS 1915+105. Part of this energy will heat the gas that will produce X-ray luminosity, and the remaining will accelerate the charged particles to relativistic velocities producing a violent radio ejecta. But what is the fraction that goes to each part, and what is the potential acceleration mechanism in the reconnection site that can produce the superluminal radio ejections?

The evaluation of the exact fraction of the magnetic energy that goes to each of the processes above is a task out of the scope of the present work, as it depends on the detailed modelling of the conditions of the disk-corona and of the efficiency of each process. However, assuming that most of the magnetic energy released through reconnection in the corona goes to heat the plasma, Liu et al. (2002) have computed the average value of the Compton parameter ($`y`$) for an optically thin corona with an original distribution of thermal electrons, $`y=\frac{4kT_c}{m_ec^2}n_c\sigma _Tl`$, which was weighted over all the disk radii. According to their results, for a disk near the critical accretion ($`\dot{M}\dot{M}_{Edd}`$), $`y<\mathrm{\hspace{0.17em}0.8}`$, which implies a soft (and steep) X-ray spectrum, just as required by GRS 1915+105 observations during the flare state, with a spectral index $`\alpha _{rx}=\left(\frac{9}{4}+\frac{4}{y}\right)^{1/2}\frac{3}{2}>\mathrm{\hspace{0.17em}1.2}`$. This result is qualitatively consistent with the present model, as it predicts, for a magnetized disk-corona a predominantly soft X-ray spectrum during nearly critical accretion.

Since observations do indicate that both, heating and acceleration of the plasma to relativistic velocities are occurring, as suggested by the observed disk X-ray luminosity of the order of $`10^{39}`$ erg s<sup>-1</sup>, and by the radio flares which in turn suggest a jet mechanical power also of the order of $`10^{39}`$ erg s<sup>-1</sup> at least in the nuclear region (see, e.g, Falcke & Biermann 1999), we will here assume that a substantial fraction of $`\dot{W}_B`$ (Eq. 12) goes to accelerate the particles.

Let us now try to answer the second part of the question above, i.e., what could be the acceleration mechanism in the reconnection region able to produce the relativistic ejecta and the radio flare?

## 4  The Acceleration Mechanism during the Flare State

As stressed by de Gouveia Dal Pino & Lazarian (2001), the particular mechanism of particle acceleration during reconnection events is still unclear in spite of numerous attempts to solve the problem (see LaRosa et al 1996, Litvinenko 1996). Cosmic rays from the Sun confirm that the process is sufficiently efficient in spite of the apparent theoretical difficulties for its explanation.

Particles could, in principle, be accelerated at once to the relativistic energies by the large induced electric field within the reconnection region (e.g., Bruhwiller & Zweibel 1992, Haswell, Tajima, & Sakai 1992, Litvinenko 1996; de Gouveia Dal Pino & Lazarian 2000, 2001), but one-shot acceleration mechanism can make it difficulty to explain how particles will afterwards produce a power-law particle and synchrotron spectrum. Another possibility is that particles are accelerated by a Fermi process in the reconnection site. In this regard, we discuss below two viable sorts of this process.

A schematic representation of a reconnection region is shown in Figure 3. The upper and lower parts of the magnetic flux move towards each other with the velocity $`v_{rec}`$. As a result, charged particles rays in the upper part of the reconnection zone ”see” the lower part of the magnetic flux to approach them with a velocity $`2v_{rec}`$. In this case an acceleration process analogous to the first-order Fermi acceleration of cosmic rays in magnetized shocks may take place (e.g., Longair 1992) and we can write

$$\frac{\mathrm{\Delta }E}{E}=\frac{2v_{rec}}{c}\mathrm{cos}\theta $$
(14)
For a particle probability distribution $`p(\theta )=2\mathrm{sin}\theta \mathrm{cos}\theta d\theta `$, the average energy gain per crossing of the reconnection region is <sup>3</sup><sup>3</sup>3 We note that since the particle Larmour radius around the magnetic field in the reconnection zone, $`r_B=\gamma m_ec^2/(ZeB_X)\mathrm{\hspace{0.17em}9}\times 10^6`$ cm $`\gamma _{3.5}B_{X,7}<<\mathrm{\Delta }R_X`$ (where $`\gamma _{3.5}=\gamma /3.5`$ is the Lorentz factor for electrons with velocities $`0.92c`$), the particles will remain confined time enough in the reconnection zone to be efficiently accelerated through the Fermi mechanism before escaping from the system.

$$\frac{\mathrm{\Delta }E}{E}=\frac{v_{rec}}{c}_0^{\pi /2}2\mathrm{cos}^2\theta \mathrm{sin}\theta d\theta =\frac{4}{3}\frac{v_{rec}}{c}$$
(15)
where $`\theta `$ is the pitch angle between the particle velocity and the magnetic field. The round trip of the particles between the upper and lower magnetic fluxes will produce an increment of energy twice as large. Thus the ratio of the particle energy, $`E`$, after such a round trip to its original energy, $`E_0`$, will be

$$\delta =\frac{E}{E_0}=1+\frac{8}{3}\frac{v_{rec}}{c}$$
(16)
Now, assuming that the escape probability of a particle from the acceleration zone is similar to the one computed for a shock front, $`4v_{rec}/c`$, then probability to stay within the reconnection region is simply

$$P=1\frac{4v_{rec}}{c},$$
(17)
and the electron particle spectrum is obtained from (see e.g., Longair 1992)

$$N(E)dEE^{1+\frac{lnP}{ln\delta }}dEE^{5/2}dE$$
(18)
This equation indicates that the particle spectrum produced by first-order Fermi acceleration within the reconnection site is steeper than that normally produced in shocks (for which a similar simplified treatment that ignores losses and non-linearities provides $`N(E)E^2`$), so that the two processes may be, in principle, distinguished. Now, the Synchrotron radiation from electrons with a power-law distribution $`N(E)E^p`$, has a power-law frequency radio spectrum $`S_\nu \nu ^{\alpha _r}`$, where the spectral index is given by $`\alpha _r=(p\mathrm{\hspace{0.17em}1})/2\mathrm{\hspace{0.17em}0.75}`$. Though crude, this estimated value is consistent with the observed radio spectrum during the flares ($`\alpha _r0.6`$).

Finally, the particle acceleration mechanism discussed above does not disregard the possibility that the relativistic fluid may be also produced behind shocks. As in the Sun, plasmoids formed by reconnection of the field lines lines near the Y-zone may violently erupt and cause the formation of a shock front. This raises the alternative possibility of a first-order Fermi acceleration process of the ionized gas behind the shock front and the development of a power-law spectrum with particle spectral index of the order of two ($`N(E)E^2`$) as remarked above and, as such, a corresponding synchrotron spectrum with a spectral index also of the order of the observed one, $`S_\nu \nu ^{0.5}`$.

After a radio flare, while the system is slowly recovering its original magnetic field configuration, the plateau phase must resume. During this recovering phase, the magnetic field intensity in the X-region will be smaller and also the mass accretion rate. In this state, the X-ray spectrum is observed to be hard and flat. As we have seen in ¶3.2, the model of Lin et al. (2002) predicts a harder X-ray spectrum when the accretion rate in the magnetized disk is very sub-Eddington for which a spectral index $`\alpha _{rx}<\mathrm{\hspace{0.17em}1.2}`$ should be expected.

## 5 Conclusions and Discussion

We have here investigated the origin of the large scale superluminal ejections observed in the galactic microquasar GRS 1915+105 during radio flare events and proposed a very simple model in which they could be due to violent magnetic reconnection episodes in the corona just above the inner edge of the magnetized accretion disk that surrounds the central $`10M_{}`$ black hole. In this region, when a large scale magnetic field is built by a turbulent dynamo process, and there is approximate balance between the disk magnetic and gas$`+`$radiation pressures ($`\beta 1`$), which implies a strong magnetic field intensity $`B_X7\times 10^8`$ G, reconnection between the magnetic lines that arise from the disk with those of the magnetosphere of the BH may become violent. During this regime, the vertical magnetic flux generates a wind that will substantially remove angular momentum from the disk therefore increasing the disk mass accretion to a value near the critical one (with $`\dot{M}\mathrm{\hspace{0.17em}10}^{19}`$ g s<sup>-1</sup>).

Whenever the system reaches the condition above and the magnetic field lines of the BH magnetosphere and those of the disk have opposite polarization, it may undergo an episode of violent reconnection that lasts for very short time $`3\times 10^4`$ s. A substantial part of the total magnetic energy released by reconnection ($`2\times 10^{39}`$ erg s<sup>-1</sup>) heats the coronal gas to $`T_c5\times 10^8`$ K that produces a steep, soft X-ray spectrum with the expected luminosity from observations. Assuming that the remaining released magnetic energy goes to accelerate the particles to relativistic velocities ($`vv_Ac`$, where $`v_A`$ is the Alfvén speed) in the reconnection site through first-order Fermi processes, we have found that a power-law electron distribution $`N(E)E^{5/2,2}`$, and a synchrotron radio power-law spectrum ($`S_\nu \nu ^{0.75,0.5}`$) are produced with radio spectral indices which are comparable to the one observed during the flares.

After each violent reconnection event, the system must resume a sub-critical disk mass accretion with smaller magnetic field intensity, and a slower reconnection rate, thus recovering the more quiescent plateau state, as described in ¶2 to 4.

Although until the present, GRS 1915 + 105 has been the only system among the galactic accreting black holes where large scale ($``$ 100 AU) superluminal ejections have been observed (maybe because this system is the only one that meets the necessary conditions as investigated here), the magnetized-disk-corona model here discussed could be applicable, in principle, also to other X-ray binary sources, or even to the microquasar extragalactic counterparts, i.e., quasars with evidence of superluminal ejection.

###### Acknowledgements.

This paper has benefited from many valuable comments from an anonymous referee. E.M.G.D.P. has been partially supported by grants of the Brazilian Agencies FAPESP (grant 1997/13084-3) and CNPq.


format: text

text:
# Numerical Analysis of the Capacities for Two-Qubit Unitary Operations

## I Introduction

A nonlocal operation is one which operates on distinct subsystems, and can not be expressed as a tensor product of operations on the individual subsystems. Such operations arise via an interaction Hamiltonian between the subsystems. Nonlocal operations may be used to create entanglement between subsystems, and also to perform classical communication. In fact, it is not possible to achieve either of these tasks without an interaction between the subsystems. In characterizing quantum operations, it is thus important to determine the capacities for creating entanglement or performing communication. These capacities provide a characterization of the strength of the operations nielsen .

Shared classical information may be considered to be the classical equivalent of entanglement. Therefore it is reasonable to consider the process of classical communication to be the classical equivalent of entanglement creation. One may therefore expect that there is a close relationship between the capacities for these two tasks. In fact, it has been shown bhls ; berry1 ; berry2  that there are a number of inequalities between the various capacities for classical communication and entanglement creation. In this paper we make a direct numerical comparison between these capacities.

The capacities of operations for creating entanglement have been studied extensively cirac ; zanardi ; durvid ; kradur ; kraus ; leifer ; childs ; kraus2 . It is relatively straightforward to determine the entanglement capability for infinitesimal operations durvid . For finite operations, most results are restricted to numerical results for certain classes of two-qubit operations leifer . We study the same classes of operations here, and the results we present for the entanglement reproduce those given in Ref. leifer , except for some data points where we believe our results are more accurate.

The capacities for classical communication were initially considered for simple operations, such as the CNOT and SWAP operations eisert ; collins . It is more difficult to consider the classical communication for general operations, because general operations can not be used to perform perfect communication. Ref. bhls  introduced asymptotic capacities, where the average communication when the operation is performed a very large number of times is considered. When the operation is performed a large number of times, it is possible to use error correcting techniques to reduce the probability of error to be arbitrarily small.

For the purposes of numerical investigation, it is not practical to use asymptotic capacities, as the number of variables that would need to be optimized over to obtain a reasonable approximation is extremely large. A far more practical type of capacity is that based on the Holevo information of ensembles. For these capacities, only the states and probabilities in the initial ensemble need be optimized over in the numerical analysis. In addition, there are connections between the capacities defined via the change in Holevo information and the asymptotic capacities bhls .

Although it is feasible to calculate capacities based on the Holevo information, it is far more computationally difficult than evaluating the entanglement capacities. It would therefore be useful if it were possible to use the entanglement capacities to estimate the capacities based on the Holevo information. In this paper we numerically study the relationship between these capacities, and show that although they are not equal, the values we have calculated are quite close.

In Sec. II we explain the definitions of the various capacities. In Sec. III we give numerical results for the communication capacity based on the Holevo information obtained when the initial ensemble has zero Holevo information, and compare these capacities to the entanglement that may be created from initial states that have zero entanglement. Then in Sec. IV we give results for the increase in Holevo information for general initial ensembles, and compare these capacities to the increase in entanglement for arbitrary initially entangled states. We conclude in Sec. V, and in the Appendix give a detailed explanation of the numerical techniques used to calculate the results presented.

## II Definitions

First we summarize the definitions of the various capacities. Throughout this paper we divide the system into two subsystems, $`A`$ and $`B`$, and denote the Hilbert spaces by $`_A`$ and $`_B`$. The party in possession of subsystem $`A`$ will be referred to as Alice and the party in possession of subsystem $`B`$ will be referred to as Bob. The subsystems $`A`$ and $`B`$ are divided into further subsystems:

$$_A=_{A_U}_{A_{\mathrm{anc}}},_B=_{B_U}_{B_{\mathrm{anc}}}.$$
(1)
The operation $`U`$ acts only upon $`_{A_U}_{B_U}`$, and the Hilbert spaces $`_{A_{\mathrm{anc}}}`$ and $`_{B_{\mathrm{anc}}}`$ are ancillas. The ancillas may have dimension that is arbitrarily large but finite.

There are two main ways of defining capacities for entanglement. The first is the entanglement that may be obtained when the initial state is unentangled 111Subsystems A and B are not entangled, but there may be entanglement within the respective subsystems.

$$E_U\underset{|\varphi _A_A,|\chi _B_B}{sup}E(U|\varphi _A|\chi _B).$$
(2)
The quantity $`E(\mathrm{})`$ is the entropy of entanglement $`E(|\mathrm{\Psi })=S[\mathrm{Tr}_A(|\mathrm{\Psi }\mathrm{\Psi }|)]`$, where $`S(\rho )=\mathrm{Tr}(\rho \mathrm{log}\rho )`$. Throughout we employ logarithms to base 2, so the entanglement is expressed in units of ebits. The second definition is the maximum increase in entanglement when the initial state may be an arbitrary pure entangled state.

$$\mathrm{\Delta }E_U\underset{|\mathrm{\Psi }_A_B}{sup}\left[E(U|\mathrm{\Psi })E(|\mathrm{\Psi })\right].$$
(3)
It is also possible to define the capacity in terms of the entanglement of formation and allow mixed states; we do not separately consider this case, because allowing mixed states does not alter the capacity bhls .

Another alternative definition of the entanglement capacity is based on the average entanglement that may be obtained in the limit that the operation is performed an extremely large number of times and the initial state is unentangled bhls . It is shown in Ref. bhls  that this capacity is equal to the maximum increase in entanglement as defined in Eq. (3). Thus the numerical results presented also apply to the asymptotic entanglement capacity.

In the numerical search, it is not possible to consider ancilla spaces with arbitrarily large dimension. It is known nielsen  that the ancilla spaces need have dimension no larger than the Hilbert spaces $`_{A_U}`$ and $`_{B_U}`$ for the capacity with initially unentangled states. It has also been found numerically leifer  that the same is true for the capacity where initially entangled states are allowed.

In the results presented below, we use a fixed dimension for the ancilla spaces. For simplicity we use equal dimensions on the two ancilla spaces $`_{A_{\mathrm{anc}}}`$ and $`_{B_{\mathrm{anc}}}`$. We use a superscript on the capacity when it is necessary to indicate the dimension of the ancilla space used. For example, $`\mathrm{\Delta }E_U^{(4)}`$ is the maximum change in entanglement when the ancilla spaces are each of dimension 4. When we refer to multiple results with different ancilla dimensions we will use a superscript $`()`$. We will omit the superscript in the case of $`E_U`$ when the dimensions of the ancilla spaces are at least as large as those of $`_{A_U}`$ and $`_{B_U}`$, because this is known to be sufficient to obtain the capacity for arbitrarily large ancilla.

The classical communication capacities that we will consider are based upon the Holevo information of ensembles. An ensemble is a set of states $`\{|\mathrm{\Phi }_i_{AB}\}`$ that are supplied with probabilities $`p_i`$. Each state $`|\mathrm{\Phi }_i_{AB}`$ is a pure state shared between Alice and Bob, and Alice chooses the index $`i`$. The ensemble is denoted by $`=\{p_i,|\mathrm{\Phi }_i_{AB}\}`$. We also define the ensemble of reduced density matrices possessed by Bob as

$$𝖤=\mathrm{Tr}_A=\{p_i,\rho _i\},$$
(4)
where $`\rho _i=\mathrm{Tr}_A|\mathrm{\Phi }_i_{AB}\mathrm{\Phi }_i|`$. The Holevo information of the ensemble $`𝖤`$ is given by

$$\chi (𝖤)=S\left(\underset{i}{}p_i\rho _i\right)\underset{i}{}p_iS(\rho _i).$$
(5)
From the Holevo-Schumacher-Westmoreland theorem holevo ; shuwes , the Holevo information gives the average communication that may be performed from Alice to Bob by coding over multiple states.

Similarly to the case for entanglement, we may define capacities based on the maximum change in Holevo information. One definition that we will use is the maximum final Holevo information when the initial ensemble has zero Holevo information. For the initial ensemble, we have an initial state $`|\psi _{AB}`$, and Alice encodes $`i`$ by applying a local operation $`V^{(i)}`$. For the capacity, the supremum is taken over the initial state $`|\psi _{AB}`$, the encoding operations $`V^{(i)}`$ and the probabilities $`p_i`$:

$$\chi _U=\underset{p_i,V^{(i)},|\psi _{AB}}{sup}\chi (p_i,\mathrm{Tr}_AUV^{(i)}|\psi _{AB}),$$
(6)
where $`\mathrm{Tr}_A|\varphi \mathrm{Tr}_A|\varphi \varphi |`$. Note that the notation we are using here differs from Ref. bhls , where the symbol $`\mathrm{\Delta }\chi _U^{(1,\mathrm{})}`$ was used for this capacity.

We may also define the maximum change in Holevo information when the initial ensemble is arbitrary:

$$\mathrm{\Delta }\chi _U=\underset{}{sup}\chi (\mathrm{Tr}_AU)\chi (\mathrm{Tr}_A).$$
(7)
Here we are using the notation conventions

$`U`$ $`\{p_i,U|\mathrm{\Phi }_i\},`$ (8)
$`\mathrm{Tr}_X`$ $`\{p_i,\mathrm{Tr}_X(|\mathrm{\Phi }_i)\}.`$ (9)
This capacity is equivalent to the capacity $`\mathrm{\Delta }\chi _U^{(1,)}`$ defined in Ref. bhls . As shown in Ref. bhls , this capacity is equal to the average entanglement-assisted communication that may be performed from Alice to Bob. Therefore this quantity may be interpreted as the asymptotic communication capacity, just as $`\mathrm{\Delta }E_U`$ may be interpreted as the asymptotic entanglement capability.

One may also interpret $`\chi _U`$ in terms of asymptotic capacities. The capacity $`\chi _U`$ gives the Holevo information after a single application of the operation $`U`$. This communication can not actually be performed for a single ensemble; it is necessary to code over multiple states to perform this average communication. Therefore $`\chi _U`$ may be interpreted as the asymptotic communication capacity if the communication protocol is limited to the relatively simple scheme where coding is performed over multiple final states. This is equivalent to restricting all the applications of $`U`$ to be performed at the same time (on input states that are not entangled with each other), rather than allowing the output of one application of $`U`$ to be used as part of the input to another application of $`U`$, as in the general case.

As in the case of the entanglement, it is not possible to use ancilla spaces of arbitrarily large dimension in the numerical search. In addition, it is not possible to use arbitrarily large numbers of states in the ensemble. In the results presented below, we perform calculations for restricted ensembles where there is a fixed number of states in the ensemble and a fixed dimension for the ancilla spaces (the ancilla spaces are again taken to be of equal dimension). We use superscripts on the capacities to indicate the number of states in the ensemble and the dimension of the ancilla spaces. For example, $`\mathrm{\Delta }\chi _U^{(2,4)}`$ is the maximum change in Holevo information for two states in the ensemble and ancillas each of dimension 4. We use a superscript $`()`$ to refer to multiple capacities with different ancilla dimensions or ensemble sizes. It must be emphasized that our use of superscripts in this paper differs from that in Ref. bhls .

## III Capacities for zero initial Holevo information

It is clear that the capacity $`\chi _U`$ is an analogous quantity for communication to $`E_U`$ for entanglement. Similarly $`\mathrm{\Delta }\chi _U`$ is analogous to $`\mathrm{\Delta }E_U`$ for entanglement. In this section we perform a direct numerical comparison between the two capacities $`\chi _U`$ and $`E_U`$. In the next section we compare the capacities $`\mathrm{\Delta }\chi _U`$ and $`\mathrm{\Delta }E_U`$.

In this paper we concentrate on two-qubit operations. It is not possible to perform calculations for the entire range of two-qubit operations. To make the problem feasible, we only consider a limited number of examples of two-qubit operations. In particular, we consider operations of the form:

$`U_1(\alpha )`$ $`=U_d(\alpha ,0,0),`$ (10)
$`U_2(\alpha )`$ $`=U_d(\alpha ,\alpha ,0),`$ (11)
$`U_3(\alpha )`$ $`=U_d(\alpha ,\alpha ,\alpha ),`$ (12)
where

$$U_d(\alpha _1,\alpha _2,\alpha _3)=\mathrm{e}^{i(\alpha _1\sigma _1\sigma _1+\alpha _2\sigma _2\sigma _2+\alpha _3\sigma _3\sigma _3)}.$$
(13)
Here $`\sigma _i`$ are the Pauli sigma operators. The operations $`U_1`$, $`U_2`$ and $`U_3`$ correspond to the CNOT, DCNOT and SWAP families of operations considered in Ref. leifer .

In order to consider the complete range of two-qubit operations in the case of entanglement, it would only be necessary to consider operations of the form (13), with $`\pi /4\alpha _1\alpha _2\alpha _30`$ kraus . This derivation relies on the fact that any two-qubit operation may be simplified to one of the form (13) with $`\pi /4\alpha _1\pm \alpha _2\alpha _30`$ using local operations durvid ; kraus ; makhlin ; hammer . In addition to using local operations, the derivation in Ref. kraus  relies on the fact that the entanglement capabilities of $`U`$ and $`U^{}`$ are identical (which implies that all the $`\alpha _i`$ may be taken to be positive).

Similarly, for the Holevo information, $`\chi (\mathrm{Tr}_A)=\chi (\mathrm{Tr}_A^{})`$ and $`\chi (\mathrm{Tr}_AU)=\chi (\mathrm{Tr}_AU^{}^{})`$. Thus the capacities of $`U`$ and $`U^{}`$ to increase the Holevo information are identical. Therefore, in order to obtain results for the complete range of two-qubit operations, it would only be necessary to consider operations of the form (13) with $`\pi /4\alpha _1\alpha _2\alpha _30`$ for both entanglement and Holevo information. This restriction on the values of the $`\alpha _i`$ defines a three dimensional region of values. In this paper we do not consider the entire region; however, the operations $`U_1`$, $`U_2`$ and $`U_3`$ which we consider form three lines on the boundaries of this region.

Sophisticated numerical maximization techniques were used to find the values of $`E_U`$ and $`\chi _U^{()}`$; the details are given in the Appendix. The numerical results for $`E_{U_1}`$ and $`\chi _{U_1}^{()}`$ are shown in Fig. 1. The values of $`E_{U_1}`$ were determined in intervals of $`\pi /400`$ for $`\alpha `$, and are shown by the solid line. An ancilla dimension of 2 is sufficient to obtain the asymptotic capacity $`E_U`$ for all two-qubit operations nielsen . In all the calculations presented for $`E_U`$ in this section, an ancilla dimension of 2 was used; we have therefore omitted the superscript on this capacity. In addition, it was found that, for $`U_1`$, the same capacity was obtained without ancilla, in agreement with Ref. kraus .

The values of $`\chi _{U_1}^{()}`$ were determined in intervals of $`\pi /40`$ for $`\alpha `$, with ancilla dimensions of 1, 2 and 4, and ensembles with 2 or 4 members. An ancilla dimension of 1 is equivalent to no ancilla; we will take the case of no ancilla to be an ancilla dimension of 1 for convenience in the notation. Just as in the case of the entanglement, these results indicate that the maximum final Holevo information is obtained without ancillas. In addition, only two states in the ensemble are required. We find no improvement in using up to four states in the ensemble, and ancillas of dimension up to four. It was found that the maximum increase in the Holevo information was obtained with equal probabilities when there were two states in the ensemble. When four states in the ensemble were used, in most cases two of the probabilities approached zero whereas the other two approached $`1/2`$.

In addition, note that there is no difference between the results obtained for the maximum final entanglement and the maximum final Holevo information. It was found that the two capacities agreed to precision better than one part in $`10^{14}`$. Our results do not prove that $`E_{U_1}=\chi _{U_1}`$, because it is not possible to test ensembles with arbitrarily large numbers of states or allow the ancilla dimension to be arbitrarily large. Nevertheless, our results strongly indicate that $`E_{U_1}`$ is equal to $`\chi _{U_1}`$.

Numerical results for $`E_{U_2}`$ and $`\chi _{U_2}^{()}`$ are shown in Fig. 2. As for $`E_{U_1}`$, the value of $`E_{U_2}`$ was determined in intervals of $`\pi /400`$ for $`\alpha `$, and an ancilla dimension of 2 was used in order to obtain the asymptotic entanglement capability $`E_{U_2}`$. The values of $`\chi _{U_2}^{()}`$ were determined in intervals of $`\pi /40`$ for $`\alpha `$, and with ancilla dimensions of 2, 3, 4 and 5. In each case the ensemble consists of four states.

In contrast to the case for the entanglement, there are improvements in using ancilla dimensions higher than 2 when considering the Holevo information. For an ancilla dimension of 2, the value of $`\chi _{U_2}^{(4,2)}`$ is equal to $`E_{U_2}`$ for the lower values of $`\alpha `$, but is less than $`E_{U_2}`$ for the upper values of $`\alpha `$. For an ancilla dimension of 3, the values of $`\chi _{U_2}^{(4,3)}`$ are larger, but still less than $`E_{U_2}`$ towards the right side of the plot. For an ancilla dimension of 4, $`\chi _{U_2}^{(4,4)}`$ is equal to $`E_{U_2}`$. When the ancilla dimension is increased to 5, the capacity $`\chi _{U_2}^{(4,5)}`$ is no larger than $`\chi _{U_2}^{(4,4)}`$. These results indicate that an ancilla dimension of 4 may be sufficient to achieve the capacity $`\chi _{U_2}`$, and that this capacity is equal to $`E_{U_2}`$.

Similar calculations were performed with equal probabilities for each of the states in the ensemble. In most cases the same results were obtained as when arbitrary probabilities were allowed. However, for one data point for $`\chi _{U_2}^{(4,2)}`$ (indicated by the circle in Fig. 2), a larger Holevo information was obtained when unequal probabilities were allowed. The numerically derived optimal ensemble had three members with non-negligible probabilities: two probabilities were approximately $`41.39\%`$, and one was $`17.22\%`$.

Results for $`E_{U_3}`$ and $`\chi _{U_3}^{()}`$ are shown in Fig. 3. Again the value of $`E_{U_3}`$ was determined in intervals of $`\pi /400`$ for $`\alpha `$, and with an ancilla dimension of 2. The values of $`\chi _{U_3}`$ were determined for ancilla dimensions of 2, 3, 4 and 5, and with four states in the ensemble.

As in the case of $`\chi _{U_2}^{(4,2)}`$, the values of $`\chi _{U_3}^{(4,2)}`$ are significantly below $`E_{U_3}`$ for most of the plot. When the ancilla dimension is increased to 3, the values of $`\chi _{U_3}^{(4,3)}`$ are noticeably larger, and closer to $`E_{U_3}`$. In addition, three of the values, at $`\alpha =\pi /40`$, $`2\pi /40`$ and $`3\pi /40`$, are actually larger than $`E_{U_3}`$ (see the subplot in Fig. 3). The difference is small, less than $`0.02`$, but it is sufficient to demonstrate that $`E_{U_3}`$ is not equal to $`\chi _{U_3}`$. For an ancilla dimension of 4, $`\chi _{U_3}^{(4,4)}`$ is not smaller than $`E_{U_3}`$ for any of the data points. In addition, $`\chi _{U_3}^{(4,4)}`$ is larger than $`E_{U_3}`$ for the same three values of $`\alpha `$ as for an ancilla dimension of 3. In fact, $`\chi _{U_3}^{(4,4)}`$ is equal to $`\chi _{U_3}^{(4,3)}`$ for these three data points. When the ancilla dimension is increased to 5, $`\chi _{U_3}^{(4,5)}`$ is no larger than $`\chi _{U_3}^{(4,4)}`$, again indicating that an ancilla dimension of 4 is sufficient to achieve the capacity $`\chi _{U_3}`$.

To summarize, we have found that, for the operation $`U_1`$, ancillas do not increase the capacity $`\chi _U`$, just as in the case for the entanglement. For the cases $`U_2`$ and $`U_3`$, the capacity $`\chi _U`$ increases with ancilla dimension up to an ancilla dimension of 4, but is unchanged when the ancilla dimension is further increased to 5. This is in contrast to the case for the entanglement, where the maximum entanglement capacity $`E_U`$ is obtained for an ancilla dimension of 2 nielsen .

The results indicate that $`\chi _U`$ is equal to $`E_U`$ for the operations $`U_1`$ and $`U_2`$, though it is possible that $`\chi _U^{()}`$ is increased for larger ancilla dimensions or ensemble sizes. For the operation $`U_3`$, it is possible to obtain slightly higher values of $`\chi _U^{()}`$, demonstrating that $`\chi _U`$ is not equal to $`E_U`$ for this operation. In all cases tested, we find that $`\chi _UE_U`$. In only a small number of cases have we found that $`\chi _UE_U`$, and in these cases the differences found are only small, suggesting that $`E_U`$ is an excellent, and efficient, estimator of $`\chi _U`$.

## IV Capacities for arbitrary initial ensembles

Next we consider the capacities $`\mathrm{\Delta }\chi _U`$ and $`\mathrm{\Delta }E_U`$. These capacities are more general, in that arbitrary initial states or ensembles are allowed. Analytic results for the relation between these capacities have been derived in Ref. berry2 . In this reference it is shown that, for two-qubit operations, $`\mathrm{\Delta }\chi _U\mathrm{\Delta }E_U`$. If $`\mathrm{\Delta }E_U`$ may be achieved with a particular ancilla dimension, then an increase in Holevo information equal to $`\mathrm{\Delta }E_U`$ may be achieved with the same ancilla dimension, and with four members in the ensemble.

In principle it is possible that there is no ancilla dimension that achieves $`\mathrm{\Delta }E_U`$, and instead $`\mathrm{\Delta }E_U`$ is approached in the limit of large ancilla dimension. However, in practice it has been found that, for two-qubit operations, it appears to be possible to achieve $`\mathrm{\Delta }E_U`$ with an ancilla dimension of 2 leifer . This means that it should be possible to achieve an increase in Holevo information equal to $`\mathrm{\Delta }E_U`$ with an ancilla dimension of 2.

The capacities $`\mathrm{\Delta }\chi _{U_1}^{()}`$ and $`\mathrm{\Delta }E_{U_1}^{()}`$ are shown in Fig. 4. Each capacity was determined in steps of $`\pi /40`$ in $`\alpha `$. It was found that the entanglement capacity $`\mathrm{\Delta }E_{U_1}^{()}`$ did not increase beyond that for no ancilla as the ancilla dimension was increased up to 5, in agreement with the result given in Ref. leifer . Throughout this section we use the superscript asterisk when referring to results for $`\mathrm{\Delta }E_U^{()}`$, because it has not been proven that the asymptotic capacity is achieved for an ancilla dimension of 2.

The capacities $`\mathrm{\Delta }\chi _{U_1}^{()}`$ without ancillas and with ancillas of dimension 2 are shown in Fig. 4. In both cases these capacities are for ensembles with two states. It is found that, even without ancilla, the capacity $`\mathrm{\Delta }\chi _{U_1}^{()}`$ is greater than the values calculated for $`\mathrm{\Delta }E_{U_1}^{()}`$. The only cases where there is equality are the trivial cases where $`\alpha =0`$ or $`\pi /4`$. These results indicate that there are operations for which there is the strict inequality $`\mathrm{\Delta }\chi _U>\mathrm{\Delta }E_U`$.

In addition, the capacity $`\mathrm{\Delta }\chi _{U_1}^{()}`$ is slightly increased by adding an ancilla. This is not so easily seen in Fig. 4; to make this difference visible, the differences between the capacities $`\mathrm{\Delta }\chi _{U_1}^{()}`$ with ancilla and the capacity with no ancilla $`\mathrm{\Delta }\chi _{U_1}^{(2,1)}`$ are plotted in Fig. 5. It can be seen that there is a small but significant increase in $`\mathrm{\Delta }\chi _{U_1}^{()}`$ when an ancilla is allowed.

In addition, there is a further improvement in using an ancilla dimension of 3 rather than an ancilla dimension of 2. There are further increases in the capacity as the ancilla dimension is increased to 4 and 5, as shown in Fig. 6. These results indicate that the true asymptotic capacity $`\mathrm{\Delta }\chi _{U_1}`$ is not actually achieved for any particular ancilla dimension. However, each increase in the capacity with the ancilla dimension is smaller than the previous, indicating that the results calculated here should give a good approximation of $`\mathrm{\Delta }\chi _{U_1}`$.

Calculations were also performed for ensembles with 4 members, and ancilla dimensions up to 4. It was found that, for all ancilla dimensions tested, there was no increase in the capacity when the number of states in the ancilla was increased. In addition, it was found that there was no increase in the capacity with up to 8 states in the ensemble and no ancilla.

The results for $`\mathrm{\Delta }\chi _{U_2}^{()}`$ and $`\mathrm{\Delta }E_{U_2}^{()}`$ are shown in Fig. 7. In each case shown, ensembles with four states were used. In the case without ancilla, it was found that $`\mathrm{\Delta }\chi _{U_2}^{(4,1)}`$ and $`\mathrm{\Delta }E_{U_2}^{(1)}`$ were equal. When an ancilla is included, there is a significant increase in both $`\mathrm{\Delta }\chi _{U_2}^{()}`$ and $`\mathrm{\Delta }E_{U_2}^{()}`$. In particular, these have a maximum of 2, rather than 1 as in the case without ancilla.

Note also that the value of $`\mathrm{\Delta }E_{U_2}^{()}`$ is increased when an ancilla is added for each of the values of $`\alpha `$ except the trivial points at $`\alpha =0`$ and $`\pi /4`$. In contrast, for the data shown in Ref. leifer  there was no visible increase in $`\mathrm{\Delta }E_{U_2}^{()}`$ when the ancilla was included for another three data points (at $`\alpha =\pi /40`$, $`2\pi /40`$ and $`3\pi /40`$). We suspect that this is because the optimization found the local maximum corresponding to the solution with no ancilla, rather than the global maximum.

We found that using ancilla dimensions above 2 up to an ancilla dimension of 5 did not increase $`\mathrm{\Delta }E_{U_2}^{()}`$, in agreement with Ref. leifer . When the ancilla was included, $`\mathrm{\Delta }\chi _{U_2}^{()}`$ was slightly greater than $`\mathrm{\Delta }E_{U_2}^{()}`$, just as in the case of the operation $`U_1`$. In addition, it was found that $`\mathrm{\Delta }\chi _{U_2}^{()}`$ was further increased as the ancilla dimension was increased beyond 2 (see Fig. 8). In this case the difference is somewhat greater, being around $`0.02`$ rather than $`10^4`$, but the values still appear to be converging for large ancilla dimension. It is expected that the results shown are a good approximation of the true value of $`\mathrm{\Delta }\chi _{U_2}`$.

The results for the third operation, $`U_3`$, are shown in Fig. 9. All results here are for ensembles with four states. In this case it was found that, if the ancilla had dimension 2, the values of $`\mathrm{\Delta }\chi _{U_3}^{(4,2)}`$ and $`\mathrm{\Delta }E_{U_3}^{(2)}`$ were identical. In other respects the results were similar to those for the operation $`U_2`$. The value of $`\mathrm{\Delta }E_{U_3}^{()}`$ was not increased by increasing the ancilla dimension above 2, as for the operations $`U_1`$ and $`U_2`$. The value of $`\mathrm{\Delta }\chi _{U_3}^{()}`$ was increased for larger ancilla dimensions, so for these larger ancilla dimensions $`\mathrm{\Delta }\chi _{U_3}^{()}`$ was not equal to $`\mathrm{\Delta }E_{U_3}^{()}`$.

The difference between $`\mathrm{\Delta }\chi _{U_3}^{()}`$ for higher ancilla dimensions and $`\mathrm{\Delta }\chi _{U_3}^{(4,2)}`$ is shown in Fig. 10. In this case the difference between $`\mathrm{\Delta }\chi _{U_3}^{(4,4)}`$ and $`\mathrm{\Delta }\chi _{U_3}^{(4,2)}`$ is almost $`0.05`$, which is larger than for both $`U_1`$ and $`U_2`$. Nevertheless, the maximum difference between $`\mathrm{\Delta }\chi _{U_3}^{(4,4)}`$ and $`\mathrm{\Delta }E_{U_3}^{(2)}`$ was still comparable with the results for $`U_1`$ and $`U_2`$.

To summarize, our results strongly indicate that $`\mathrm{\Delta }\chi _U`$ is strictly greater than $`\mathrm{\Delta }E_U`$ for most two-qubit operations, rather than simply greater than or equal to, we stated in Refs. berry1 ; berry2 . In addition, our results show that the maximum change in Holevo information is not obtained for ancilla dimensions of two, as appears to be the case for the entanglement. The indications are that there is no ancilla dimension that is sufficiently large that $`\mathrm{\Delta }\chi _U`$ is achieved, so $`\mathrm{\Delta }\chi _U`$ would only be obtained asymptotically in the limit of large ancilla dimension.

## V Conclusions

Our results demonstrate that, for two-qubit operations, there are close relationships between the capacities based on the Holevo information and the capacities based on the entanglement. In most cases, the largest values calculated for $`\chi _U`$ (with the largest ancilla dimension and ensemble size) were equal to $`E_U`$. In some cases values of $`\chi _U`$ were calculated that are slightly above $`E_U`$, but the difference was very small. In addition, the maximum values calculated for $`\mathrm{\Delta }\chi _U`$ were very close to those for $`\mathrm{\Delta }E_U`$. These results indicate that the communication capacities $`\chi _U`$ and $`\mathrm{\Delta }\chi _U`$, which are very difficult to calculate, may be estimated from the entanglement capacities $`E_U`$ and $`\mathrm{\Delta }E_U`$, respectively, which may be calculated far more rapidly.

We have also found that, in all cases tested (except the trivial cases where $`\alpha =0`$ or $`\pi /4`$), $`\mathrm{\Delta }\chi _U`$ is greater than $`\mathrm{\Delta }E_U`$, rather than just greater than or equal to, as stated in Refs. berry1 ; berry2 . As $`\mathrm{\Delta }\chi _U`$ is equal to the asymptotic communication capacity in a single direction bhls , and the asymptotic bidirectional communication capacity is no larger than $`2\mathrm{\Delta }E_U`$ for two qubit operations berry2 , our results show that there is an inevitable trade-off involved in performing bidirectional communication. That is, it is not possible to perform an average communication of $`\mathrm{\Delta }\chi _U`$ in both directions at the same time.

Our results also demonstrate that larger capacities $`\chi _U^{()}`$ and $`\mathrm{\Delta }\chi _U^{()}`$ are obtained when the dimension of the ancillas is increased above 2, in contrast to the case of the entanglement capacities, for which an ancilla dimension of 2 is sufficient. This has been shown analytically for $`E_U`$ in Ref. nielsen , and found numerically for $`\mathrm{\Delta }E_U`$ in Ref. leifer  (our numerical results also support this).

It must be emphasized that there are limitations on the conclusions that can be drawn from the results due to their numerical nature. The capacities found here were found by numerical maximization, which has the drawback that it is possible that, in some cases, a local maximum may have been found, rather than the global maximum. However, the calculations presented here were performed very carefully, and in some cases repeated many times with different random numbers, in order to avoid local maxima. Therefore we are reasonably confident that the results presented here are very close to the global maxima.

A more serious issue is that it is not possible to perform calculations for arbitrarily large ancilla dimensions. Although it has been proven that an ancilla dimension of 2 is sufficient for $`E_U`$, this has not been proven for $`\mathrm{\Delta }E_U`$. Here we have found that $`\mathrm{\Delta }E_U^{()}`$ is unchanged as the ancilla dimension is raised from 2 to 5. This result strongly indicates that the values obtained for $`\mathrm{\Delta }E_U^{()}`$ are the correct asymptotic values. Nevertheless, it is, in principle, possible that larger values of $`\mathrm{\Delta }E_U`$ may be achieved for larger ancilla dimensions than those tested. In order to show conclusively that there is a strict inequality between $`\mathrm{\Delta }\chi _U`$ and $`\mathrm{\Delta }E_U`$, it would be necessary to prove that the asymptotic entangling capacity is achieved for an ancilla dimension of 2.

In the cases of $`\chi _U`$ and $`\mathrm{\Delta }\chi _U`$, it is more difficult to estimate the capacity because the capacity increases with ancilla dimension. For $`\mathrm{\Delta }\chi _U`$ the capacity only increases by a small amount as the ancilla dimension is increased from 2 to 5, indicating that the results are a good approximation of the capacity for arbitrary ancilla dimension. For $`\chi _{U_2}`$ and $`\chi _{U_3}`$, there are significant increases in the capacity with the ancilla dimension up to an ancilla dimension of 4, but there is no increase as the ancilla dimension is increased to 5. This result indicates that an ancilla dimension of 4 is sufficient, though further results for higher ancilla dimensions would be required to be confident of this.

In addition there is the problem that it is not possible to perform calculations for arbitrarily large ensemble sizes. For the operation $`U_1`$ it was found that the capacity was not increased as the ensemble size was increased above 2, so it is likely that an ensemble size of 2 is sufficient for this operation. For the operations $`U_2`$ and $`U_3`$ it was found that the capacities $`\chi _U^{()}`$ and $`\mathrm{\Delta }\chi _U^{()}`$ increased with the ensemble size up to an ensemble size of 4. It is possible that the capacities may be larger for larger ensemble sizes; this is a topic for future research.

*

## Appendix A Numerical Methods

The numerical techniques used to search for the maxima are similar to simulated annealing simulated . For the case of the entanglement capacity $`E_U^{()}`$, a vector of complex numbers, $`(\psi _i)^{(dim_A+dim_B)}`$, represents the initial tensor product state. The first $`dim_A`$ numbers are the coefficients representing Alice’s local state, $`|\varphi _A`$, and the last $`dim_B`$ numbers are the coefficients representing Bob’s local state, $`|\chi _B`$. The total state is simply $`|\mathrm{\Psi }=|\varphi _A|\chi _B`$. The initial values of $`\psi _i`$ were selected using a Gaussian distribution followed by normalization. Throughout this section, all Gaussian distributions for complex variables are real Gaussians multiplied by random phases (with a uniform distribution).

At each step another vector of complex random numbers, $`(\mathrm{\Delta }\psi _i)`$, was selected via a Gaussian distribution. If the state represented by $`(\psi _i+\mathrm{\Delta }\psi _i)/𝒩`$ (where $`𝒩`$ is a normalization constant) gave a larger final entanglement after application of $`U`$, then the coefficients $`(\psi _i)`$ were replaced with $`(\psi _i+\mathrm{\Delta }\psi _i)/𝒩`$. This technique is equivalent to applying simulated annealing with a temperature of zero, because in no case were coefficients chosen that gave a smaller final entanglement. It was found that this technique provided very rapid convergence, and using non-zero temperatures did not provide better convergence.

Initially the standard deviation in the Gaussian distribution used for the increments, $`\sigma `$, was chosen to be 1. The value of $`\sigma `$ was halved each time there were 1000 consecutive increments tested with none providing a larger final entanglement. This process was terminated when $`\sigma `$ fell below $`10^9`$. At this stage the progressive changes in the final entanglement were on the order of $`10^{15}`$ or less.

In the case of the capacity $`\chi _U^{()}`$, a vector of complex random numbers, $`(\psi _i)^{dim(_A_B)}`$, was used to represent one initial state in the ensemble, $`|\psi _{AB}`$. Initial values of $`\psi _i`$ were selected using a Gaussian distribution then normalizing. The remainder of the initial states in the ensemble were obtained by local unitary operations, $`V^{(k)}`$. These unitary operations were represented by matrices of complex numbers, $`V_{ij}^{(k)}`$. The initial values of $`V_{ij}^{(k)}`$ were selected by generating complex random numbers with a Gaussian distribution, then using Gram-Schmidt orthormalization on the row vectors. The probabilities were represented by a vector of real numbers $`(p_k)`$. Initial values were selected using a uniform distribution then normalizing.

At each step, complex random numbers for the increments $`\mathrm{\Delta }\psi _i`$ and $`\mathrm{\Delta }V_{ij}^{(k)}`$ were chosen using a Gaussian distribution with standard deviation $`\sigma `$, and real random numbers for the increments $`\mathrm{\Delta }p_k`$ were chosen using a uniform distribution from $`\sigma /2`$ to $`\sigma /2`$. These increments were used to create a new ensemble with state $`(\psi _i+\mathrm{\Delta }\psi _i)/𝒩`$ and probabilities $`p_k(1+\mathrm{\Delta }p_k)/(1+_kp_k\mathrm{\Delta }p_k)`$. The new operations for the new ensemble were obtained by adding the increments $`\mathrm{\Delta }V_{ij}^{(k)}`$ to $`V_{ij}^{(k)}`$ then applying Gram-Schmidt orthonormalization to the row vectors.

In this case the temperature of the simulated annealing was not taken to be zero. The new ensemble was selected if the new value for the Holveo information after application of $`U`$, $`\chi _{\mathrm{new}}`$, was larger than the previous value, $`\chi _{\mathrm{old}}`$. If the new value was lower, the new ensemble was selected if

$$R<\mathrm{e}^{(\chi _{\mathrm{new}}\chi _{\mathrm{old}})/\tau },$$
(14)
where $`R`$ is a uniform random number between zero and one, and $`\tau `$ is a tolerance equivalent to the temperature. The value of the tolerance $`\tau `$ was initially taken to be $`10^6`$. The value of the final Holevo information was checked every 10000 iterations; if it had decreased, then the tolerance was divided by 2, and if it had increased, then the tolerance was multiplied by $`1.1`$. This adjustment was found to provide reasonably rapid convergence to the final value.

There were two alternative schemes used to adjust the standard deviation, $`\sigma `$, in the Gaussian distribution used for the increments to the ensemble. The first was similar to that for the entanglement, except that $`\sigma `$ was halved if 10000 alternatives were tested with no increase in the Holevo information. The other scheme was to adjust $`\sigma `$ such that approximately 20% of the alternatives were accepted. Only the first scheme was used in the initial part of the calculation (for the first $`10^6`$ or so increments tested). For later parts of the calculation both alternatives were tried.

The numerical techniques used to calculate $`\mathrm{\Delta }E_U^{()}`$ and $`\mathrm{\Delta }\chi _U^{()}`$ were similar to those used for $`E_U^{()}`$ and $`\chi _U^{()}`$, though there are minor differences. In the case of $`\mathrm{\Delta }E_U^{()}`$, the initial state, $`|\psi _{AB}`$ was taken to be a general entangled state, rather than a tensor product of two local states. This state was represented by a vector of coefficients $`(\psi _i)^{dim(_A_B)}`$.

The case of $`\mathrm{\Delta }\chi _U^{()}`$ is rather simpler than that for $`\chi _U^{()}`$. Rather than it being necessary to condider a set of unitary operations and probabilities, the ensemble was simply represented by a set of complex coefficients $`\psi _i^{(k)}`$ for all of the states in the ensemble. Rather than separately storing probabilities, these states were not normalized, and the normalizations of these states were taken to be the probabilities.

###### Acknowledgements.

This research has been supported by an Australian Department of Education Science and Training Innovation Access Program Grant to support collaboration in the European Fifth Framework project QUPRODIS, and by Alberta’s informatics Centre of Research Excellence (iCORE).