text:
# Cosmic Rays and Neutrinos from GRBs: Predictions versus Acceleration Modeling

## 1 Introduction

Gamma Ray Bursts (hereafter GRBs) are unique high energy phenomena in astrophysics because of their possibility to manifest all the interesting ”astroparticles” processes, such as generation of high energy gamma rays, ultra high energy cosmic rays, high energy neutrinos and gravitational waves (see Dermer (2001) and Mészáros (2002)). The ”fireball” model (Rees & Mészáros, 1992) has been successful in explaining the afterglow stage, and its more elaborated form with the addition of internal shocks (Rees & Mészáros (1994), Paczyński & Xu (1994)) during the ”free” expansion has been successful in explaining the light curve and gamma spectra (Fishman & Meegan (1995), Beloborodov et al. (2000)). Because of the highly relativistic dynamics necessarly involved in the GRB phenomenon, the generation of very high energy particles is expected through strong shocks and strong magnetic perturbations. Indeed GRB population could be a main source of UHE-cosmic rays which could be generated by the external shocks (Vietri, 1995) or by the internal shocks (Waxman, 1995), and a flux of neutrinos produced by the collisions of the UHE-cosmic rays with the GRB gamma photon is reasonably expected (Waxman & Bahcall, 1999).

In this paper, our intend is to look at the sensitivity of the predictions of cosmic ray and neutrino generation to the description of the Fermi acceleration process together with opacity effects for protons. We have analysed the consequences of two assumptions, one is the so called ”Bohm scaling” assumption, the other is what we called ”Kolmogorov scaling”. The former consists in the statement that the Fermi acceleration time is proportional to the Larmor time of the accelerated particle, with a constant proportionality factor larger than unity (often chosen between 1 and 10 in the litterature). The latter consists in taking into account that the ratio between these two characteristic times depends on the rigidity of the particle through a law governed by the turbulence spectrum, as is confirmed by numerical works (Casse et al., 2001) , where no Bohm scaling has been found. In astrophysical media such as the solar wind, the interstellar medium, the turbulence spectrum is likely consistent with the Kolmogorov law.

The Bohm scaling is very convenient to make a first investigation of the high energy physics performances of the objects, and our prejudices are often grounded on this. However, in this paper, we will show that this scaling leads to unrealistic results. The analysis based on the Kolmogorov scaling is more reliable, and it turns out that it leads to very different conclusions, as we will present with details in the paper. The analysis is developed both in the radiation free stage and in the radiative stage. We especially emphasise the regime when the fireball is opaque to pp-collisions during the beginning of the internal shocks stage because of its baryonic load, which should occur for a large fraction of GRBs population. We make unusual predictions about the performance of cosmic ray generation, and about neutrino emission that could have two components, namely a non-thermal pp-neutrino emission and a p$`\gamma `$-neutrino emission.

The paper is organised as follows. In section 2, we briefly present the dynamical description we need for the estimation of the opacity effects, not only for photons but also for protons. In section 3, we analyse the consequence of the fast acceleration regime governed by the Bohm scaling. In section 4, we present the properties of the more progressive acceleration process governed by the Kolmogorov scaling and predict a reasonable non-thermal pp-neutrino emission with its spectrum. We investigate, more briefly in section 5, the radiation free stage, where UHE-cosmic rays are expected to be accelerated and p$`\gamma `$-neutrinos generated. We end the paper with a discussion that summaries our conclusions about cosmic ray and neutrino generations and the sensitivity of the predictions to the acceleration model.

## 2 Preliminary considerations

### 2.1 Dynamics of the fireball

In this subsection, we summarize all the results we need for this paper that describe the expansion of the fireball (Mészáros et al., 1993). The wind flow is considered to be a set of discrete shells which are successively emitted with an energy $`E_s=E/N_s`$, where $`N_s`$ is the total number of shells. The duration, $`t_w`$, of this wind flow provides with an interval of shell number, namely, $`1N_sct_w/r_0`$ where $`r_0`$ is the size of the central object.

At the very beginning of the expansion of a shell, the pressure is supposed to be dominated by the radiative pressure. The temperature, $`T`$, of the plasma which is mainly composed by electron-positron pairs is equal to the photon temperature. Considering a shell is initially spherical with a radius, $`r_0`$, and with an energy emitted in $`\gamma `$-rays equal to $`E_s`$, we have

$`{\displaystyle \frac{E_s}{(4/3)\pi r_0^3}}=aT^4,`$
where $`a`$ is the Stefan constant. Thus, the plasma temperature is

$`T=\left({\displaystyle \frac{3E_s}{4\pi ar_0^3}}\right)^{1/4}.`$
For an energy, $`E`$, in the neighborhood of $`10^{51}`$ ergs, a number of shells, $`N_s`$, of 20 and a radius, $`r_0`$, of $`10^7`$ cm, the temperature is below 10 MeV.

When a shell starts, the energy $`E_s`$ is very upper to the baryon mass energy. We can define the ratio, $`\eta `$, between these two energies which is

$`\eta ={\displaystyle \frac{E_s}{(M_b/N_s)c^2}}={\displaystyle \frac{E}{M_bc^2}}1`$
where $`M_b`$ is the total baryonic mass ejected.

In the observer frame, the shell thickness, $`\mathrm{\Delta }r`$, is supposed to remain constant and equal to $`r_0`$ until the broadening radius $`r_b`$ (Goodman (1986), Mészáros et al. (1993)). Beyond this radius, defining the Lorentz factor, $`\mathrm{\Gamma }`$, of the baryonic matter, the thickness becomes

$`\mathrm{\Delta }r(r/c)\mathrm{\Delta }vr/2\mathrm{\Gamma }^2`$ (1)
In the same frame, we also define a radius, $`r_s`$, where the kinetic energy of baryonic matter reaches its saturation value. At this moment, the Lorentz factor $`\mathrm{\Gamma }=\mathrm{\Gamma }_{max}`$ is close to $`\eta `$.

In the co-moving frame, according to the Lorentz transformation, the shell thickness is given by

$`\mathrm{\Delta }R=\mathrm{\Gamma }\mathrm{\Delta }r.`$
If we consider an adiabatic expansion of the shell, conservation of entropy in the co-moving frame for a radius lower to $`r_s`$ is given by

$`r^2\mathrm{\Gamma }r_0T^3=constant.`$ (2)
At the same time, conservation of energy is such that

$`r^2\mathrm{\Gamma }^2r_0T^4=constant.`$ (3)
From (2) and (3), we deduce two laws of evolution for $`rr_s`$ which are $`\mathrm{\Gamma }(r)r`$ and $`T(r)r^1`$. Then, the saturation radius $`r_s`$ can be defined by

$`r_s=\eta r_0.`$
For $`\eta `$ of the order of 300, we obtain $`r_s3\times 10^9`$ cm.

Beyond the radius $`r_s`$, the equation (2) is always valid but the Lorentz factor of the shell remains constant. In this case, the temperature is such that $`T(r)r^{2/3}`$. Considering $`\mathrm{\Delta }rr_0`$, the radius $`r_b`$, according to (1), satisfies $`r_b\mathrm{\Gamma }^2r_0`$. At last, $`\mathrm{\Gamma }`$ is supposed to reach its saturation value, $`\eta `$, around $`r_s`$ well before $`r_b`$. In fact, $`\eta `$ is the average value of the saturation bulk Lorentz factor, and we have to bear in mind that deviations from this average is expected to generate internal shocks. The radius $`r_b`$ is so defined by

$`r_b\eta ^2r_0.`$
For $`\eta `$ of the order of 300, we obtain $`r_b9\times 10^{11}`$ cm.

To conclude, all the different parameters of a shell we need, can be summarized by the next following expressions :

In the observer frame, the thickness of a shell is

$`\mathrm{\Delta }r\{\begin{array}{ccc}& r_0& \mathrm{for}rr_b\hfill \\ & r/\mathrm{\Gamma }_{max}^2r/\eta ^2& \mathrm{for}rr_b\hfill \end{array}.`$
The Lorentz factor of a shell is such that

$`\mathrm{\Gamma }(r)\{\begin{array}{ccc}& r/r_0& \mathrm{for}rr_s\hfill \\ & \mathrm{\Gamma }_{max}\eta & \mathrm{for}rr_s\hfill \end{array}.`$
In the co-moving frame, the thickness of a shell becomes equal to

$`\mathrm{\Delta }R\{\begin{array}{ccc}& \mathrm{\Gamma }r_0r& \mathrm{for}rr_s\hfill \\ & \mathrm{\Gamma }_{max}r_0\eta r_0& \mathrm{for}r_srr_b\hfill \\ & r/\mathrm{\Gamma }_{max}r/\eta & \mathrm{for}rr_b\hfill \end{array}.`$ (9)
The temperature of the pair electron-positron plasma is such that

$`T(r)\{\begin{array}{ccc}& 10\times \left(\frac{T(r_0)}{10MeV}\right)\left(\frac{r_0}{r}\right)MeV& \mathrm{for}rr_s\hfill \\ & 10\times \left(\frac{T(r_0)}{10MeV}\right)\left(\frac{r_0}{r_s}\right)\left(\frac{r_s}{r}\right)^{2/3}MeV& \mathrm{for}rr_s\hfill \end{array},`$
and the co-volume of a shell, $`V_c(r)`$, is

$`V_c(r)\{\begin{array}{ccc}& \mathrm{\Omega }r^3& \mathrm{for}rr_s\hfill \\ & \mathrm{\Omega }\eta r_0r^2& \mathrm{for}r_srr_b\hfill \\ & (\mathrm{\Omega }/\eta )r^3& \mathrm{for}rr_b\hfill \end{array},`$ (14)
where $`\mathrm{\Omega }`$ ($`4\pi /500`$) is the opening angle of the emission.

The collision date (if any…), $`t_c`$, between two shells of Lorentz factor $`\gamma _1`$ and $`\gamma _2`$, which started at two different times separated by $`\mathrm{\Delta }t_0`$ is such that $`t_c2\frac{\gamma _1^2\gamma _2^2}{|\gamma _1^2\gamma _2^2|}\mathrm{\Delta }t_0`$, which leads to collisions at various distances $`r_c`$ such that $`r_c/r_b\mathrm{\Delta }t_0/t_0`$ (see Daigne & Mochkovitch (1998)). For a long GRB having a duration $`t_w10`$ s, the maximum collision radius reaches the deceleration radius $`r_d`$. The light curve duration could be due to either the duration of the wind $`t_w`$ or to the time spread of the shell at the most remote collision, namely, $`\mathrm{\Delta }t_d=t_d/\eta ^2`$. It turns out that this latter time is comparable to the duration of the wind of a long GRB, namely, $`110`$ s. The observed millisecond variations should come from internal collisions located at a few $`r_b`$ where $`\mathrm{\Delta }r_br_0`$.

Anyway after a few $`r_b`$, all the shells mix up and form a single jet. That is the reason why volume and thickness of a shell can really be defined only before the broadening radius. After this radius, in the equations (9) and (14), they are undervalued.

### 2.2 Conversion of the GRB energy into cosmic rays

During the evolution, the energy of the fireball is shared between several forms : the thermal energy $`E_{th}`$, the magnetic energy $`E_m`$, the bulk kinetic energy, the cosmic ray energy $`E_{}`$ and the energy in the form of hydromagnetic perturbations $`E_m^{}`$. The energy share $`E_m^{}`$ is the reservoir for particle acceleration.

Until the saturation radius $`r_s`$, the ratio of the pair thermal energy over the fireball energy, $`(aT^4V_c)/E_s`$, starts with a value close to unity and then decays as $`1/r`$. In this stage, the proton population is a tiny contribution both in number of particles and in energy. Assuming that the protons dominate the baryonic load, we get the total number of protons

$$N_p\frac{M_b}{m_p}0.67\times 10^{51}\left(\frac{\eta }{300}\right)^1\left(\frac{E}{10^{51}erg}\right).$$

Thus, the density of protons by shell, $`n_p(r)`$, in the co-moving frame and before the broadening radius, is such that

$`n_p(r)10^4\left({\displaystyle \frac{t_0}{1ms}}\right)\left({\displaystyle \frac{t_w}{10s}}\right)^1{\displaystyle \frac{N_p}{V_c(r)}},`$
and the maximum value of $`n_p(r)`$ is equal to $`N_p/V_c(r)`$.
For $`t_0=r_0/c=1`$ ms, we can write

$`n_p(r)\{\begin{array}{ccc}& 2.7\times 10^{27}\left(\frac{t_w}{10s}\right)^1\left(\frac{\eta }{300}\right)^1\left(\frac{\mathrm{\Omega }/4\pi }{2\times 10^3}\right)^1\left(\frac{E}{10^{51}erg}\right)\left(\frac{r}{r_0}\right)^3cm^3& \\ & 1.0\times 10^{20}\left(\frac{t_w}{10s}\right)^1\left(\frac{\eta }{300}\right)^2\left(\frac{\mathrm{\Omega }/4\pi }{2\times 10^3}\right)^1\left(\frac{E}{10^{51}erg}\right)\left(\frac{r}{r_s}\right)^2cm^3& \end{array}`$
respectively for $`rr_s`$ and $`r_srr_b`$.
Considering $`\eta =300`$, even if the GRB has a long duration about 10 s, the density of protons, at $`r=r_s`$, is larger than $`10^{20}`$ cm-3. We note that, up to the saturation radius, the density of a shell is comparable to the density of a solid.

The corresponding ratio of the thermal energy of protons, $`3T/(2\eta m_pc^2)`$, starts with a very low value ($`10^4`$) and decays like $`1/r`$ as well. The cosmic ray component is supposed to develope out of the proton thermal component and would grow up to several percents of $`E`$ to account for a significant contribution of the UHE cosmic ray generation in the Universe. Assume that a fraction $`\xi _{}`$ of the proton population is injected into the cosmic ray component (i.e. $`N_{}=\xi _{}N_p`$) and that they reach a mean energy $`\overline{ϵ}(r)=\overline{\gamma }(r)m_pc^2`$ in the co-moving frame. Then, the contribution to the energy of the fireball, in the observer frame, is given by

$$E_{}(r)=\frac{\xi _{}\overline{\gamma }(r)\mathrm{\Gamma }(r)}{\eta }E$$

When the radius $`r`$ reaches the saturation radius $`r_s`$, the energy of the cosmic ray component is equal to $`\xi _{}\overline{\gamma }(r)E/(1+4\xi _{}\overline{\gamma }(r)/3)`$. For a cosmic ray spectrum in $`ϵ^2`$, $`\overline{\gamma }=\mathrm{log}(\gamma _{max})`$, which clearly shows that the goal of converting about $`10`$ percents of the fireball energy into cosmic rays energy is achieved when a sizeable fraction of the protons are injected in the cosmic ray population.

Regarding the magnetic energy of the fireball, two points of view can be considered; either the magnetic field behaves like in a jet and the poloidal component decreases like $`1/r^2`$ thus its pressure decreases like the relativistic thermal pressure in $`1/r^4`$; whereas the toroidal component decreases more slowly; or the magnetic field of the shell disconnects from the central source; which is very likely. In this latter case, the magnetic energy of the shell is conserved as long as the dissipation is negligible (e.g. $`BV_c^{1/2}`$) and decays when the Fermi processes and/or reconnections become efficient in accelerating particles. Initially, the intensity of magnetic field is usually considered as of the order of the equipartition value, which means that the magnetic energy is a sizeable fraction of the fireball energy ($`E_m(r_0)=\xi _mE`$).

### 2.3 The importance of pp-collisions during the primeval stage

We begin with defining the radius $`r_{}`$ where a shell becomes optically thin with respect to Compton scattering. It can easily be checked that a typical shell width $`\mathrm{\Delta }R`$ becomes smaller than the flow transverse radius after a short while, when $`r>\eta \sqrt{\pi /4\mathrm{\Omega }}r_0`$ which is comparable to $`r_s`$. It will turn out that the photosphere is located at a much larger distance for large enough $`\eta `$ (see (18)) and therefore the opacity of a shell is determined by its width. The Compton opacity is $`\tau _{}=\sigma _Tn_e\mathrm{\Delta }R`$ where $`\sigma _T`$ is the Thompson cross section. We assume $`n_en_pn_b`$ and thus the co-moving baryon density is related to the wind mass flux $`\dot{M}_w`$ by

$$n_b=\frac{\dot{M}_w}{\mathrm{\Gamma }\mathrm{\Omega }r^2m_pc}$$

Moreover for $`r>r_s`$, the GRB energy is converted into relativistic kinetic energy such that the kinetic luminosity (kinetic energy flux) of the wind $`L_w=E/t_w=\eta \dot{M}_wc^2`$. For $`r>r_s`$, $`\mathrm{\Gamma }\eta `$ and the Compton opacity is thus

$$\tau _{}=\frac{\sigma _TL_w}{\mathrm{\Omega }m_pc^3\eta ^2}\frac{\mathrm{\Delta }R}{r^2}.$$

We can define a critical value for $`\eta `$ such that the photospheric radius is located at $`r_b`$, where shock acceleration starts. This critical value $`\eta _{}`$ is given by

$$\eta _{}\left(\frac{\sigma _TL_w}{\mathrm{\Omega }m_pc^3r_0}\right)^{1/5}570\times \left(\frac{E}{10^{51}erg}\right)^{1/5}\left(\frac{t_w}{1s}\right)^{1/5}\left(\frac{\mathrm{\Omega }/4\pi }{2\times 10^3}\right)^{1/5}$$

The photospheric radius is such that

$`r_{}=\{\begin{array}{ccc}r_b\left(\eta _{}/\eta \right)^5& \mathrm{for}\eta \eta _{}& \\ r_b\left(\eta _{}/\eta \right)^{5/2}& \mathrm{for}\eta \eta _{}& \end{array}.`$ (18)
These simple formulae illustrate the requirement of a large value of the asymptotic bulk Lorentz factor $`\eta `$ in order to observe an optically thin X-ray spectrum. For $`\eta `$ larger than the critical value, the internal shocks generate only a non-thermal spectrum just after $`r_b`$, whereas, for $`\eta <\eta _{}`$, the internal shocks start accelerating particles in an optically thick plasma. Assuming that the GRB energy varies within two orders of magnitude, the wind time by two orders of magnitude and the solid angle by one order of magnitude, the possible values of $`\eta _{}`$ extend within an interval of one order of magnitude. However it is worth mentioning that this usual approach of the opacity issue is an underestimate of the opacity effects, because the coalescence of shells and ultimately their dissolution into a single jet increase the opacity.

The possibility for the internal shock scenario to partially operate in an optically thick regime is important when one consider pp-collisions. Indeed, any relativistic protons can collide with other protons (or neutrons) producing pions and thus neutrinos. The cross section, $`\sigma _{pp}`$, for this production is constant and equal to $`2.7\times 10^{26}`$ cm2 when the kinetic energy of protons is larger than 1 GeV. The sheets become thin to pp-collisions always before the photosphere since $`\tau _{pp}=n_p\sigma _{pp}\mathrm{\Delta }R=1`$ at $`r_{pp}`$ such that $`r_{pp}=r_{}\sqrt{\sigma _{pp}/\sigma _T}<r_{}`$ for $`r_{pp}<r_b`$ and $`r_{pp}=r_{}\sigma _{pp}/\sigma _T<r_{}`$ for $`r_{pp}>r_b`$. In this latter case, $`r_{pp}0.04r_{}`$; which occurs when $`\eta \eta _{}/2`$. In this paper, we intend to emphasize the importance of this energy limitation to proton acceleration and estimate the resulting neutrino emission.

## 3 Consequences of a fast acceleration regime

In this section, we analyse the consequences of the usual assumption of a Bohm scaling for the Fermi process for both the energy distribution cut off and the depletion of the energy reservoir for acceleration.

### 3.1 Bohm scaling and energy losses

We first consider a Fermi acceleration process that is assumed to follow the so called Bohm’s scaling namely characterized by its time scale proportional to the Larmor time, i.e. :

$$t_{acc}=\kappa _0t_L=\kappa _0\frac{10^4}{Z}\left(\frac{ϵ}{1GeV}\right)\left(\frac{B}{1G}\right)^1s,$$

where we take $`\kappa _0=10`$, as often assumed in astroparticle physics.

In the early phase, the cut off of the proton energy distribution can be caused by either the synchrotron loss or the pp-collisions; which is obtained by equating the acceleration time and the loss time. The synchrotron time $`t_{syn}`$ is given by

$$t_{syn}^1=\frac{4Z^4\sigma _T\gamma }{3mc}\left(\frac{m_e}{m}\right)^2W_m,$$

where $`W_m`$ is the density of magnetic energy.

This expression leads to the cut off energy :

$$ϵ_{syn}2.4\times 10^{11}\frac{1}{\sqrt{\kappa _0}}\left(\frac{B}{1G}\right)^{1/2}GeV.$$

We consider now the limitation caused by pp-collisions. The cut off energy $`ϵ_{pp}`$ is such that $`(n_p\sigma _{pp}c)\kappa _0t_L=1`$, and thus

$$ϵ_{pp}\frac{10^4}{\kappa _0}\left(\frac{B}{1G}\right)\left(\frac{\eta ^6}{\eta _{}^5}\right)\left(\frac{r_0}{c}\right)GeV.$$

At the broadening radius ($`r_b`$), where the magnetic field could be as high as $`10^7`$ G, these two limitations are comparable :

$$ϵ_{syn}2.4\times 10^7\left(\frac{\kappa _0}{10}\right)^{1/2}\left(\frac{B}{10^7G}\right)^{1/2}GeV,$$

and

$$ϵ_{pp}3\times 10^9\left(\frac{\kappa _0}{10}\right)^1\left(\frac{B}{10^7G}\right)\left(\frac{\eta }{300}\right)\left(\frac{\eta }{\eta _{}}\right)^5GeV$$

During the expansion the ratio between these two limitations evolves such that

$$\frac{ϵ_{pp}}{ϵ_{syn}}1.2\times 10^3\left(\frac{\kappa _0}{10}\right)^{1/2}\left(\frac{\eta }{300}\right)\left(\frac{\eta }{\eta _{}}\right)^5\left(\frac{B}{10^7G}\right)^{3/2}\left(\frac{r}{r_b}\right)^2$$

which suggests that for any decrease of the magnetic field in $`r^\alpha `$ with $`\alpha `$ between 1 and 2, the most severe limitation would be due to synchrotron losses. However we will show in the next subsection that such efficient Fermi acceleration would dissipate the magnetic energy very rapidly.

### 3.2 The depletion of the energy reservoir for acceleration

Such an efficient acceleration produces a strong depletion of the energy reservoir for particle acceleration. The previous results would make sense only if the depletion of the energy reservoir for particle acceleration is slow. This is exactly the purpose of this subsection.

Let $`E_m^{}`$ be the energy in the hydromagnetic perturbations (magnetic and kinetic energy) involved in the acceleration of particles. We assume that this is a sizeable fraction $`\xi _m^{}`$ of the total magnetic energy and that even the large scale magnetic field in the shell is mostly disconnected from the central source after a short while and thus is rather tangle. Therefore we assume that the total magnetic energy dissipates at the same rate than $`E_m^{}`$. A minimum dissipation is assumed by considering its energy loss by particle acceleration only. According to (3.1), the acceleration power is $`Q^+=\frac{N_{}}{\kappa _0t_{L0}^{}}m_pc^2`$, where $`t_{L0}^{}`$ is the Larmor time for a proton of 1 GeV in the co-moving frame. The minimum depletion is then governed by $`\dot{E}_m^{}=Q^+`$, which reads :

$$\dot{E}_m=\frac{E}{t_{}^{}}\frac{\overline{B}}{B_{}},$$
(19)
where we define a characteristic time $`t_{}^{}`$ at a given radius $`\stackrel{~}{r}`$ where the magnetic field has an intensity $`B_{}`$, such that

$$t_{}^{}(\eta \kappa _0\xi _m^{}/\xi _{})t_{L0}^{}(B_{}).$$

Clearly, this time measures the rate of depletion at the considered radius $`\stackrel{~}{r}`$ and must be compared to the co-moving dynamical time $`t^{}(\stackrel{~}{r})=\stackrel{~}{r}/\mathrm{\Gamma }c`$. Let us make an estimate at $`r_b`$, the late stage actually; for $`B=10^7`$ G and $`ϵ=1`$ GeV, $`t_{}^{}3\times 10^8`$ s, whereas $`t^{}(r_b)10^2`$ s. This perturbation burning out within a few nanoseconds is not realistic at all. Before this stage, the depletion time is much shorter than the fireball dynamical time by a factor proportional to $`(Br)^1`$ that goes at best like $`r^{1/2}`$ for $`Br^{3/2}`$ corresponding to magnetic energy conservation and the perturbations are burnt out almost immediately, their energy being radiated by photons and neutrinos. The observation of these perturbations in the light curve of gamma emissions suggests a much slower depletion; which we will analyse in the next section.

## 4 The interest of a progressive acceleration regime

The Bohm scaling is often used in astroparticle physics, which sometimes can provide some rough estimate of the high energy cut off. But the true scaling (Casse et al., 2001) depends on the turbulence spectrum. We will assume that the perturbations are distributed according to the Kolmogorov law. This will significantly change the performance of the Fermi process and the conclusions about cosmic ray generation.

### 4.1 Kolmogorov scaling and energy losses

In the co-moving frame, the acceleration time depends on the speed $`\beta _{}c`$ of magnetic perturbations that scatter particles on both sides of the shocks. Considering the same average velocity, we have the relation

$`t_{acc}={\displaystyle \frac{t_s}{\beta _{}^2}},`$
where the scattering time, $`t_s`$, can be expressed like

$`t_s=(\eta _b\rho ^{\beta 1}\omega _L)^1`$
with $`\eta _b=\frac{<\delta B^2>}{\overline{B}^2}`$ and $`\rho =\frac{r_L}{l_c}1`$ where $`ł_c`$ is the mean coherence length of the magnetic perturbations whose spectrum is supposed to be a power law of index $`\beta `$. According to the Kolmogorov theory, one can take $`\beta =5/3`$ and, with $`\kappa _0=1/(\beta _{}^2\eta _b)`$, we get

$`\kappa =\kappa _0\rho ^{2/3}.`$
Considering $`l_c\mathrm{\Delta }R`$,

$`\rho ={\displaystyle \frac{r_L}{\mathrm{\Delta }R}}{\displaystyle \frac{ϵ}{ZeB\mathrm{\Delta }Rc}}`$
and we have

$`\kappa 2.2\times 10^3\left({\displaystyle \frac{\kappa _0}{10}}\right)\left({\displaystyle \frac{ϵ}{1GeV}}\right)^{2/3}\left({\displaystyle \frac{B}{1G}}\right)^{2/3}\left({\displaystyle \frac{\mathrm{\Delta }R}{1cm}}\right)^{2/3}.`$
The synchrotron limitation becomes

$`ϵ_{syn}1.2\times 10^{19}\left({\displaystyle \frac{\kappa _0}{10}}\right)^{3/4}\left({\displaystyle \frac{B}{1G}}\right)^{5/4}\left({\displaystyle \frac{\mathrm{\Delta }R}{1cm}}\right)^{1/2}GeV.`$ (20)
For $`r=r_b`$ ($`10^{12}`$ cm), $`B=10^7`$ G and $`\mathrm{\Omega }=4\pi /500`$, $`ϵ_{syn}3.8\times 10^5`$ GeV.

As for the pp-collision limitation,

$`ϵ_{pp}\left({\displaystyle \frac{\kappa _0}{10}}\right)^3\left({\displaystyle \frac{\eta }{300}}\right)^3\left({\displaystyle \frac{\eta }{\eta _{}}}\right)^{15}\left({\displaystyle \frac{B}{10^7G}}\right)\left({\displaystyle \frac{\mathrm{\Delta }R}{10^{12}cm}}\right)^2\left({\displaystyle \frac{r}{r_b}}\right)^6GeV.`$
For $`rr_b`$, $`\mathrm{\Delta }R=r/\eta `$, defining $`r_b\eta _{}^2r_0`$, we obtain

$$ϵ_{pp}2.0\times 10^5\left(\frac{\kappa _0}{10}\right)^3\left(\frac{B}{10^7G}\right)\left(\frac{\eta _{}}{570}\right)\left(\frac{\eta }{\eta _{}}\right)^8\left(\frac{r}{r_b}\right)^4GeV.$$
(21)
This energy increases up to a maximum value reached at $`r_{pp}`$ (when it’s larger than $`r_b`$), namely,

$$ϵ_{pp}3.5\times 10^2\left(\frac{\kappa _0}{10}\right)^3\left(\frac{B}{10^7G}\right)\left(\frac{\eta _{}}{570}\right)\left(\frac{\eta _{}}{\eta }\right)^4GeV.$$
(22)
As long as $`r<r_{pp}`$, the energy limitation is clearly due to pp-collisions rather than synchrotron losses. These results show, first, that it is not possible to accelerate protons beyond 100 GeV before the fireball becomes thin for protons at $`r_{pp}`$, for $`\eta >\eta _{}`$; second, that, for $`\eta <\eta _{}/2`$, the proton energy could increase above 100 GeV. For instance, if $`\eta =\eta _{}/3`$, $`r_{}240r_b`$, $`r_{pp}10r_b`$, $`B`$ weakens by a factor $`(r_{pp}/r_b)^\alpha 3.2\times 10^2`$ with $`\alpha =3/2`$ and thus compensates the increase due to $`\eta `$. Therefore the energy remains of order 100 GeV. We have to examine how this regime changes the rate of perturbation energy depletion.
In the opposite situation where the protons are not accelerated beyond GeV in the opaque stage, for $`\eta >\eta _{}`$ and/or $`B<10^5`$ G at $`r_b`$, nucleosynthesis is possible as shown by Lemoine (Lemoine, 2002).

### 4.2 The depletion of the energy reservoir for acceleration

Taking into account the variation of $`\kappa `$ in the acceleration power, we define a new characteristic time $`t_{}^{}`$ which is quite different from the last one :

$$t_{}^{}(\eta \kappa _0<\rho ^{2/3}>^1\xi _m^{}/\xi _{})t_{L0}^{}(B_{}).$$

For a power law distribution, not harder than $`ϵ^2`$, $`<\rho ^{2/3}>\rho _0^{2/3}`$, $`\rho _0`$ being the rigidity for $`1`$ GeV. For bell type distribution of standard deviation $`\overline{ϵ}`$, $`<\rho ^{2/3}><\rho (\overline{ϵ})^{2/3}>`$. Compared to the Bohm regime, the Kolmogorov regime increases this time by a factor $`\rho _0^{2/3}`$. At $`r_b`$ we get a more extended time of the order of $`0.1`$ s, which is much more realistic than in the Bohm scaling case. If $`r_{pp}`$ is significantly larger, the depletion behaviour is changed such that $`t_{}^{}/t^{}`$ decreases instead of increasing like in Bohm regime. Indeed in Kolmogorov regime, this ratio scales like $`B^{1/3}r^{5/9}`$, which decreases when $`B`$ decreases less slowly than $`r^{5/3}`$.

For a more detailed estimate, we integrated the differential equation (19) which, after some algebra, leads to the following result :

$$\frac{E_m}{E_m(r_b)}=\left[1\frac{E}{E_m(r_b)}\frac{t^{}(r_b)}{5t_{}^{}(r_b)}\left(1\left(\frac{r}{r_b}\right)^{5/6}\right)\right]^6.$$

In fact, the perturbations are burnt out in a fraction $`\delta t^{}`$ of $`t_{}^{}(r_b)`$, which can be obtained directly from the differential equation (19) :

$$\delta t^{}=t_{}^{}(r_b)\frac{E_m(r_b)}{E},$$

which confirms the previous statement of a reasonable depletion time.

### 4.3 The pp-neutrino emission

In this subsection, we give an evaluation of the number of emitted pp-neutrinos and we calculate their energy spectrum.

The pp-collisions produce neutrinos after some reactions which are

$`p+p`$ $`D+\pi ^+`$
$`p+p+a(\pi ^++\pi ^{})+b\pi ^0`$
$`p+n+\pi ^++a(\pi ^++\pi ^{})+b\pi ^0`$
$`2n+2\pi ^++a(\pi ^++\pi ^{})+b\pi ^0`$
The production of mesons $`\pi ^{}`$ and $`\pi ^+`$ gives neutrinos et $`\mu `$-mesons through the decay reactions

$`\pi ^{}\mu ^{}+\overline{\nu }_\mu `$
$`\pi ^+\mu ^++\nu _\mu `$ (23)
and other neutrinos are produced after

$`\mu ^{}e^{}+\overline{\nu }_e+\nu _\mu `$
$`\mu ^+e^++\nu _e+\overline{\nu }_\mu .`$ (24)
Let’s write energy of different particles

$`m_{\pi ^\pm }`$ $``$ $`140MeV/c^2`$
$`m_\mu `$ $``$ $`105MeV/c^2.`$
$`m_{e^\pm }`$ $``$ $`0.5MeV/c^2`$
Thus, neutrinos coming from (23) have a minimum energy equal to

$`ϵ_{\nu _\mu }{\displaystyle \frac{m_{\pi ^\pm }^2m_\mu ^2}{2m_{\pi ^\pm }}}30MeV`$
and decay of $`\mu `$-mesons (24) gives neutrinos having a minimum energy which varies from 25 to 50 MeV.
Let $`f_\nu (ϵ_\nu )`$ be the energy distribution of neutrinos, normalised such that $`f_\nu (ϵ_\nu )𝑑ϵ_\nu =n_\nu `$, the number density of neutrinos; and $`f_{}(\gamma )`$ the Lorentz factor distribution of relativistic protons, normalised such that $`f_{}(\gamma )𝑑\gamma =n_{}`$, the number density of relativistic protons. The kinetic equation for neutrinos can be written in a simplified way as follows :

$$\frac{}{t}f_\nu +c\stackrel{}{n}.f_\nu =\xi _\nu \nu _{pp}f_{}(\gamma )\delta (ϵ_\nu \gamma _cϵ_0)𝑑\gamma ,$$
(25)
where $`ϵ_0`$ is the average energy of the neutrinos generated by the pp-collision in the center of mass frame (above estimated), $`\xi _\nu `$ the average number of produced neutrinos at each collision and $`\gamma _c`$ the Lorentz factor of the collision frame. The energy of the emitted neutrinos remain close to the minimum values because the pions do not take all the energy available above the threshold, but just a little excess above it, the remaining energy being kept by the proton or the neutron. Therefore the number $`\xi _\nu `$ is just a few. The Lorentz factor $`\gamma _c\sqrt{\gamma }`$, where $`\gamma `$ is the Lorentz factor of the relativistic proton. So we can easily deduce that, for $`f_{}\gamma ^s`$, the energy flux of neutrinos is in $`ϵ_\nu ^{2s+2}`$, between $`ϵ_0`$ and $`ϵ_0\sqrt{ϵ_{pp}/m_pc^2}`$, in the co-moving frame. For $`s=2`$, a neutrino spectrum in $`ϵ_\nu ^2`$ is expected between 5 GeV and 150 GeV, typically, for the observer frame seing the GRB flow coming front.

Assuming an isotropic neutrino emission in the co-moving frame, the power emitted is

$$\dot{E}_\nu ^{}=\mathrm{\Omega }\xi _\nu ϵ_0\nu _{pp}n_{}r^2d(\mathrm{\Delta }R),$$
(26)
where the integral over the width extends down to the lower bound at $`r_{pp}`$, corresponding to $`\tau _{pp}=1`$. Thus $`E_\nu ^{}=t_w^{}\mathrm{\Omega }\xi _\nu ϵ_0cn_{}r_{pp}^2`$. The quantity $`t_w^{}\mathrm{\Omega }cn_{}r_{pp}^2`$ is the number of relativistic protons $`\xi _{}N_p`$ which have flown through the ”proto”-sphere during the GRB event. Finally, we get a simple formula giving the amount of neutrino energy in the co-moving frame :

$$E_\nu ^{}=\xi _{}\xi _\nu \frac{ϵ_0}{m_pc^2}\frac{E}{\eta }.$$
(27)
For an observer seing the GRB shells coming front, the neutrino energy emission is multiplied by the bulk Lorentz factor, thus $`E_\nu =\xi _{}\xi _\nu (ϵ_0/m_pc^2)E`$ which reasonably leads to $`E_\nu 10^310^2E`$. The number of emitted neutrinos by a GRB is simply $`\xi _\nu N_{}`$. A neutrino telescope of collecting surface $`A`$ can receive $`N_\nu A/4\pi D^2`$ from a GRB exploding at a distance $`D`$. For example, a GRB occuring at 100 Mpc can provide with $`10^5`$ neutrinos crossing a $`km^2`$-detector. However, the number of events recorded by the detector is obtained by multiplying this number by the detection probability which is quite low.

In the same process, $`\pi ^0`$-decay occurs and generates gamma photons. An excess of these photons could be seen in the gamma spectrum; they superimpose a spectrum in $`ϵ_\gamma ^{2s+2}`$ on the synchrotron spectrum due to the electrons. Because we focused on proton acceleration and neutrino emission, we disregarded the interesting issue of the neutron component which can decouple under some conditions (Derishev et al., 1999).

## 5 The stage of UHE Cosmic Rays acceleration

Because of their magnetic field intensity and the size of their wind, the GRBs are considered as possible accelerators of UHE-cosmic rays (Vietri (1995), Waxman (1995)) : for $`r>r_b`$, particles could reach the maximal energy beyond which they are no longer confined, $`ϵ_{max}=ZeBr`$, namely

$`ϵ_{max}=7.7\times 10^{21}\left({\displaystyle \frac{B}{10^7G}}\right)\left({\displaystyle \frac{r}{10^{12}cm}}\right)eV,`$
for a proton ($`Z=1`$).

### 5.1 Proton energy in the transparent stage

We examine the loss limitation suffered by the protons after crossing the ”proto”-sphere. After the equation (20), we have already given the value of the maximum energy due to the synchrotron loss at $`r_b`$; thus, it scales like :

$$ϵ_{syn}3.8\times 10^5\left(\frac{B(r_b)}{10^7G}\right)^{5/4}\left(\frac{r}{r_b}\right)^{(5\alpha 2)/4}GeV.$$
(28)
It will turn out that this value is above the threshold of the p$`\gamma `$-process, but is far below the range expected for getting UHE-Cosmic Rays. For $`\alpha =3/2`$, this range would be reached at a few 100 $`r_b`$.

The expansion loss can be significant. Indeed, the cut off energy is obtained by setting $`t_{acc}=t^{}=t/\mathrm{\Gamma }`$ in the co-moving frame, and the limitation, $`ϵ_{exp}`$, is such that

$$ϵ_{exp}10^5\left(\frac{\kappa _0}{10}\right)^3\left(\frac{\eta }{300}\right)\left(\frac{B(r_b)}{10^7G}\right)\left(\frac{r}{r_b}\right)^{1\alpha }GeV.$$
(29)
Whereas the synchrotron loss diminishes with distance, the expansion limitation becomes more and more severe (see figure 1). Moreover, it can easily be checked that the escape time due to transverse diffusion is always much longer than the expansion time.

### 5.2 The p$`\gamma `$-neutrino emission

A signature of VHE-protons acceleration is expected with the neutrino production resulting from the photo-production of pions, which is efficient to produce neutrinos through the $`\mathrm{\Delta }resonance`$ (the so-called GZK effect) :

$`\stackrel{(2/3)}{}p+\pi ^0p+\gamma +\gamma `$
$`p+\gamma `$ $`\mathrm{\Delta }^+`$
$`\stackrel{(1/3)}{}n+\pi ^+\mathrm{}p+e^++e^{}+\nu _e+\overline{\nu }_e+\nu _\mu +\overline{\nu }_\mu `$
The threshold of the process is such that $`ϵ_pϵ_\gamma \frac{1}{2}m_nm_\pi \frac{1}{4}m_\pi ^2`$, which leads to a threshold for the proton Lorentz factor $`\gamma _{th}10^4`$, the energy of the photons being in the keV range in the co-moving frame. As previously seen, $`\gamma _{max}`$ is at least ten times larger than $`\gamma _{th}`$. The GRBs are likely opaque to this process until the shells reach the radius $`r_{p\gamma }`$ such that

$$r_{p\gamma }=\left(\frac{\sigma _{p\gamma }E_\gamma }{\mathrm{\Omega }\overline{ϵ}_\gamma }\right)^{1/2},$$
(31)
where $`\overline{ϵ}_\gamma `$ is the average energy of the target photons in the observer frame (typically 1 MeV); which puts this p$`\gamma `$-”proto”-sphere at a rather large radius of $`1.5\times 10^{15}`$ cm for $`E_\gamma =0.1E`$.

In the collision frame, the energy of a pion, generated sufficiently above the threshold, is $`ϵ_\pi ^{}ϵ_p^{}ϵ_\gamma ^{}\gamma ϵ_\gamma (1\mathrm{cos}\theta )`$, where $`\gamma `$ is the Lorentz factor of the proton. The emitted neutrinos have an energy which is a fraction $`\alpha _0`$ ($`5\%`$) of the pion energy. Similarly to the statistical treatment of the pp-collisions, we write a simplified (with delta approximation instead of a function smoothed by angle averaging) kinetic equation for the isotropic distribution of the neutrinos generated by p$`\gamma `$-collisions in the co-moving frame, which is :

$$\frac{}{t}f_\nu +c\stackrel{}{n}.f_\nu =\xi _\nu c\sigma _{p\gamma }_{\gamma >\gamma _{th}}f_{}(\gamma )𝑑\gamma f_\gamma (ϵ_\gamma )\delta (ϵ_\nu \alpha _0\gamma ^2ϵ_\gamma )𝑑ϵ_\gamma .$$
(32)
Integrating both sides of the equation over the co-volume and the proper time, we easily get the number of emitted neutrinos :

$$N_\nu =\xi _\nu N_{}^>=\xi _\nu N_{}\left(\frac{1}{\gamma _{th}}\frac{1}{\gamma _{max}}\right),$$
(33)
where $`N_{}^>`$ is the number of cosmic rays above the p$`\gamma `$-threshold energy, the last result being obtained for a $`\gamma ^2`$ proton distribution. The integration, like in the case of pp-collisions, starts at the p$`\gamma `$-”proto”-sphere. The total energy radiated by neutrinos in the co-moving flow is also derived easily from the kinetic equation :

$$E_\nu ^{}\alpha _0\xi _\nu N_{}\mathrm{log}\left(\frac{\gamma _{max}}{\gamma _{th}}\right)\overline{ϵ}_\gamma ^{}.$$
(34)
The logarithm factor is obtained for a $`\gamma ^2`$ distribution of protons and can easily be modified. This result can also be rewritten in the observer frame as follows :

$$E_\nu \alpha _0\gamma _{th}\overline{ϵ}_\gamma \mathrm{log}\left(\frac{\gamma _{max}}{\gamma _{th}}\right)N_\nu .$$
(35)
The energy radiated in the form of p$`\gamma `$-neutrinos is $`E_\nu 10^7E`$. Numerical computation of such spectra has been done (see Mücke et al. (1998)). The kinetic equation also provides with the neutrino spectrum. Assuming $`f_{}`$ is a power law distribution and that the energy distribution of target photons is $`ϵ_\gamma f_\gamma ϵ_\gamma ^\alpha `$, we get a power law energy spectrum for the neutrinos, namely

$$ϵ_\nu ^2\frac{dN_\nu (>ϵ_\nu )}{dϵ_\nu }ϵ_\nu ^{\frac{s1}{2}},$$

in an energy range depending on the energy range of the protons since $`ϵ_\nu \alpha _0\gamma ^2ϵ_\gamma `$. These spectra are given in the co-moving frame. For the observer, they are Doppler beamed with the bulk Lorentz factor of the relativistic wind. However, the number of events is so low that it is still “virtual” to talk about a spectrum…

## 6 Discussion

In order to account for the non-thermal and highly variable gamma emission of GRBs, the fireball model and the internal shock model have been designed with baryonic load parameter $`\eta `$ that has been supposed large enough to get a relativistic wind achieving a large Lorentz factor $`\mathrm{\Gamma }\eta `$. However, it turns out that opacity effects could easily be significant at the beginning of the emission. We stress that point through a discussion involving a critical value $`\eta _{}`$ of the baryon load parameter. Indeed one emphasis of this paper is to analyse the opacity of the GRB to relativistic protons with respect to pion production by pp-collisions. Actually, the opacity condition relative to pp-collision is not far from the Compton opacity condition. Therefore it is reasonable to think that a significant fraction of the GRBs experience proton Fermi acceleration with efficient pp-collisions revealed by neutrino radiation.

During a stage of pp-opacity, we have shown that Fermi acceleration leads to conclusions about the pp-process that are very sensitive to the choice of the efficiency law of the acceleration. We have shown that the Bohm scaling assumption leads to an efficient proton acceleration that would be limited by synchroton loss, whereas the correct law governed by the turbulence spectrum (Kolmogorov law was used) leads to very different estimates and the pp-collisions process turns out to be the main limitation of the proton energy. Moreover, the excessive efficiency of the Bohm scaling would make the acceleration to deplete its energy reservoir in a time too short to maintain the gamma emission. The most important surprise raised in estimating the energy limitation in the radiation free stage. Indeed, the expansion limitation turns out to be drastic with the Kolmogorov law and maintains the proton energy below a few $`10^5`$ GeV.

The analysis of the paper indicates that a double neutrino emission can be expected with many GRBs, namely a stage of pp-neutrino emission followed by a stage of p$`\gamma `$-neutrinos. The number of emitted pp-neutrinos gives the amount of relativistic protons; and the number of p$`\gamma `$-neutrinos gives the number of protons above the threshold, which is about $`10^4`$ less. We proposed an analytical shape of the neutrinos spectra, as well as the photon spectrum generated by the $`\pi ^0`$-decay.
The magnetic field intensity is an important parameter that controls the proton acceleration and the synchrotron losses. The neutrino emissions significatively depend on its value at the crucial distance $`r_b`$. We took a high but still reasonable value of $`10^7G`$ at this distance; if we take less, the synchrotron limitation is less important, but the expansion limitation, which controls the highest energy of the protons in the GRB, becomes more severe, and also the pp-neutrino emission becomes less energetic. If we unduly take more, synchrotron losses dominate over pp-collisions and the acceleration is more efficient against the expansion losses, however not sufficient to get UHE-cosmic rays…This paper does not exclude the possibility of UHE-cosmic ray generation in GRBs. It simply states that its achievement with Bohm scaling is not reliable and leads to observational inconsistencies and that its achievement with Kolmogorov scaling is impossible…We think that there is another possibility (Pelletier (1999), Pelletier & Kersalé (2000)) that deserves a detailed investigation that we will present in a forthcoming paper.

Acknowledgement : The authors are grateful to Frédéric Daigne, Gilles Henri, Martin Lemoine and Eli Waxman for fruitful discussions.


format: text

text:
# Coulomb excitation of 11Be reexamined

( 28 November 2002)

## Abstract

The study of Coulomb excitation of 11Be to the first excited state in intermediate energy collisions with heavy targets is presented. The existing experimental data are reanalised by including the probability of projectile survival into the calculation of Coulomb excitation cross sections. The survival probabilities are calculated using a recently developed global optical model potential tailored in line with the double folding model. The extracted $`B(E1)`$ values for the transition $`\frac{1}{2}^+\frac{1}{2}^{}`$ in 11Be are found to be slightly larger than those obtained so far using the $`b_0`$-recipe in Coulomb excitation theory.

The Coulomb excitation process has been proven to be a useful tool to study nuclear structure (see, e.g., Ref.Ber01 ). It is an attractive test of the different nuclear models as the Coulomb interaction is very well known. Recently a number of experiments has been done to study the structure of light neutron rich nuclei, such as 11Be Ann95 ; Nak97 ; Fau97 , 15C, 17C, 19C and several oxygen isotopes, using the Coulomb excitation process Pra02 .

The isotope 11Be is one of the most studied halo nuclei nowadays. The reason why 11Be is an attractive system to study by Coulomb excitation lies in the peculiar structure of the spectrum of this nucleus. Both the $`\frac{1}{2}^+`$ ground state and the only bound excited state $`\frac{1}{2}^{}`$ are weakly bound and have energy difference of 0.32 MeV, thus making the $`E`$1 transition in 11Be the fastest known between bound states.

There have been several experiments that measured the excitation of the first excited state in 11Be and extracted the corresponding $`B(E1)`$ value. The adopted value ca 0.116 $`e^2`$fm2 was obtained by Millener et al. in Ref.Mil83  by averaging the results of three experiments on the lifetime of the excited state using a Doppler-shift technique.

The experiment done in GANIL Ann95  studied the inelastic scattering of 11Be on a lead target at a bombarding energy of 45 MeV/u. The obtained cross section was only 40$`\%`$ of the one predicted by calculations for pure Coulomb excitation. To explain this large difference it was proposed that the higher order effects may contribute, but calculations done in Ref.Ber95  showed that the cross section falls by only 4$`\%`$ if coupling to continuum was taken into account. The inclusion of monopole and quadrupole modes of nuclear excitation also performed in the work Ber95  resulted in the slight increase of the cross section less than 2$`\%`$.

Similar experiment was done by Nakamura et al. in Ref.Nak97  at $`E`$ = 64 MeV/u. The $`B(E1)`$ value extracted from the cross section assuming pure Coulomb excitation was comparable with the result of the Millener analysis. Later the inelastic scattering of 11Be was measured at MSU Fau97  on lead and gold targets for energies 60 and 58 MeV/u, respectively. The extracted $`B(E1)`$ value confirmed that of Nak97  and agreed, at least marginally, with the lifetime experiments. For summary of this see Table 2.

In the analysis performed in the mentioned papers, the authors used the formalism of pure Coulomb excitation and excluded nuclear processes approximately by using a low impact parameter cutoff $`b_0`$. In this work we show that the inclusion of the full fledged survival probability $`|S(b)|^2`$ in conjunction with $`b_0`$, though slightly improves the agreement with the data, it does allow one to perform a more realistic smooth cutoff based analysis.

In the calculation of the Coulomb excitation cross section, the important ingredient of the model is the minimum value of the impact parameter, $`b_0`$, from which the integration of the excitation amplitude is performed. This $`b_0`$ is roughly determined by the sum of the projectile and target radii. 11Be, being a halo nucleus, has a very diffuse structure, which makes the definition of the minimum impact parameter obscure.

To improve upon the calculation we multiply the probability of the Coulomb excitation at a given impact parameter by the survival probability $`|S(b)|^2`$ calculated for the system under consideration. In this work we take the S-matrix from the optical model. The calculation of the elastic S-matrices requires the knowledge of the densities of the projectile and the target. The halo density of 11Be was calculated in the framework of the particle-rotor model, which includes the excitation of the rotational 2+ state of the core 10Be Nun96 .

As was mentioned above, the contribution of nuclear excitation is small. This is so since nuclear effects are limited to a very small impact parameter region around the grazing value. In fact, as shown in Refs.Mot94 ; Ber94 ; Shy96  the Coulomb excitation was found to be by far the dominant piece of the cross section. Further, in Ref.Fau97  an experiment was done with light targets to study the importance of nuclear excitation. The results were 4.0 mb for carbon target and 1.7 mb for beryllium target. Since the nuclear cross section corresponds to the area of a ring around grazing $`b_0`$, it scales with $`A`$ as $`A^{1/3}`$. Accordingly, we find for the lead target approximately 10 mb nuclear excitation contribution implying a mere 3$`\%`$ effect. Therefore in the following we ignore the nuclear excitation effect.

The semiclassical model is usually used to describe Coulomb excitation at intermediate and high energies. This model assumes the straight line trajectory for the projectile and treats quantum mechanically the absorption of radiation by the nucleus. The formalism of this method was presented in Ref.Ald75  and later was extended to the relativistic case in Ref.Win79 . For the high energies the first order perturbation theory is a good approximation to calculate the amplitudes for Coulomb excitation. In the first order perturbation theory the process of Coulomb excitation can be described as emission and absorption of virtual photons Fer24 . We included only $`E1`$ multipolarity in our analysis since for the transition $`\frac{1}{2}^+\frac{1}{2}^{}`$ in 11Be the dipole multipolarity is the dominant one Ber88 .

Using the formalism of virtual photons the Coulomb excitation cross section has the following form:

$$\sigma _C=_0^{\theta _0}\frac{d\sigma _C}{d\mathrm{\Omega }}𝑑\mathrm{\Omega }=2\pi _{b_0}^{\mathrm{}}P(b)b𝑑b$$
(1)
with

$$P(b)=\frac{16\pi ^3}{9\mathrm{}c}B(E1)N_{E1}(\omega ,b),$$
(2)
and $`b_0=a\mathrm{cot}\frac{\theta _0}{2}`$, where $`a`$ is half the distance of closest approach for head-on Coulomb collision. $`N_{E1}(\omega ,b)`$ is the number of virtual photons given by

$$N_{E1}(\omega ,b)=\frac{Z_1^2\alpha }{\pi ^2}\left(\frac{\xi }{b}\right)^2\left(\frac{c}{v}\right)^2\left[K_1^2\left(\xi \right)+\frac{1}{\gamma ^2}K_0^2\left(\xi \right)\right],$$
(3)
where $`Z_1e`$ is the charge of the target, $`\xi =\frac{\omega b}{\gamma v}`$, $`\omega =E_{ex}/\mathrm{}`$, $`E_{ex}`$ is the excitation energy of the state in the projectile, $`v`$ is the relative energy, $`\gamma `$ is the relativistic factor and $`\alpha `$ is the fine structure constant. The behavior of the integrand in the Eq.(1) is determined by the impact parameter dependence of the modified Bessel functions of the zero and first order $`K_0`$ and $`K_1`$.

Coulomb recoil was taken into account using the method of Winther and Alder Win79  replacing $`b`$ in the expression for $`P(b)`$ by $`b^{}=b+\pi a/2\gamma `$. Further, we multiply $`P(b)`$ of Eq.(1) by the survival probability $`|S(b)|^2`$, calculated from the optical potential NLM3Y of Ref.Cha02  and then write

$$\stackrel{~}{P}(b)=P(b)|S(b)|^2.$$
(4)
The optical potential was calculated from the double folding model, using projectile and target densities and the effective nucleon-nucleon interaction. The effect of the exchange-related non-locality is taken fully into account. The local-equivalent, energy-dependent potential can be written as Cha02

$$V(r,E)=V_N(r,E)iW(r,E),$$

where

$$V_N(r,E)=V_F(r)e^{4v^2/c^2},$$

where $`c`$ is the speed of light and $`v`$ is the local relative velocity between the two nuclei

$$v^2(r,E)=\frac{2}{\mu }(EV_C(r)V_N(r,E)).$$
(5)
The local folding potential $`V_F(r)`$ is obtained in the usual way. For the Coulomb part, $`V_C(r)`$ we employ a similar folding prescription using the proton densities. The imaginary part of the optical potential was found to be

$$W(r,E)0.78V_N(r,E).$$
(6)
The determination of the number 0.78 is discussed in Ref.Alv02 . This number was obtained by fitting calculated elastic scattering cross sections with experimental data for the following systems: 12C+12C,16O,40Ca,90Zr,208Pb; 16O+208Pb, 40Ar+208Pb and for the wide range of energies. The calculations showed that the number 0.78 is approximately system-independent.

We employ the same potential here by using the appropriate halo+core density of 11Be. For the densities of the heavy targets and the two-parameter Fermi distributions were used

$$\rho (r)=\frac{\rho _0}{1+\mathrm{exp}\frac{rR_0}{a}},$$
(7)
where $`a=0.56`$ fm, $`R_0=1.31A^{1/3}0.84`$ fm and $`\rho _0`$ is determined by the normalization condition. These densities were obtained in Ref.Cha02  with the aim of providing a global description of the nuclear interaction, based on an extensive study involving charged distributions extracted from electron scattering experiments and theoretical densities calculated through the Dirac-Hartree-Bogoliubov model. For the density of the core of 11Be, 10Be, we used the Eq.(7) with $`a=0.5`$ fm and $`R_0=1.75`$ fm, which gives $`R_{rms}`$(10Be) = 2.3 fm.

We use the existing particle + core excitation model to describe the halo structure of a light halo nucleus such as 11Be Nun96 . The inclusion of the core excitation allows for the coupling between the collective degrees of freedom of the core and the orbital motion of the single neutron. In this calculation we assumed a quadrupole-deformed core (with deformation parameter $`\beta _2`$) and included only the ground state 0+ and the first excited rotational state 2+. The interaction between the valence particle and the core was described by a deformed Woods-Saxon potential. The spin-orbit interaction, $`V_{so}`$, is the standard undeformed one. To reproduce the adopted $`B(E`$1) of 11Be we adjusted the model parameters obtained in Ref.Nun96  and these parameters are presented in Table 1. The particle-core model with these parameters gives a radius for 11Be of 3.0 fm, the $`B(E1)=0.116`$ $`e^2`$fm2.

In Fig.1 we show the calculated squared moduli of the elastic S-matrices for 11Be + 208Pb at 59.7 MeV/u. It is seen that using the density with the halo reduces the survival probability.

In Ref.Nak94  the projectile survival probability was determined from the angle dependence of the Coulomb dissociation cross section. The very small number of experimental points (3 in the grazing region) with their rather large error bars is certainly insufficient to make meaningful comparison. In Ref.Nak94  the survival probability was approximated by the function with the shape $`1/(1+\mathrm{exp}((bb_0)/a))`$ with parameters $`b_0=12.3(1.2)`$ fm and $`a=0.9(0.6)`$ fm. In fact, one should include the size of the error bars in the presentation of the experimental survival probability. For the purpose of completeness we show in Fig.2 the experimentally determined survival probability of Nakamura et al. presented in the shaded area together with ours shown as the full line. It is clear that our calculation of the survival probability agrees reasonably well with Nakamura’s one.

The experimental data on Coulomb excitation of 11Be is summarized in Table 2. The fifth column shows the values of $`B(E1)`$ extracted from the corresponding Coulomb excitation cross section by assuming pure Coulomb excitation and choosing appropriate $`b_0`$, that is without including projectile survival probability. In the analysis of Nakamura et al.Nak97  Eq.(1) was used with $`b_0=12.3`$ fm which was determined from the impact parameter dependence observed in the Coulomb breakup of 11Be Nak94 . Fauerbach et al. in Ref.Fau97  used the similar formalism to extract $`B(E1)`$ values from the experimental cross sections with $`b_0`$ obtained from the aperture of the experimental setup.

In our calculation we use $`\stackrel{~}{P}(b)`$ in Eq.(1) and start the integration using the experimental value of $`b_0`$, which is determined by the aperture of experimental setup. Thus for the case of, e.g., 11Be + 208Pb at 59.7 MeV/u we integrated the probability shown by dotted line starting with $`b_0=11.5`$ fm.

The theoretical cross sections calculated in this work are presented in columns 6, 7 and 8 of Table 2. First we assumed the $`B(E1)`$ values from the previous analysis Nak97 ; Fau97  and calculated Coulomb excitation cross sections without taking into account Coulomb recoil and using density of Eq.(7). We obtained the cross sections slightly larger than the experimental for 11Be + 208Pb at 64 and 45 MeV/u as we integrated the probability from $`b_0`$ obtained from experimental aperture and not from grazing impact parameter as was done in Nak97 . Taking into account Coulomb recoil reduced the calculated cross section by approximately 2$`\%`$. The further reduction (approximately 2$`\%`$) in the calculated cross section was obtained by using the halo densities of 11Be.

In Fig.3 we show the impact parameter dependence of the Coulomb excitation probabilities $`P(b)`$ and $`\stackrel{~}{P}(b)`$ for the case of 11Be + 208Pb at 59.4 MeV/u. The inclusion of projectile survival probability reduces the Coulomb excitation probability for impact parameters from 11.5 fm and up to 17 fm. This results in the reduction of the cross section by 4$`\%`$ compared to the result of the pure Coulomb excitation with $`b_0=11.5`$ fm. For the case of 11Be + 197Au the the reduction of the calculated cross section was 5$`\%`$.

Finally, we extracted the revised values of $`B(E1)`$ using $`\stackrel{~}{P}(b)`$ and taking Coulomb recoil into account. We found no difference between the values of $`B(E1)`$ extracted in the work Nak97  and results of our analysis. The revised $`B(E1)`$ values obtained from the experimental cross sections of the work Fau97  are approximately 4$`\%`$ larger than those obtained in Fau97 . The revised $`B(E1)`$ values are presented in the last column of the Table 2.

In conclusion, a better treatment of nuclear absorption exemplified by the use of an appropriate survival probability in the calculation of the Coulomb excitation cross section of 11Be leads to slightly smaller cross-section and then a slightly (about 4$`\%`$) larger $`B(E1)`$ value. We expect similar effect in the Coulomb dissociation cross-section of 11Be and other halo nuclei.

This work was supported in part by FAPESP and the CNPq.


text:
# Nonequilibrium stabilization of charge states in double quantum dots

## Abstract

We analyze the decoherence of charge states in double quantum dots due to cotunneling. The system is treated using the Bloch-Redfield generalized master equation for the Schrieffer-Wolff transformed Hamiltonian. We show that the decoherence, characterized through a relaxation $`\tau _r`$ and a dephasing time $`\tau _\varphi `$, can be controlled through the external voltage and that the optimum point, where these times are maximum, is not necessarily in equilibrium. We outline the mechanism of this nonequilibrium-induced enhancement of lifetime and coherence. We discuss the relevance of our results for recent charge qubit experiments.

The loss of quantum coherence is a central paradigm of modern physics. It not only governs the transition between the quantum-mechanical and the classical world, but has recently also gained practical importance in the context of engineering quantum computing devices. Decoherence naturally occurs in small quantum systems coupled to macroscopic heat baths. A huge class of such baths generates Gaussian noise and can hence be mapped on an ensemble of harmonic oscillators as in the spin-boson model Leggett . This can even apply, if the fundamental degrees of freedom of the bath are fermions, as it is, e.g., the case if the bath is a linear electrical circuit EPJ ; ZimanyiPRL , which is producing Gaussian Johnson-Nyquist noise.

In this Rapid Communication, we study a generically different system: a double quantum dot coupled to electronic leads. Such systems are studied as realizations of quantum bits Blick ; UdoFrank ; Hayashi . The position (either left or right dot) of an additional spin-polarized electron is used as the computational basis of a charge qubit as realized in Ref. Hayashi . For another proposal of a charge qubit in semiconductors, see Ref. Hollenberg .

Our system simultaneously couples to two distinct reservoirs of real fermions. Other than oscillator bath models, this allows for the application of a voltage between these reservoirs as a new parameter for controlling decoherence. The voltage creates nonequilibrium between the baths, which to the best of our knowledge has not been studied yet in the literature on open quantum systems.

We study the dynamics of the reduced density matrix and identify the usual two modes of decoherence, dephasing and relaxation: Dephasing is the loss of phase information, manifest as the decay of coherent oscillations. This corresponds to the time evolution of the off-diagonal elements of the reduced density matrix in the energy basis. Relaxation is the process during which a quantum system exchanges energy with the environment and ends up in a stationary state. This is described through the evolution of the diagonal density matrix elements. We are going to show that, surprisingly, the charge states can be stabilized by external nonequilibrium, i.e., the relaxation time is longest at a well-defined finite voltage. We will show, that this working point is also very favorable in terms of dephasing but competes with another local maximum at zero voltage. Our theory should also have applications in other systems.

We consider a double quantum dot system with an appreciable tunnel coupling between the dots allowing for coherent molecular states in these systems Oosterkamp . The computational basis is formed by the position states of an additional spin-polarized electron Hayashi ; Hollenberg . A superposition can be created by variation of the interdot coupling. In order to stabilize the charge, the coupling of the dots to the two leads is driven to weak values and the dot is tuned to the Coulomb blockade regime Wilfred  where sequential tunneling is suppressed through the addition energy. Even then, the system couples to the environment through cotunneling Nazarov , the correlated exchange of two electrons with the external leads which ends up in a state with the same total charge as the initial one.

The relevant Hilbert space is spanned by four states written $`|i,j`$ denoting $`i,j`$ additional electrons on the left and right dot, respectively. $`|1,0`$ and $`|0,1`$ define the computational basis as they are energetically accessible, the closest virtual intermediate states for cotunneling are $`|v_0=|0,0`$ and $`|v_2=|1,1`$. This model applies if all relevant energy scales of the system ($`\epsilon _{\mathrm{as}}`$ and $`\gamma `$, see below) are much smaller than the charging energies to the next virtual levels ($`\epsilon _{v_2}`$ and $`\epsilon _{v_0}`$, also below), which in turn have to be smaller than the orbital excitation of the individual dots. This can be realized in small dots.

The total Hamiltonian of this system can be written as $`H=H_0+H_1`$ where $`H_0=H_{\mathrm{sys}}+H_{\mathrm{res}}`$ describes the energy spectrum of the isolated double-dot through $`H_{\mathrm{sys}}=\epsilon _{\mathrm{as}}(a^La^La^Ra^R)\epsilon _{v_0}\widehat{n}_{v_0}+\epsilon _{v_2}\widehat{n}_{v_2}+\gamma (a^La^R+a^Ra^L)`$ and the two electronic leads $`H_{\mathrm{res}}=_\stackrel{}{k}\epsilon _\stackrel{}{k}^Lb_\stackrel{}{k}^Lb_\stackrel{}{k}^L+_\stackrel{}{k}\epsilon _\stackrel{}{k}^Rb_\stackrel{}{k}^Rb_\stackrel{}{k}^R`$. The sum over the dot states only runs over the restricted Hilbert space described above, the $`a^{L/R}`$ act on the lowest additional electron state on either dot. The double-dot is characterized by the asymmetry energy $`\epsilon _{\mathrm{as}}=\epsilon _l\epsilon _r`$ between the individual dots and the interdot tunnel coupling $`\gamma `$. The virtual states $`|v_2`$ and $`|v_0`$ are separated from the system by energy differences $`\epsilon _{v_2}`$ (upper virtual level) and $`\epsilon _{v_0}`$ (lower virtual level). The tunneling part $`H_1=t_c_{\stackrel{}{k},n}(a_n^Lb_\stackrel{}{k}^L+a_n^Lb_\stackrel{}{k}^L)+t_c_{\stackrel{}{k},m}(a_m^Rb_\stackrel{}{k}^R+a_m^Rb_\stackrel{}{k}^R)`$ describes the coupling of each dot to its lead and will be treated as a perturbation. $`t_c`$ represents the tunnel matrix elements between the dots and the leads. It can be absorbed in a tunneling rate $`\mathrm{}\mathrm{\Gamma }=2\pi t_c^2N(ϵ_F)`$. This has to be chosen small such that the Kondo temperature is low $`T_\mathrm{K}T`$ and perturbation theory holds, e.g., $`\mathrm{\Gamma }=10^9`$ Hz. Figure 1 shows a sketch of the setup.

From now on, we use the basis of molecular states obtained by diagonalizing $`H_{\mathrm{sys}}`$ with splitting $`2\delta =2\sqrt{\epsilon _{\mathrm{as}}^2+\gamma ^2}`$. In order to treat cotunneling by leading-order perturbation theory, we rewrite $`H_1`$ using a Schrieffer-Wolff transformation Schrieffer . This removes the transitions to the virtual states and generates an effective Hamiltonian containing indirect transition terms between the molecular states. A more detailed description of our calculation is given in Refs. UdoFrank,  and Diplom, . The final Hamiltonian is of the form $`H=H_0+H_1^{}`$ where

$`H_1^{}`$ $`=`$ $`{\displaystyle \underset{c,d}{}}\alpha _c^{}\alpha _d[{\displaystyle \underset{Y,Y^{},\stackrel{}{k},\stackrel{}{k^{}}}{}}H_{\stackrel{}{k},\stackrel{}{k^{}},c,d}^{Y,Y^{}}b_\stackrel{}{k}^Yb_\stackrel{}{k^{}}^Y^{}+`$ (1)
$`+{\displaystyle \underset{Y,Y^{},\stackrel{}{k},\stackrel{}{k^{}}}{}}H_{\stackrel{}{k},\stackrel{}{k^{}},c,d}^{Y,Y^{}}b_\stackrel{}{k}^Yb_\stackrel{}{k^{}}^Y^{}],`$
where $`Y`$ and $`Y^{}`$ denote right or left lead, the $`\alpha `$s describe molecular states and the $`H_{\stackrel{}{k},\stackrel{}{k^{}},c,d}^{Y,Y^{}}`$ are given through 2nd order perturbation theory, i.e., they are of $`O(\mathrm{\Gamma }^2)`$. Note, that $`H_1^{}`$ conserves the particle number because it acts upon the double-dot by injecting and extracting an electron in a single step. The terms with $`YY^{}`$ transfer charge between different reservoirs. Note that Eq. (1) is a simple and generic Hamiltonian connecting a quantum system to two distinct particle reservoirs and is potentially relevant for systems other than quantum dots as well.

We study the open system dynamics in the case of a time-independent Hamiltonian with a fully general initial reduced density matrix. We use the well-established and controlled Bloch-Redfield Agyres , which has been demonstrated to work down to low temperature for certain models Hartmann . It involves a Born approximation in $`H_1^{}`$, i.e., it captures all cotunneling processes in lowest nonvanishing order. The Redfield equations Weiss  for the elements of the reduced density matrix $`\rho `$ in the eigenstate basis of $`H_{\mathrm{sys}}`$ (i.e., the molecular basis) read

$$\dot{\rho }_{nm}(t)=i\omega _{nm}\rho _{nm}(t)\underset{k,l}{}R_{nmkl}\rho _{kl}(t),$$
(2)
where $`\omega _{nm}=(E_nE_m)/\mathrm{}`$ and the Redfield tensor elements $`R_{nmkl}`$ are composed of golden rule rates describing different cotunneling processes, which are essentially independent due to the low symmetry of the system. Each contribution has a typical cotunneling structure UdoFrank ; Diplom . An overview of the most important processes is given below. $`n`$, $`m`$, $`k`$ and $`l`$ can be either $`+`$ (excited molecular state) or $``$ (molecular ground state) with according energies $`E_\pm `$. This type of perturbative analysis is only valid above the Kondo temperature $`T_K`$Haldane , which can be easily driven to low values by pinching off the tunneling barrier to the leads.

From the formal solution of Eq. (2) we can identify the relaxation and dephasing rates as

$`\mathrm{\Gamma }_r`$ $`=`$ $`\mathrm{Re}(R_{++++}+R_{})={\displaystyle \frac{1}{\tau _r}}`$ (3)
$`\mathrm{\Gamma }_\varphi `$ $`=`$ $`\mathrm{Re}(R_{++})=\mathrm{Re}(R_{++})={\displaystyle \frac{1}{\tau _\varphi }}.`$ (4)
The transition frequencies $`\omega _{\mathrm{nm}}`$ are weakly shifted.

Figure 2 shows a choice of processes entering the Redfield tensor. All processes contribute to dephasing, because the phase of an electron, which is injected from the reservoirs, is always random. Figures (a) and (b) illustrate relaxation processes. Only (b) contributes to the current, i.e., in general the relaxation rate must not be confused with the cotunneling current. In (c) and (d) two pure dephasing processes are presented, only (d) contributes to the current flow. In general, processes without current can emerge, if the cotunneling processes take place between a single lead and the two-state system (TSS). The observable current is then given by the difference of current-carrying processes in forward and backward direction. We have evaluated the rates entering, Eqs. (3)and (4) using the $`H_{\stackrel{}{k},\stackrel{}{k^{}},c,d}^{Y,Y^{}}`$s. Due to the high number of terms, details are not shown and will be given elsewhere Diplom .

We now turn to the discussion of our results, starting with the relaxation time $`\tau _r`$.

We observe in Fig. 3 that for an asymmetric TSS, i.e., for $`\epsilon _{\mathrm{as}}0`$, there is a pronounced peak of the relaxation time at $`V=\mathrm{sgn}(\epsilon _{\mathrm{as}})2\delta `$, where $`V=\mu _R\mu _L`$, i.e. the sign has to be chosen with opposite polarity to the asymmetry energy. This means in particular that the relaxation is minimal far away from equilibrium. This is the central result of our paper. It is most clearly visible for $`T2\delta `$, but obviously still dominates the calculated result for temperatures $`T2\delta `$, as it can be seen in the insets of Figures 3 and 4. In order to remain in the cotunneling regime, the voltages are still quite small as compared to the excitation energy to the next charge states $`\epsilon _{v_2}`$ and $`\epsilon _{v_0}`$, but on the order of the molecular level splitting, i.e. $`|V|2\delta \epsilon _{v_2},\epsilon _{v_0}`$. For quantum computation, achieving a maximum relaxation time is e.g. appreciable during read-out EPJ .

Although surprising, it can be understood from the analysis of the different rates, that $`V=0`$ does not necessarily imply the lowest relaxation rate. At $`V=0`$ there is no net current, i.e. no net exchange of particles in the ensemble average, however, this is achieved by the cancellation of finite currents of equal size in forward and backward direction. These currents are rather small UdoFrank  such that current heating is reduced to a minimum. To $`\tau _\mathrm{r}`$, Eq. (3), such current-carrying processes contribute with equal sign - the system relaxes no matter to which reservoir. On top of this, one also has to take into account the aforementioned current-less relaxation channels.

The appearance of the peaks as preferred stable points in Fig. 3 can be understood based on the analysis of the current-carrying processes, \[e.g. Figs. 2 (b) and (d)\] as schematically shown in Fig. 5. At low voltages, $`|V|<2\delta `$, the system relaxes into a thermal state close to the ground state. Relaxation takes place by spontaneous emission of energy into the environment and creation of an electron-hole pair in the leads. This pair can recombine through the electrical circuit which fixes the electrochemical potentials. This leads to electrical current. As the voltage is increased away from $`V=0`$, emission processes which lead to a current against the polarity of the source are suppressed, the others are enhanced, see Figs. 5 (I) and (II). Depending on the asymmetry of the double dot, i.e. on the weight of the excited state on the left and the right dot, this leads to an enhancement or a suppression of the rate. At $`|V|2\delta `$, the emission processes against the source are completely blocked: the dot relaxation does not provide enough energy to overcome the electromotive force. The rate vanishes linearily as a function of voltage reflecting the size of the available phase space for cotunneling, see Fig. 5.

At higher voltages, $`|V|2\delta `$, inelastic cotunneling Silvano  sets in, see Figs. 5 (III) and (IV): The source provides enough energy to even excite the double dot, creating a nonequilibrium steady-state population of the molecular levels. Hence, inelastic cotunneling provides a way for the dot to absorb energy from the environment even at low temperature. This process can be experimentally identified by a sharp increase of the current Silvano ; UdoFrank .

Hence, at $`V=\pm 2\delta `$, three of the four processes depicted in Fig. 5 vanish at low temperatures, whereas at $`V=0`$ only two vanish. The linear voltage-dependence of the rates leads to the rather sharp cusps seen in Fig. 3. This behavior is smeared out at higher temperatures by thermal fluctuations. The peak height is set by the remaining processes: Energy emission with the source and currentless relaxation, Fig. 2 (a). As explained above, the relative weight of the former strongly depends on the weight of the excited molecular state on the individual dots and thus is responsible for the strong asymmetry of the peaks in Fig. 3 for different asymmetry energies.

Finally, we analyze the properties of the dephasing time $`\tau _\varphi `$ as a function of the bias voltage.

The total dephasing rate contains relaxing as well as flipless (“elastic”) processes. We hence observe in Fig. 4 a peak structure at $`V=\mathrm{sgn}(\epsilon _{as})2\delta `$ as in the relaxation time, Fig. 3, and a similar peak at $`V=0`$. The latter can be understood from the suppression of flipless processes (energy exchange 0)in an analogous way to the relaxation peak in Fig. 3 (energy exchange $`2\delta `$). At low asymmetry energy $`ϵ_{\mathrm{as}}<\gamma `$, the dephasing time at $`V=0`$ is longest. At high asymmetry $`ϵ_{\mathrm{as}}>\gamma `$ and at the nonequilibrium working point $`V=\mathrm{sgn}(ϵ_{\mathrm{as}})2\delta `$, $`\tau _\varphi `$ is even longer. In general, this indicates the existence of two preferable working points for quantum computation: One in equilibrium, the other again far from equilibrium. As also already seen in the inset of Fig. 3, the voltage dependence at higher temperature is here smeared out and the peaks merge.

A measurement of the relaxation and dephasing times should be feasible either by a time-resolved measurement of $`\sigma _z(t)`$, e.g., through a single-electron transistor or point contact Elzerman , the saturation broadening method Marlies  or resonance schemes such as proposed in Ref. \[Engel, \] for spins.

Note that parts of the double-dot literature focus on decoherence through phonons or photons (see Refs. Qin, ; Oosterkamp, ; Hazelzet, ; Brandes, ), whereas we focus on the cotunneling, which becomes relevant when phonons are suppressed by a cavity Eva . If the spin in a dot is used as qubit Elzerman ; Loss , cotunneling serves as an indirect contribution to decoherence.

To conclude, we have studied the decoherence of charge states in a double quantum dot due to cotunneling. We have shown that decoherence can be controlled through a bias voltage $`V`$ (and thus creating a nonequilibrium situation) between the two fermionic baths. In particular, the optimum working point for read-out and potentially also for operation of the qubit can be in an out-of-equilibrium situation at a voltage $`V=\mathrm{sgn}(\epsilon _{\mathrm{as}})2\delta `$. We have given a consistent physical interpretation of our findings in terms of stability and phase space. This effect of stabilization through nonequilebrium should potentially be significant for other qubit candidates as well.

We thank J. von Delft, L. Borda, J. König, M.J. Storcz, A.W. Holleitner, A.K. Hüttel and E.M. Weig for clarifying discussions. Work supported by SFB 631 of the DFG and in part by the National Security Agency (NSA) and Advanced Research and Development Activity (ARDA) under Army Research Office (ARO) contract number P-43385-PH-QC.


text:
# Entanglement generation in continuously coupled parametric generators

## I Introduction

Quantum entanglement and its consequences puzzled physicists since the famous paper by Einstein, Podolsky and Rosen Einstein35 . It was recognized in recent years that the entanglement is not only an intriguing and mind boggling feature of quantum mechanics but also a crucial and extremely useful resource for information processing. The various protocols relying on the entanglement include quantum teleportation Bennett93 , entanglement swapping Pan98 , dense coding Bennett92 , quantum cryptography Ekert91  and quantum computing Barenco95 . The efficiency of quantum information processing significantly depends on the degree of entanglement of the quantum state shared by the parties involved in a given protocol. It is therefore highly desirable to establish reliable sources of pure strongly entangled states.

Quantum optics provides a natural and convenient framework for the generation of entangled states, their manipulation and measurement. The above-mentioned protocols have been originally established for the quantum systems with two-dimensional Hilbert spaces – the qubits. The polarization entangled pairs of photons generated by means of the spontaneous parametric downconversion proved to be a very good source of entangled qubits and they were employed in a large number of experimental demonstrations of quantum information processing.

Recently, however, a significant attention has been paid to the quantum information processing (QIP) with systems whose Hilbert space is infinite dimensional. A typical example of such a system is a single mode of the optical field whose Hilbert space is spanned by the infinite (but countable) number of the Fock states $`|n`$. It is in the spirit of this latter approach that one considers the quadrature operators $`\widehat{x}`$, $`\widehat{p}`$ and one speaks about QIP with continuous variables (CV). Similarly as in the case of qubits, the parametric downconversion provides a source of CV entanglement. Here, the relevant entangled state is the so-called two-mode squeezed vacuum (sometimes also called the twin beam) that can be either prepared by nondegenerate downconversion or via mixing of two single-mode squeezed states on a balanced beam splitter. This last-mentioned approach has been in fact used in the experiment on the teleportation of continuous variables Furusawa98 . The whole setup for the generation of the CV entangled state is certainly nontrivial and quite complicated. A natural question arises whether it could be possible to integrate the generation of the single mode squeezed states and the mixing on a beam splitter into a single small and compact device that would serve as a source of the entanglement.

In this paper we show that the nonlinear optical couplers that have been thoroughly investigated during recent years both theoretically and experimentally (see Perina00 ; Fiurasek01  and references therein) can serve exactly for that purpose. We theoretically investigate the generation of squeezing and entanglement in the nonlinear optical coupler consisting of a pair of nonlinear waveguides operating by means of the degenerate parametric downconversion. The waveguides are placed close together so that the two modes can overlap and exchange energy via linear coupling. We study the dependence of the generalized squeezing as well as CV entanglement on the coupler’s parameters. The entanglement is quantified by the logarithmic negativity Vidal02 , which is an easily computable entanglement monotone. Finally, the optimum phase matching condition that yields the maximal entanglement is analyzed.

## II The model and the equations of motion

The device under study is schematically depicted in Fig. 1. The nonlinear coupler consists of two linearly coupled waveguides where the degenerate parametric downconversion takes place. A ‘blue’ photon with frequency $`2\omega `$ may split into a pair of identical ‘red’ photons and also an inverse process may occur. Assuming a strong coherent pumping of the modes with frequency $`2\omega `$ (second-harmonic modes) and a linear coupling between the modes with frequency $`\omega `$ (fundamental modes) in the lossless and nondispersive media, we can describe such a device by the following effective interaction momentum operator Perina96

$$\widehat{G}_{\mathrm{int}}=\mathrm{}\left(g_A^{}\widehat{A}_F^2+g_B^{}\widehat{B}_F^2+2g_L\widehat{A}_F\widehat{B}_F^{}+\text{h.c.}\right),$$
(1)
where $`\widehat{A}_F(\widehat{A}_F^{})`$ is annihilation (creation) operator of the fundamental mode in the first waveguide and $`\widehat{B}_F(\widehat{B}_F^{})`$ is annihilation (creation) operator of the fundamental mode in the second waveguide. The constant $`g_A(g_B)`$ is a product of the strong coherent amplitude of the second-harmonic mode in the first (second) waveguide and the nonlinear coupling constant, which is proportional to the quadratic susceptibility $`\chi _A^{(2)}`$ ($`\chi _B^{(2)}`$). The constants $`g_A`$ and $`g_B`$ determine the efficiency of the downconversion processes. Linear coupling constant denoted as $`g_L`$ describes the energy exchange between the waveguides by means of coupling via the evanescent waves. The symbol $`\mathrm{}`$ denotes the reduced Planck constant and h.c. represents Hermitian conjugate terms.

We shall work in the interaction picture where the spatial evolution of an operator $`\widehat{O}`$ is governed by the Heisenberg-like equation (see Perina00 ; Fiurasek01  and references therein)

$$i\mathrm{}\frac{d\widehat{O}}{dz}=[\widehat{G}_{\mathrm{int}},\widehat{O}].$$
(2)
Here $`z`$ is a spatial coordinate along the direction of propagation, $`\widehat{G}_{\mathrm{int}}`$ is an interaction momentum operator and $`[,]`$ denotes a commutator. Using the momentum operator (1), we obtain the Heisenberg equations of motion in the interaction picture of the form:

$`{\displaystyle \frac{d\widehat{A}_F}{dz}}`$ $`=`$ $`2ig_A\widehat{A}_F^{}+2ig_L^{}\widehat{B}_F,`$
$`{\displaystyle \frac{d\widehat{B}_F}{dz}}`$ $`=`$ $`2ig_B\widehat{B}_F^{}+2ig_L\widehat{A}_F.`$ (3)
For our purposes it is convenient to deal with the quadrature components $`\widehat{x}`$, $`\widehat{p}`$ rather than with the creation and annihilation operators. The quadratures can be directly measured via balanced homodyne detection. Recall that the quadrature components can be expressed as linear combinations of the creation and annihilation operators as follows:

$`\widehat{x}_A`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{2}}\left(\widehat{A}_F^{}+\widehat{A}_F\right),\widehat{p}_A=i{\displaystyle \frac{\sqrt{2}}{2}}\left(\widehat{A}_F^{}\widehat{A}_F\right),`$
$`\widehat{x}_B`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{2}}\left(\widehat{B}_F^{}+\widehat{B}_F\right),\widehat{p}_B=i{\displaystyle \frac{\sqrt{2}}{2}}\left(\widehat{B}_F^{}\widehat{B}_F\right).`$
We further express the (generally complex) coupling constants in terms of their amplitudes and phases,

$$g_j=|g_j|\text{exp(}i\phi _j\text{)}\text{for}j=L,A,B.$$
(5)
With these definitions at hand, we can derive from Eq. (II) the Heisenberg equations for the quadrature operators. We write them down in the compact matrix form,

$$\frac{d\widehat{\xi }}{dz}=𝐌\widehat{\xi },$$
(6)
where

$$𝐌=\left(\begin{array}{cccc}S_A& C_A& S_L& C_L\\ C_A& S_A& C_L& S_L\\ S_L& C_L& S_B& C_B\\ C_L& S_L& C_B& S_B\end{array}\right),\widehat{\xi }=\left(\begin{array}{c}\widehat{x}_A\\ \widehat{p}_A\\ \widehat{x}_B\\ \widehat{p}_B\end{array}\right),$$
(7)
and

$$S_j=2|g_j|\mathrm{sin}\left(\phi _j\right),C_j=2|g_j|\mathrm{cos}\left(\phi _j\right).$$
(8)
Eq. (6) represents a system of linear differential equations with constant coefficients that can be solved very simply by determining the eigenvalues and eigenvectors of the matrix $`𝐌`$. Explicit analytical expressions for the solution can be found in Perina96 . Here we write the solution in a formal matrix notation

$$\widehat{\xi }(z)=𝐒(z)\widehat{\xi }(0),𝐒(z)=e^{𝐌z}.$$
(9)
This evolution represents a linear canonical transformation of the quadrature operators. All possible linear canonical transformations of $`N`$ modes form a group of symplectic (also called Bogolyubov) transformations $`Sp(2N,R)`$. In our case, for each $`z`$ the matrix $`𝐒(z)`$ is an element from the four dimensional representation of the symplectic group $`Sp(4,R)`$.

One of the most important features of the symplectic transformations is that they transform Gaussian states again onto Gaussian states. The Gaussian states are those whose Wigner function is a Gaussian. Such states are thus fully characterized by the first and second moments of the quadrature operators. Moreover, the first moments (mean values) of the quadratures can be always set to zero via appropriate displacements applied locally to each mode. Since we are interested in entanglement properties of the state generated in the coupler, we can focus solely on the evolution of the second moments.

## III The generalized squeeze variance

It is convenient to arrange the second moments into the covariance matrix $`𝐕`$ whose elements are defined as follows Simon94 :

$$V_{jk}=\frac{1}{2}\left(\mathrm{\Delta }\widehat{\xi }_j\mathrm{\Delta }\widehat{\xi }_k+\mathrm{\Delta }\widehat{\xi }_k\mathrm{\Delta }\widehat{\xi }_j\right).$$
(10)
Here $`\mathrm{\Delta }\widehat{\xi }_j=\widehat{\xi }_j\widehat{\xi }_j`$ and $`\widehat{\xi }_j=\text{Tr}(\widehat{\rho }\widehat{\xi }_j)`$. In our analysis of the coupler we shall assume that all processes are spontaneous, i.e., both fundamental modes are initially in the vacuum state. The covariance matrix corresponding to this input state is proportional to the identity matrix,

$$𝐕(0)=\frac{1}{2}\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right).$$
(11)
The evolution of the covariance matrix corresponding to the symplectic transformation (9) is given by the following formula:

$$𝐕(z)=𝐒(z)𝐕(0)𝐒^T(z),$$
(12)
where $`T`$ stands for the transposition.

Quantum entanglement is one of the manifestations of the nonclassical features of the quantum states. The CV quantum state is said to have a classical analogue if its density matrix can be written as a convex mixture of projectors on coherent states. Such states cannot be entangled and they can exhibit only classical correlations. Since the nonclassicality is a necessary prerequisite for the entanglement Buzek02  it is of interest to investigate whether and how much nonclassical are the states generated in the nonlinear coupler. This issue has been addressed in many recent papers devoted to the analysis of the quantum-statistical properties of light propagating in nonlinear couplers Perina00 ; Fiurasek01 . Various kinds of nonclassical behaviour have been predicted, such as the generation of light with sub-Poissonian photon number statistics, and single-mode or two-mode squeezed light.

For Gaussian states, the nonclassicality criterion is particularly simple because these states are nonclassical if and only if (iff) they are squeezed. More formally, we say that the (generally $`N`$-mode) Gaussian state is nonclassical iff there exists a quadrature $`\widehat{X}`$ that is a linear combination of the quadratures $`\widehat{x}_j`$ and $`\widehat{p}_j`$ of the $`N`$ involved modes such that the variance of $`\widehat{X}`$ is below the coherent-state level $`1/2`$. The smallest variance obtained as minimum over the variances of all possible linear combinations of the quadratures $`\widehat{x}_j`$ and $`\widehat{p}_j`$ is called the generalized squeeze variance $`\lambda `$ Simon94 . This variance can be calculated as the lowest eigenvalue of the covariance matrix (12)

$$\lambda (z)=\mathrm{min}\{\mathrm{eig}[𝐕(z)]\}.$$
(13)
The Gaussian state is nonclassical and squeezed iff $`\lambda <1/2`$.

In our discussion of the coupler’s behavior we distinguish two main regimes of operation. In the first regime, the linear coupling dominates over the nonlinear interactions

$$2|g_L|>|g_A+g_B|$$
(14)
and the coupler operates below the threshold. If the opposite inequality holds then the nonlinear interactions in the waveguides become dominant and the coupler operates above the threshold. The threshold condition ($`2|g_L|=|g_A+g_B|`$) was determined by analyzing the eigenvalues of the matrix $`𝐌`$ (see also Ref. Fiurasek00 ). These eigenvalues are complex below the threshold, which leads to oscillatory dynamics typical for the linear coupling. Above the threshold, all eigenvalues are real and the dynamics of the coupler is reminiscent of a pure amplification process.

The behavior of the coupler also significantly depends on the phase matching condition. It can be inferred from the analytical expressions for the symplectic transformation $`𝐒(z)`$ that there is a single effective phase difference

$$\mathrm{\Delta }\varphi =\varphi _A\varphi _B+2\varphi _L$$
(15)
that controls the dynamics of the coupler.

In figure 2 we plot the dependence of the generalized squeeze variance $`\lambda `$ on the length $`z`$ of the symmetric coupler ($`|g_A|=|g_B|`$) that operates below the threshold. The three curves correspond to three different choices of the phase difference $`\mathrm{\Delta }\varphi `$. For $`\mathrm{\Delta }\varphi =0`$ we observe a periodic evolution of $`\lambda (z)`$. In this specific case, the eigenvalues of $`𝐌`$ are purely imaginary, the dynamics is purely oscillatory, and the quantum state of the two modes periodically returns to the initial vacuum state. For the two other choices $`\mathrm{\Delta }\varphi =\pi /2`$ and $`\mathrm{\Delta }\varphi =\pi `$ the parametric amplification sets on and the squeeze variance exponentially decreases with growing length of the coupler. The squeezing is fastest for $`\mathrm{\Delta }\varphi =\pi `$. For comparison, we plot in figure 3 the function $`\lambda (z)`$ for symmetric coupler operating above the threshold. It follows that in this case the phase difference almost does not influence the dynamics of squeezing and the three curves shown in Fig. 3 corresponding to three different phase shifts $`\mathrm{\Delta }\varphi `$ almost coincide.

For the asymmetric coupler ($`|g_A||g_B|`$) below the threshold, the eigenvalues of $`𝐌`$ have nonzero both imaginary and real parts giving rise to periodically modulated exponential decrease of $`\lambda (z)`$. The modulations smooth out as $`\mathrm{\Delta }\varphi `$ goes to $`\pi `$. Above the threshold $`\lambda (z)`$ exhibits qualitatively the same behavior as in the symmetric case.

## IV Separability and entanglement

Having discussed the squeezing of light in the coupler, we now turn our attention to the generation of entangled states. A state $`\widehat{\rho }_{AB}`$ of two subsystems $`A`$ and $`B`$ is entangled iff $`\widehat{\rho }_{AB}`$ cannot be written as a convex mixture of product states

$$\widehat{\rho }_{AB}\underset{j}{}p_j\widehat{\rho }_{A,j}\widehat{\rho }_{B,j},p_j>0,$$
(16)
where $`\widehat{\rho }_{A,j}`$ and $`\widehat{\rho }_{B,j}`$ are states of subsystems $`A`$ and $`B`$. Since we assume that the fundamental modes are initially in pure vacuum states and the evolution is unitary, the two-mode state is pure for all $`z`$. It holds that every pure bipartite state $`|\psi _{AB}`$ that is not a product state is entangled. While it is easy to decide whether a pure state is separable or not, this problem is much more difficult for the general mixed bipartite state. Several separability criteria have been proposed in the literature. Among them, the most powerful and important is the Peres-Horodecki (PH) criterion Horodecki95 ; Peres96 ; Horodecki96 . According to this criterion, a quantum state $`\widehat{\rho }_{AB}`$ is entangled if the partially transposed density matrix

$$(\widehat{\rho }_{AB})_{m\alpha ,n\beta }^{T_A}(\widehat{\rho }_{AB})_{n\alpha ,m\beta }$$
(17)
is nonpositive, i.e., it has at least one negative eigenvalue. In general, this is only a sufficient condition on entanglement and there exist entangled states with positive partially transposed density matrix (the so-called bound entangled states). However, in case of two-mode bipartite Gaussian states this criterion turns out to be both sufficient and necessary Duan00 ; Simon00 .

The PH criterion is qualitative, it tells us whether the state is entangled or separable, but it does not quantify the amount of the entanglement present in that state. For this purpose, various entanglement measures have been proposed. The PH criterion suggests that a function of the negative eigenvalues of the partially transposed density matrix could be a good entanglement measure. Indeed, such measure can be constructed and it is called the logarithmic negativity Vidal02

$$E_N(\widehat{\rho })=\mathrm{log}_2[1+2𝒩(\widehat{\rho })].$$
(18)
Here the negativity $`𝒩(\widehat{\rho })`$ is the sum of the absolute values of the negative eigenvalues $`\mu _j`$ of the partially transposed matrix $`\widehat{\rho }_{AB}^{T_A}`$

$$𝒩(\widehat{\rho })=\underset{j}{}|\mu _j|.$$
(19)
The logarithmic negativity $`E_N`$ has the great advantage that it can be easily computed for any bipartite Gaussian state. In fact, analytical formula for $`E_N`$ has been derived for the general mixed two-mode Gaussian state. The details can be found in Vidal02 , here we only mention the main results and formulas.

It is convenient to decompose the covariance matrix as follows

$$𝐕=\left(\begin{array}{cc}𝐀& 𝐂\\ 𝐂^T& 𝐁\end{array}\right),$$
(20)
where $`𝐀`$ and $`𝐁`$ are the covariance matrices of the modes on Alice’s and Bob’s sides, respectively, and $`𝐂`$ contains the intermodal correlations. We must calculate the symplectic spectrum $`(c_1,c_2)`$ of $`𝐕`$ that can be obtained by solving the biquadratic equation (only positive roots are considered)

$$\zeta ^4+\left(\text{det}𝐀+\text{det}𝐁2\text{det}𝐂\right)\zeta ^2+\text{det}𝐕=0.$$
(21)
The logarithmic negativity is a function of the symplectic eigenvalues

$$E_N=\underset{j=1}{\overset{2}{}}F(c_j),$$
(22)
where the function $`F(c)`$ reads

$$F(c)=\{\begin{array}{cc}0& \text{for}\mathrm{\hspace{0.17em}2}c1\\ \mathrm{log}_2(2c)& \text{for}\mathrm{\hspace{0.17em}2}c<1.\end{array}$$
(23)
The evolution of the logarithmic negativity is plotted in Fig. 4. The parameters are the same as in Fig. 2, so the coupler is below the threshold. As could have been expected, for $`\mathrm{\Delta }\varphi =0`$ we get a purely oscillatory evolution and $`E_N`$ periodically attains its maximum. The oscillations are clearly present also for $`\mathrm{\Delta }\varphi =\pi /2`$ but the general trend is that $`E_N`$ grows with $`z`$. This is in agreement with the fact that for this detuning the state becomes more and more squeezed with growing $`z`$, cf. Fig. 2. Remarkably, $`E_N(z)=0`$ for all $`z`$ when $`\mathrm{\Delta }\varphi =\pi `$. For this particular phase difference, the coupler behaves like two decoupled single mode squeezers that do not interact at all Fiurasek00 . Figure 5 shows the behavior of $`E_N(z)`$ when the coupler operates above the threshold. We can see that the oscillatory evolution is replaced with a monotonic growth of $`E_N(z)`$ with $`z`$. However, the amount of entanglement significantly depends on $`\mathrm{\Delta }\varphi `$ even when the coupler is above the threshold. This is in marked difference with the generalized squeeze variance that is to a large extent insensitive to the choice of $`\mathrm{\Delta }\varphi `$, see Fig. 3.

In the case of the asymmetric coupler above the threshold $`E_N`$ grows monotonously with growing $`z`$. For small $`z`$ the rapidity of the growth decreases as $`\mathrm{\Delta }\varphi `$ goes from $`0`$ to $`\pi `$. Below the threshold $`E_N`$ behaves similarly as in the symmetric case. The only difference occurs for $`\mathrm{\Delta }\varphi =\pi `$ when $`E_N`$ changes periodically and no decoupling of the waveguides takes place.

The results depicted in Figs. 4 and 5 indicate that for each $`z`$ and a given set of parameters, we may tune the phase difference such as to maximize the entanglement. We have performed numerical calculations and determined the optimal phase difference $`\mathrm{\Delta }\varphi _{\mathrm{opt}}`$ and the corresponding maximal achievable $`E_N`$. The results for the coupler operating below the threshold are presented in Fig. 6. The $`\mathrm{\Delta }\varphi _{\mathrm{opt}}`$ exhibits an interesting behavior. Initially, it is optimal to set $`\mathrm{\Delta }\varphi _{\mathrm{opt}}=0`$ and this is optimal up to certain length $`z_0`$. At this distance, the optimal phase difference becomes to deviate from $`0`$ and $`\mathrm{\Delta }\varphi _{\mathrm{opt}}`$ oscillates with $`z`$. The figure 6 suggests that for large $`z`$, $`\mathrm{\Delta }\varphi _{\mathrm{opt}}`$ approaches some fixed asymptotic value. We note that similar results were obtained also for the coupler above the threshold.

## V Conclusions

In this paper we have explored the controllable compact source of squeezing and entanglement formed by a nonlinear optical coupler composed of two nonlinear waveguides operating by means of spontaneous degenerate parametric downconversion. Firstly, the behavior of the generalized squeeze variance in the below the threshold and above the threshold regimes has been analyzed. Secondly, by calculating the logarithmic negativity the possibility of the entanglement generation in the coupler has been demonstrated. In particular, it was shown that the amount of the entanglement can be controlled via the effective phase difference (15). Finally, with the help of the numerical calculation the optimal effective phase difference providing maximum entanglement at the output of the coupler has been found. Our theoretical analysis thus clearly illustrates the potential utility of the coupler for controlled generation of CV entanglement.

## VI Acknowledgments

This work was supported by Project LN00A015 and Research Project No. CEZ: J14/98 of the Czech Ministry of Education and by the EU grant under QIPC Project IST-1999-13071 (QUICOV).