text:
# Optimal approximate reversal of quantum operations on a single qubit

## Abstract

We demonstrate how insights gained from reformulating the problem of quantum teleportation into one of reversing quantum operations, and designing optimum completely positive maps for teleportation, can enable one to explore optimal approximate reversal of quantum operations on a single qubit. In particular, we show that the optimal approximate reversal of a generalized depolarizing channel can be achieved using only unitary transformations. We also show that for a quantum channel, which reveals some information about the input state, extremal completely positive maps and not unitary transformations yield optimal approximate reversal.

The formalism of quantum operations, described in detail by Kraus Kraus , describes the most general possible state change in quantum mechanics. In this formalism there is an input state $`\rho _{in}`$ and an output state $`\rho _{out}`$, which are connected by a map

$$\rho _{in}\rho _{out}=\frac{(\rho _{in})}{\mathrm{tr}[(\rho _{in})]}.$$

The map is determined by a quantum operation $``$, a linear, trace-decreasing map that preserves complete positivity. The most general form for $``$ can be shown to be Kraus

$$(\rho _{in})=\underset{k}{}A_k\rho _{in}A_k^{},$$

$$\underset{k}{}A_k^{}A_kI.$$
(1)
The Kraus operators $`A_k`$ completely specify the quantum operation $``$. In the special case of unitary evolution $`U`$ experienced by a closed quantum system, there is only one term in the operator-sum representation (1), $`A_1=U`$, leaving us with

$$\rho _{in}\rho _{out}=U\rho _{in}U^{}.$$

Unitary evolution is invertible, that is, the original input state $`\rho _{in}`$ can be recovered with certainty by subjecting the output state $`\rho _{out}`$ to the inverse unitary evolution $`U^1`$. This is not true for state change occuring in open quantum systems. The reversal of general quantum operations therefore becomes an important issue in quantum information processing (see Ref. Caves  and references therein). A quantum operation $``$ is “perfectly” reversible on a subspace $`𝒢`$ of the total state space $``$ if there exists a trace-preserving quantum operation $``$, acting on the total state space $``$, such that for all input states $`\rho _{in}`$ whose support lies in $`𝒢`$, Caves

$$\rho _{in}=\left(\frac{(\rho _{in})}{\mathrm{tr}[(\rho _{in})]}\right)=\frac{(\rho _{in})}{\mathrm{tr}[(\rho _{in})]}.$$

The necessary and sufficient conditions for a general quantum operation to be perfectly reversible were derived in Ref. Caves .

In this paper, we consider the problem of optimal approximate reversal of a quantum operation $``$ on the whole state space $``$ instead. From hereon, reversal means approximate reversal. We show how the techniques used in Ref. Rehacek  to identify and design optimum completely positive maps for quantum teleportation, can be employed to tackle our problem. The optimality of a given $``$ will be judged by the fidelity Jozsa  $`F(\rho _{in},(\rho _{in}))`$ between $`\rho _{in}`$ and $`(\rho _{in})`$ averaged over an isotropic distribution of input states $`\rho _{in}`$. We show explicitly that the action of a generalized depolarizing channel on a single qubit can be reversed optimally using only unitary transformations. We also consider a quantum operation which reveals information about the input state. In this case, optimal reversal requires more general quantum operations which are not unitary.

We begin with the observation Nielsen  that one can recast the problem of achieving optimal quantum teleportation Bennett  into one of optimal reversal of quantum operations. A general teleportation scheme involves a sender, Alice, and a receiver, Bob. Alice is in possession of two $`n`$-level quantum systems, the input system $`1`$, and another system $`2`$ arbitrarily entangled with a third $`n`$-level target system $`3`$ in Bob’s possession.

Initially the composite system $`123`$ is prepared in a state with density operator $`\stackrel{~}{\rho }_1\chi _{23}`$, where $`\stackrel{~}{\rho }_1`$ is an unknown state of the input system $`1`$, and $`\chi _{23}`$ is an arbitrary entangled state of systems $`2`$ and $`3`$. Since the systems $`1`$ and $`3`$ are identical and thus have the same state space, a one-to-one correspondence from the state space of the composite system onto itself can be established by a unitary swap operator $`U_{1(2)3}`$, which acts on product states according to

$$U_{1(2)3}(|\stackrel{~}{a}_1|b_2|c_3)=|\stackrel{~}{c}_1|b_2|a_3,$$

swapping the states of systems $`1`$ and $`3`$, while leaving system $`2`$ alone. $`U_{1(2)3}`$ obviously satisfies $`(U_{1(2)3})^2=I_{123}`$, the identity operator on the composite system, and $`U_{1(2)3}^{}=U_{1(2)3}`$. When extended to operators $`Q_{123}`$ on the composite system, the correspondence becomes

$$\stackrel{~}{Q}_{123}Q_{123}=U_{1(2)3}\stackrel{~}{Q}_{123}U_{1(2)3}^{}.$$

It follows that

$$\stackrel{~}{\rho }_1\chi _{23}=U_{1(2)3}(\stackrel{~}{\chi }_{12}\rho _3)U_{1(2)3}^{},$$
(2)
where $`\stackrel{~}{\chi }_{12}`$ is the counterpart of $`\chi _{23}`$.

To teleport the input state $`\stackrel{~}{\rho }_1`$ to Bob’s target system $`3`$, Alice performs a generalized measurement on systems $`1`$ and $`2`$. This generalized measurement is described by operators $`\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3`$, where $`\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}`$ are Kraus operators on the joint system $`12`$, $`i`$ labels the outcome of the measurement, and

$$\underset{i}{}\underset{j}{}\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}=I_{12}.$$

If Alice’s measurement has outcome $`i`$, she communicates her measurement result to Bob via a classical channel.

The state of Bob’s target system $`3`$ conditioned on Alice’s measurement result $`i`$ is given by

$$\rho _3^i=\frac{1}{p_i}\mathrm{tr}_{12}\left[\underset{j}{}(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)(\stackrel{~}{\rho }_1\chi _{23})(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)\right],$$
(3)
where

$$p_i=\mathrm{tr}_{123}\left[\underset{j}{}(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)(\stackrel{~}{\rho }_1\chi _{23})(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)\right].$$

Substituting Eq.(2) into Eq.(3) gives

$$\rho _3^i=\frac{1}{p_i}\mathrm{tr}_{12}\left[\underset{j}{}(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)U_{1(2)3}(\stackrel{~}{\chi }_{12}\rho _3)U_{1(2)3}^{}(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)\right].$$
(4)
Writing

$$\stackrel{~}{\chi }_{12}=\underset{k}{}q_k|\stackrel{~}{s}_k_{12}\stackrel{~}{s}_k|,$$

where the vectors $`|\stackrel{~}{s}_k_{12}`$ make up the complete orthonormal set of eigenvectors of $`\stackrel{~}{\chi }_{12}`$ in the joint space of $`1`$ and $`2`$, and performing the partial trace of Eq.(4) in the complete orthonormal basis $`|\stackrel{~}{P}_l_{12}`$ for the joint system $`12`$ gives

$$\rho _3^i=\frac{1}{p_i}\underset{j,k,l}{}q_k{}_{12}{}^{}\stackrel{~}{P}_l|(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)U_{1(2)3}(|\stackrel{~}{s}_k_{12}\stackrel{~}{s}_k|\rho _3)U_{1(2)3}^{}(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)|\stackrel{~}{P}_l_{12}^{}$$

$$=\underset{j,k,l}{}\left[\sqrt{\frac{q_k}{p_i}}{}_{12}{}^{}\stackrel{~}{P}_l|(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)U_{1(2)3}|\stackrel{~}{s}_k_{12}^{}\right]\rho _3\left[\sqrt{\frac{q_k}{p_i}}{}_{12}{}^{}\stackrel{~}{s}_k|U_{1(2)3}^{}(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)|\stackrel{~}{P}_l_{12}^{}\right].$$

Therefore, $`\rho _3^i`$ is related to $`\rho _3`$ by a quantum operation $`^i`$:

$$\rho _3^i=^i(\rho _3)=\underset{m}{}A_3^{im}\rho _3A_3^{im},$$
(5)
where

$$A_3^{im}\sqrt{\frac{q_k}{p_i}}{}_{12}{}^{}\stackrel{~}{P}_l|(\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)U_{1(2)3}|\stackrel{~}{s}_k_{12}^{},$$

$$\underset{m}{}A_3^{im}A_3^{im}I_3,$$
(6)
and the single index $`m`$ denotes the triple $`(j,k,l)`$.

For Bob to successfully complete the teleportation protocol, he must perform a $`i`$-dependent trace-preserving quantum operation $`^i`$:

$$^i(\rho _3^i)=\underset{n}{}B_3^{in}\rho _3^iB_3^{in},$$

$$\underset{n}{}B_3^{in}B_3^{in}=I_3,$$
(7)
such that the fidelity Jozsa  $`F^i(\rho _3,^i^i(\rho _3))`$ between the input state $`\rho _3`$ and the teleported state $`^i^i(\rho _3)`$ is optimal, that is, as close to one as possible. In other words, Bob has to determine $`^i`$ which optimally reverses $`^i`$.

Before we discuss how optimal $`^i`$ can be determined, we demonstrate that for two-level systems $`1`$, $`2`$ and $`3`$ (from hereon we consider only two-level systems), with

$$\chi _{23}=q_1|\mathrm{\Phi }^+_{23}\mathrm{\Phi }^+|+q_2|\mathrm{\Phi }^{}_{23}\mathrm{\Phi }^{}|+q_3|\mathrm{\Psi }^+_{23}\mathrm{\Psi }^+|+q_4|\mathrm{\Psi }^{}_{23}\mathrm{\Psi }^{}|$$
(8)
where $`0q_k1`$, $`_{k=1}^4q_k=1`$,

$$|\mathrm{\Phi }^\pm _{23}=\frac{1}{\sqrt{2}}(|00_{23}\pm |11_{23}),$$

$$|\mathrm{\Psi }^\pm _{23}=\frac{1}{\sqrt{2}}(|01_{23}\pm |10_{23}),$$

are the Bell states, and

$$\stackrel{~}{\mathrm{\Pi }}_{12}^{1j}=\stackrel{~}{\mathrm{\Pi }}_{12}^1=|\stackrel{~}{\mathrm{\Phi }}^+_{12}\stackrel{~}{\mathrm{\Phi }}^+|,\stackrel{~}{\mathrm{\Pi }}_{12}^{2j}=\stackrel{~}{\mathrm{\Pi }}_{12}^2=|\stackrel{~}{\mathrm{\Phi }}^{}_{12}\stackrel{~}{\mathrm{\Phi }}^{}|,$$

$$\stackrel{~}{\mathrm{\Pi }}_{12}^{3j}=\stackrel{~}{\mathrm{\Pi }}_{12}^3=|\stackrel{~}{\mathrm{\Psi }}^+_{12}\stackrel{~}{\mathrm{\Psi }}^+|,\stackrel{~}{\mathrm{\Pi }}_{12}^{4j}=\stackrel{~}{\mathrm{\Pi }}_{12}^4=|\stackrel{~}{\mathrm{\Psi }}^{}_{12}\stackrel{~}{\mathrm{\Psi }}^{}|,$$
(9)
then $`^1`$ is a generalized depolarizing channel. Here, we use $`|0`$ and $`1`$ to denote an orthonormal set of basis states for each two-level system. Eq.(8) and Eq.(9) allow us to calculate $`p_1=p_2=p_3=p_4=\frac{1}{4}`$. Substituting Eq.(8) into Eq.(2), we obtain

$$\stackrel{~}{\chi }_{12}=q_1|\stackrel{~}{\mathrm{\Phi }}^+_{12}\stackrel{~}{\mathrm{\Phi }}^+|+q_2|\stackrel{~}{\mathrm{\Phi }}^{}_{12}\stackrel{~}{\mathrm{\Phi }}^{}|+q_3|\stackrel{~}{\mathrm{\Psi }}^+_{12}\stackrel{~}{\mathrm{\Psi }}^+|+q_4|\stackrel{~}{\mathrm{\Psi }}^{}_{12}\stackrel{~}{\mathrm{\Psi }}^{}|.$$

That is, the complete orthonormal set of eigenvectors of $`\stackrel{~}{\chi }_{12}`$ in the joint space of $`1`$ and $`2`$ is given by

$$|\stackrel{~}{s}_1_{12}=|\stackrel{~}{\mathrm{\Phi }}^+_{12},|\stackrel{~}{s}_2_{12}=|\stackrel{~}{\mathrm{\Phi }}^{}_{12},|\stackrel{~}{s}_3_{12}=|\stackrel{~}{\mathrm{\Psi }}^+_{12},|\stackrel{~}{s}_4_{12}=|\stackrel{~}{\mathrm{\Psi }}^{}_{12}.$$
(10)
Eq.(6), when $`i=1`$, is in this case reduced to

$$A_3^{1kl}=\sqrt{\frac{q_k}{p_1}}{}_{12}{}^{}\stackrel{~}{P}_l|(\stackrel{~}{\mathrm{\Pi }}_{12}^1I_3)U_{1(2)3}|\stackrel{~}{s}_k_{12}^{}.$$
(11)
Making the choice

$$|\stackrel{~}{P}_1_{12}=|\stackrel{~}{\mathrm{\Phi }}^+_{12},|\stackrel{~}{P}_2_{12}=|\stackrel{~}{\mathrm{\Phi }}^{}_{12},|\stackrel{~}{P}_3_{12}=|\stackrel{~}{\mathrm{\Psi }}^+_{12},|\stackrel{~}{P}_4_{12}=|\stackrel{~}{\mathrm{\Psi }}^{}_{12},$$
(12)
and substituting Eq.(10) and Eq.(12) into Eq.(11) yields a generalized depolarizing channel $`^1`$ specified by Kraus operators

$$A_3^{111}=\sqrt{q_1}I_3,A_3^{121}=\sqrt{q_2}\sigma _3^z,A_3^{131}=\sqrt{q_3}\sigma _3^x,A_3^{141}=\sqrt{q_4}\sigma _3^y,$$
(13)
where

$$\sigma ^x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma ^y=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma ^z=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$

are the Pauli matrices.

Recently, Rehacek et al. Rehacek  found an efficient iterative algorithm for identifying quantum operations $`^i`$ on Bob’s side which optimizes the average teleportation fidelity $`F`$. In particular, for two-level systems $`1`$, $`2`$ and $`3`$, the connection between optimum $`^i`$ and extremal completely positive maps was discussed. Assume that the input states are pure, that is,

$$\stackrel{~}{\rho }_1=|\stackrel{~}{\psi }_1\stackrel{~}{\psi }|,|\stackrel{~}{\psi }_1=\mathrm{cos}\frac{\theta }{2}|0_1+e^{i\varphi }\mathrm{sin}\frac{\theta }{2}|1_1,$$

where $`\theta `$ and $`\varphi `$ are the polar and azimuthal angles respectively. No generality is lost by this assumption, since mixed states are just convex combinations of pure states. $`F`$ is obtained by averaging over all possible Alice’s measurement outcomes $`i`$ and over an isotropic distribution of input states $`\stackrel{~}{\rho }_1`$:

$$F=\underset{i}{}p_i𝑑\psi F^i(\rho _3,^i^i(\rho _3))$$

$$=\frac{1}{4\pi }\underset{i}{}\underset{n}{}\mathrm{tr}_{123}\left[_0^\pi _0^{2\pi }\mathrm{sin}\theta d\theta d\varphi \rho _3B_3^{in}(\underset{j}{}\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)(\stackrel{~}{\rho }_1\chi _{23})B_3^{in}\right]$$

$$=\frac{1}{2}+\frac{1}{12}\underset{i}{}\mathrm{tr}_3(\stackrel{}{X}^i\stackrel{}{T}^i).$$
(14)
Here Rehacek ,

$$\stackrel{}{X}^i\underset{n}{}B_3^{in}\stackrel{}{\sigma }_3B_3^{in}=M\stackrel{}{\sigma }_3+\stackrel{}{c},$$
(15)
where $`M`$ is a $`3\times 3`$ matrix, $`\stackrel{}{c}`$ is a column $`3`$-vector, and

$$\stackrel{}{T}^i\mathrm{tr}_{12}[(\underset{j}{}\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}I_3)(\stackrel{}{\sigma }_1\chi _{23})].$$
(16)
Two observations about Eq.(14) are important. First, quantum operations $`^i`$ corresponding to different Alice’s measurement outcome $`i`$ are independent, therefore each term on the right-hand side of Eq.(14) can be maximized independently. In particular, when we consider the case above, for $`i=1`$, we will be looking for the quantum operation $`^1`$ which optimally reverses the generalized depolarizing channel (13). Second, $`F`$ is a linear functional of $`\stackrel{}{X}^i`$. This means that all its maxima lie on the boundary of the set of physically allowed operators $`\stackrel{}{X}^i`$ that is determined by the constraint of complete positivity of the corresponding transformations. The set of extremal completely positive maps comprising the boundary of the convex set of all completely positive maps thus contains all Bob’s optimum transformations. The matrix $`M`$ can always be brought to a diagonal form via a unitary transformation. When in diagonal form, extremal completely positive maps can be parametrized by two angles $`u[0,2\pi )`$ and $`v[0,\pi )`$, Rehacek

$$M=\left(\begin{array}{ccc}\mathrm{cos}u& 0& 0\\ 0& \mathrm{cos}v& 0\\ 0& 0& \mathrm{cos}u\mathrm{cos}v\end{array}\right),\stackrel{}{c}=\left(\begin{array}{c}0\\ 0\\ \mathrm{sin}u\mathrm{sin}v\end{array}\right).$$

Substituting $`\stackrel{~}{\mathrm{\Pi }}_{12}^1=|\stackrel{~}{\mathrm{\Phi }}^+_{12}\stackrel{~}{\mathrm{\Phi }}^+|`$ and Eq.(8) into Eq.(16), we obtain

$$T_x^1=\frac{1}{4}t_x\sigma _3^x,T_y^1=\frac{1}{4}t_y\sigma _3^y,T_z^1=\frac{1}{4}t_z\sigma _3^z,$$
(17)
where

$$t_x=q_1q_2+q_3q_4,$$

$$t_y=q_1q_2q_3+q_4,$$

$$t_z=q_1+q_2q_3q_4.$$

Note that $`1t_x,t_y,t_z1`$. The fidelity between a pure input state $`\rho _{in}`$ and the output state $`^1^1(\rho _{in})`$, which has gone through the generalized depolarizing channel $`^1`$ followed by a “reversal channel” $`^1`$, averaged over an isotropic distribution of input states is thus given by

$$F^1=𝑑\psi F^1(\rho _{in},^1^1(\rho _{in}))$$

$$=\frac{1}{2}+\frac{1}{6}(t_x\mathrm{cos}u+t_y\mathrm{cos}v+t_z\mathrm{cos}u\mathrm{cos}v).$$
(18)
For $`t_x,t_y,t_z=0`$, we have a totally random channel and the optimum $`F^1`$ is expectedly one-half. If $`0<t_x,t_y,t_z1`$, then $`\mathrm{cos}u=\mathrm{cos}v=1`$ optimizes $`F^1`$. That is, the optimum reversal channel is the “unitary” noiseless channel specified by $`B^1=I`$. The noiseless, bit flip, phase flip, bit-phase flip, and depolarizing channels, with $`q_1>q_2,q_3,q_4`$, belong to this category. Another category is when any two elements from the set $`\{t_x,t_y,t_z\}`$ are negative but the third element is positive, then the optimum reversal channel is the “unitary” channel specified by either $`B^1=\sigma ^x`$, or $`B^1=\sigma ^y`$, or $`B^1=\sigma ^z`$. For instance, $`q_2=1`$ or $`t_x=t_y=1,t_z=1`$, gives the unitary phase flip channel, which can be optimally reversed by $`B^1=\sigma ^z`$. For these two categories, we have

$$F_{\mathrm{max}}^1=\frac{1}{2}+\frac{1}{6}(\mathrm{max}\{4q_11,4q_21,4q_31,4q_41\})$$
(19)
When only one element from the set $`\{t_x,t_y,t_z\}`$ is negative or when $`1t_x,t_y,t_z<0`$, the unitary reversal channels yield Eq.(19). In these cases, for $`|t_x|,|t_y|<|t_z|`$, we could use more general “nonunitary” trace-preserving reversal channels which optimize $`F^1`$. These yield

$$F_{optimal}^1=\frac{1}{2}\frac{1}{6}\frac{t_xt_y}{t_z}$$
(20)
However, it is clear from Eq.(19) that $`F_{optimal}^1<F_{\mathrm{max}}^1`$. Therefore, no general nonunitary trace-preserving reversal channels can do better than unitary ones.

In Rehacek , $`\chi _{23}=|\mathrm{\Psi }^{}_{23}\mathrm{\Psi }^{}|`$ and an imperfect measurement drawn from a one-parametric family of POVMs were considered. We give here the Kraus operators corresponding to the first element in Rehacek ,

$$\stackrel{~}{\mathrm{\Pi }}_{12}^{11}=\frac{\mathrm{sin}\mu }{\sqrt{2}}|11_{12}11|$$

$$\stackrel{~}{\mathrm{\Pi }}_{12}^{12}=\frac{1}{\sqrt{2}\sqrt{1+\mathrm{cos}^2\mu }}(\mathrm{cos}\mu |01_{12}|10_{12})(\mathrm{cos}\mu {}_{12}{}^{}01|{}_{12}{}^{}10|)$$
(21)
In this case, Eq.(6) reduces to

$$A^{114l}=\sqrt{\frac{1}{p_1}}{}_{12}{}^{}\stackrel{~}{P}_l|(\stackrel{~}{\mathrm{\Pi }}_{12}^{11}I_3)U_{1(2)3}|\stackrel{~}{\mathrm{\Psi }}^{}_{12}^{},$$

$$A^{124l}=\sqrt{\frac{1}{p_1}}{}_{12}{}^{}\stackrel{~}{P}_l|(\stackrel{~}{\mathrm{\Pi }}_{12}^{12}I_3)U_{1(2)3}|\stackrel{~}{\mathrm{\Psi }}^{}_{12}^{},$$
(22)
where

$$p_1=\frac{1}{4}(1\mathrm{sin}^2\mu \mathrm{cos}\theta ).$$

Adopting

$$|\stackrel{~}{P}_1_{12}=|00_{12},|\stackrel{~}{P}_2_{12}=\frac{1}{\sqrt{1+\mathrm{cos}^2\mu }}(|01_{12}+\mathrm{cos}\mu |10),$$

$$|\stackrel{~}{P}_3_{12}=\frac{1}{\sqrt{1+\mathrm{cos}^2\mu }}(\mathrm{cos}\mu |01_{12}|10_{12}),|\stackrel{~}{P}_4_{12}=|11_{12},$$
(23)
we have explicitly

$$A^{1144}=\frac{\mathrm{sin}\mu }{\sqrt{1\mathrm{sin}^2\mu \mathrm{cos}\theta }}\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right),A^{1243}=\frac{1}{\sqrt{1\mathrm{sin}^2\mu \mathrm{cos}\theta }}\left(\begin{array}{cc}\mathrm{cos}\mu & 0\\ 0& 1\end{array}\right),$$
(24)
a trace-decreasing quantum channel, which has dependence on the input state. It was shown in Rehacek  that there are values of $`\mu `$ where general nonunitary reversal channels are better than unitary ones.

In conclusion, we have demonstrated how to arrive at the generalized depolarizing channel, Eq.(13), by considering a teleportation scheme involving the entangled state, Eq.(8), and joint measurement, Eq.(9). We then show how the the techniques introduced in Ref. Rehacek  allow us to determine that unitary transformations optimally reverse the generalized depolarizing channel, Eq.(13). Same considerations are given to the teleportation scheme in Ref. Rehacek , and we arrive at a quantum channel which reveals information about the input state. In this case, optimal reversal requires more general quantum operations which are not unitary. In principle, appropriate choices of entangled state $`\chi _{23}`$ and joint measurement $`\stackrel{~}{\mathrm{\Pi }}_{12}^{ij}`$ would yield any desired quantum operation $``$. The iterative procedure in Ref. Rehacek  could then be applied to identify and design the desired optimal reversing quantum operation $``$.

The author thanks Yuri Suhov and Andrew Skeen for useful discussions. This publication is an output from project activity funded by The Cambridge MIT Institute Limited (“CMI”). CMI is funded in part by the United Kingdom Government. The activity was carried out for CMI by the University of Cambridge and Massachusetts Institute of Technology. CMI can accept no responsibility for any information provided or views expressed.


format: text

text:
# Long-range potential fluctuations and 1/𝑓 noise in hydrogenated amorphous silicon

## I introduction

When electric current flows through a resistor, strong low-frequency voltage noise superposed on the thermal Johnson-Nyquist noiseJohnson ; Nyquist  is generally observed. The spectrum of that excess noise has shape close to $`1/f`$, where $`f`$ is frequencyWeissman . In recent decades, the model, according to which the $`1/f`$ noise is produced by an ensemble of two-state systems having broadly distributed activation energies (BDAE)vanderZiel ; duPre ; DDH , has enjoyed a lot of success. However, the BDAE model addresses neither the origin of the two-state systems nor the noise mechanism, thus leaving the absolute noise intensity as an adjustable parameter. The lack of direct microscopic calculations predicting the absolute noise intensity, while being a glaring theoretical gap as such, also fuels an old but still continuing debate over the question of what actually fluctuates — the number of carriers or their mobilityWeissman .

In this work, we present a theoretical and experimental investigation of $`1/f`$ noise in hydrogenated amorphous silicon (a-Si:H). Our theory links the noise to the fluctuations of the number of carriers, predicts the absolute noise intensity, and also allows us to extract novel and useful information about the defects in this important materialStreet .

This theory can be outlined as follows: Hydrogenated amorphous silicon has a significant number of deep defects known as dangling bonds. Thermal fluctuations of the numbers of electrons occupying the dangling bonds cause long-range potential fluctuations, which give rise to fluctuations of the local densities of carriers, which, in turn, lead to resistance fluctuations, which, in the presence of current, manifest themselves as voltage noise. It is the inclusion of the long-range potential fluctuations into the above scheme that distinguishes our treatment from many similar theoretical proposals.

## II Experimental setup

The $`1/f`$ noise in a-Si:H has been studied in the literature in a variety of experimental settings (see e.g. Refs. BA ; BAV ; KK ; Verleg ; GJK ; Johanson ; Goennenwein ; BK ; Kasap ), and exhibited certain features, which depend on numerous details of each particular experiment. In this work we present a fully developed quantitative study of only one situation, which corresponds to our own experiments. Other experimental settings will only be discussed briefly in the end of the paper.

The present study is focused on an $`nin`$ film of a-Si:H, where $`n`$ stands for a 40nm-thick electron doped layer, and $`i`$ for an undoped layer of thickness $`d=0.91\mu `$m. The $`nin`$ structure is grown by plasma enhanced chemical vapor deposition (PECVD) on a highly conductive wafer of crystaline silicon. The contact layer on the top of the structure consists of a 30nm-thick film of titanium followed by 30nm-thick film of copper. The film has area $`A=0.56\text{cm}^2`$. It is thermally annealed and then protected from light. We observe and analyse the voltage noise spectra at frequencies $`f=1÷10^4`$Hz and temperatures $`T=340÷434`$K in the presence of electric current flowing perpendicular to the plane of the film. Other details of our experimental setup are described in Ref. Verleg .

## III Formulation of the theoretical problem

Now we turn to the theoretical derivation of the noise spectrum under the above experimental conditions. The central and, presumably, quite general part of this derivation is Section IV. Most of the rest is specific to a-Si:H and to the experimental setting considered.

Our goal is to compute voltage noise spectrum $`S_V(f)`$ expressed as:

$$\frac{S_V(f)}{V^2}=4_0^{\mathrm{}}C_V(t)\text{cos}(2\pi ft)𝑑t,$$
(1)
where $`V`$ is the applied voltage, and

$$C_V(t)=\frac{\delta R(t)\delta R(0)}{R_0^2}.$$
(2)
Here, $`R_0`$ is the average resistance of the film, and $`\delta R(t)`$ is the equilibrium resistance noise. As usualWeissman , the link between the resistance noise and the voltage noise is established experimentally by observing that $`S_V(f)V^2`$.

The density of states of undoped a-Si:H is shown in Fig. 1a. It is characterized by a band gap of $`1.8`$ eV between the mobility edges $`E_v`$ and $`E_c`$ in the valence and conduction bands, respectively. The proximity of the $`n`$-layers induces band bending in the undoped layer as shown in Fig 1b. Around the center of the film $`E_c(z)=E_{c0}\beta z^2`$, where the $`z`$-axis is directed perpendicular to the film, $`E_{c0}=\mu +0.63`$eV, $`\mu `$ is the chemical potential, and $`\beta =1.6`$eV/$`\mu `$m2. The difference $`E_{c0}\mu `$ was obtained experimentally as the conductivity activation energy, while the value of $`\beta `$ was found numerically (cf.BHD ).

Since $`\mu `$ is significantly closer to $`E_c`$ than to $`E_v`$, the resistivity $`\rho `$ is inversely proportional to the density of conduction electrons $`n_e`$. Keeping in mind that $`n_e\text{exp}\left(\frac{E_c\mu }{k_BT}\right)`$, we obtain $`\rho (z)=X\text{exp}\left(\frac{E_c(z)\mu }{k_BT}\right)`$, which is strongly peaked around $`z=0`$. Here $`X`$ is a proportionality coefficient, and $`k_B`$ is the Boltzmann constant. The total film resistance can then be found as $`R_0=\frac{1}{A}\rho (z)𝑑z\frac{\rho (0)}{A}\sqrt{\frac{\pi k_BT}{\beta }}`$.

To simplify the treatment we replace the actual profile of $`\rho (z)`$ by a layer of constant resistivity $`\rho _0=\rho (0)`$ spreading along the $`z`$-axis between $`z_T`$ and $`z_T`$ (see Fig. 1b). We require the resistance of this constant resistivity layer to be equal to $`R_0`$, which gives

$$z_T=\frac{1}{2}\sqrt{\frac{\pi k_BT}{\beta }}$$
(3)
(typical value: $`0.13\mu `$m).

## IV correlation function

The resistivity within the “constant resistivity layer” still fluctuates as a consequence of the fluctuations of the screened Coulomb potential $`\varphi (t,𝐫)`$ created by deep defects:

$$\varphi (t,𝐫)=\underset{i}{}\frac{\mathrm{\Delta }q_i(t)}{ϵ|𝐫𝐚_i|}\text{exp}\left(\frac{|𝐫𝐚_i|}{r_{s_i}}\right).$$
(4)
Here $`\mathrm{\Delta }q_i(t)`$ is the fluctuation of the $`i`$th defect charge with respect to its average value, $`𝐚_i`$ the position the defect, $`r_{s_i}`$ the screening radius of that defect, and $`ϵ=12`$ the dielectric constant. The defects causing the potential fluctuations may be located outside of the constant resistivity layer. \[Note: the above formula does not include large static component of the random Coulomb potential.\]

When the local potential $`\varphi (t,𝐫)`$ fluctuates, the mobility edge tracks it, i.e. $`E_c(t,𝐫)=E_{c0}+e\varphi (t,𝐫)`$, where $`e`$ is the electron charge. Since $`\mu `$ does not shift with $`e\varphi (t,𝐫)`$, the density of conduction electrons re-equilibrates following $`E_c(t,𝐫)`$ on the timescale of electron drift from the uniform resistivity layer to the $`n`$-layers (and then to the contact layers). Because of the strong band bending, the drift takes less than $`10^7`$s, i.e. the re-equilibration is effectively instantaneous on the timescales of the noise studied ($`\frac{2\pi }{f}10^4÷1`$ s). The resistivity can thus be rewritten as $`\rho (t,𝐫)=X\mathrm{exp}\left(\frac{E_c(t,𝐫)\mu }{k_BT}\right)`$. Assuming for a moment \[and proving later\] that $`|e\varphi (t,𝐫)|k_BT`$ , we expand $`\rho (t,𝐫)=\rho _0+\delta \rho (t,𝐫)`$, where $`\rho _0=X\text{exp}\left(\frac{E_{c0}\mu }{k_BT}\right)`$, and

$$\delta \rho (t,𝐫)=\frac{e\varphi (t,𝐫)}{k_BT}\rho _0.$$
(5)
For $`\delta \rho \rho _0`$, the fluctuation of the total resistance is

$$\delta R(t)=\frac{1}{A^2}_𝒱\delta \rho (t,𝐫)d^3r,$$
(6)
where $`𝒱`$ is the space inside the constant resistivity layer. Substituting $`R_0=\frac{2z_T\rho _0}{A}`$ and $`\delta R(t)`$ given by Eq.(6) into Eq.(2) and then using Eq.(5), we obtain

$$C_V(t)=\left(\frac{e}{2k_BTz_TA}\right)^2_𝒱d^3𝐫_𝒱d^3𝐫^{}\varphi (t,𝐫)\varphi (0,𝐫^{}).$$
(7)
## V Description of Defects

In order to evaluate $`\varphi (t,𝐫)\varphi (0,𝐫^{})`$, we have to describe the deep defects in the undoped layer of a-Si:Hnumbers .

There exist both theoretical arguments and experimental evidenceSJT  indicating that the concentration of defects in the undoped layer should be greater in the vicinity of the $`n`$-layers and then decay towards the center on the lengthscale of about $`0.5\mu `$m. Here, however, in order to simplify the theoretical treatment, we assume that the defect concentration across the entire undoped layer has a constant value, which we estimate asnumbers  $`n_D=6\times 10^{15}\text{cm}^3`$.

Each defect has 4 possible states: one zero-electron state $`D^+`$, two one-electron states $`D^0`$, and one two-electron state, $`D^{}`$. We assume a Gaussian probability distribution

$$p(E)=\frac{1}{\sqrt{2\pi }\mathrm{\Delta }E}\text{e}^{\frac{(EE_0)^2}{2\mathrm{\Delta }E^2}}$$
(8)
for energy $`E`$ of an electron occupying the $`D^0`$ state. Here $`E_0=\mu 0.22`$eV and $`\mathrm{\Delta }E=0.15`$eV. Energy $`E`$, as such, is associated only with the $`D^+D^0`$ transition. The $`D^0D^{}`$ transition, which requires capturing the second electron by the same defect, is characterized by energy $`E^{}=E+U`$, where $`U=0.2`$eV is the correlation energy. Therefore, the corresponding probability distribution is

$$p^{}(E^{})=\frac{1}{\sqrt{2\pi }\mathrm{\Delta }E}\text{e}^{\frac{(E^{}E_0U)^2}{2\mathrm{\Delta }E^2}}.$$
(9)
We shall treat the $`D^+D^0`$ transitions independently from the $`D^0D^{}`$ transitions, which is justified as long as exp$`\left(\frac{U}{k_BT}\right)1`$.

The noise spectrum obtained later will depend only weakly on defect parameters $`n_D`$, $`E_0`$, $`U`$ and $`\mathrm{\Delta }E`$, mainly through a weak dependence on one combination of them — the density of defect states at the chemical potential. That combination, in turn, is not very sensitive to the choice of $`E_0`$ and $`U`$.

Using the equilibrium descriptionAY  of the $`D^+`$, $`D^0`$ and $`D^{}`$ states, we obtain the mean squared charge fluctuations for the $`D^+D^0`$ transition:

$$\mathrm{\Delta }q^2(E)_+=\frac{1}{4}e^2\text{sech}^2\left(\frac{E\mu }{2k_BT}c\right);$$
(10)
where $`c=\frac{1}{2}\text{ln}2`$. For the $`D^0D^{}`$ transition:

$$\mathrm{\Delta }q^2(E^{})_{}=\frac{1}{4}e^2\text{sech}^2\left(\frac{E^{}\mu }{2k_BT}+c\right).$$
(11)
In order to escape from a deep defect, an electron should reach the mobility edge $`E_c`$. However, the activation barriers $`E_B`$ (indicated in Fig.1c) can vary as a result of the medium-range disorder of the amorphous structure (on a lengthscale of $`1÷10`$nm). We assume a Gaussian probability distribution for the values of $`E_B`$:

$$P(E_B)=\frac{1}{\sqrt{2\pi }\mathrm{\Delta }E_B}\mathrm{exp}\left(\frac{(E_BE_{B0})^2}{2\mathrm{\Delta }E_B^2}\right).$$
(12)
where $`E_{B0}`$ and $`\mathrm{\Delta }E_B`$ are to be extracted from the experimental spectra. These are the only two adjustable parameters in our treatment. They will affect the spectral shape, but not the integrated noise intensity.

The fluctuation rate for the $`D^+D^0`$ transition isRBZ

$$\frac{1}{\tau _+(E,E_B)}=\omega _0(E_B)\left[\text{e}^{\frac{E_B\mu }{k_BT}}+\frac{1}{2}\text{e}^{\frac{E_BE}{k_BT}}\right];$$
(13)
and for the $`D^0D^{}`$ transition:

$$\frac{1}{\tau _{}(E^{},E_B)}=\omega _0(E_B)\left[\frac{1}{2}\text{e}^{\frac{E_B\mu }{k_BT}}+\text{e}^{\frac{E_BE^{}}{k_BT}}\right],$$
(14)
where the attempt frequency is estimated from the detailed balance condition as

$$\omega _0(E_B)=v_{th}\sigma 𝒩_c(E_B)$$
(15)
Here, $`v_{th}=\sqrt{3k_BT/m_e^{}}`$ is the thermal velocity of conduction electrons; $`m_e^{}`$ their effective mass ($`m_e^{}0.4m_e=3.6410^{28}`$g); $`\sigma =10^{15}\text{cm}^2`$ the cross-section of electron capture by $`D^+`$ or $`D^0`$ defect; and $`𝒩_c(E_B)N_c(E_B)k_BT`$ the “concentration” of thermally accessible electronic states with energies above $`E_B`$; $`N_c(E_B)=N_c(E_{c0})\sqrt{(E_BE_{c0}+\epsilon _c)/\epsilon _c}`$ the empirical fit to the density of states above the mobility edge; $`N_c(E_{c0})=4\times 10^{21}\text{eV}^1\text{cm}^3`$; and $`\epsilon _c=0.02\text{eV}`$. The typical value of $`\omega _0(E_B)`$ is then $`10^{13}\text{s}^1`$.

Next we evaluate the screening radius. Since the screening by the conduction electrons in the undoped layer can be neglected in view of their very small concentration ($`10^{10}÷10^{13}`$cm-3), two other screening mechanisms should be considered, namely: (i) by the $`n`$-layers together with the contact layers; and (ii) by the defects in the undoped layer.

We describe the first mechanism as a perfect screening by metallic surfaces. In other words, we assume that a defect in the undoped layer is screened by the infinite set of its mirror images constructed with respect to the mirror planes located at $`z=\pm d/2`$. Although the screening law due to this mechanism is not exponential, we obtain the best possible value for $`r_s`$ entering Eq.(4) as the distance from the defect to the point, where the screened potential is factor of $`\text{e}=2.718\mathrm{}`$ smaller than the bare Coulomb potential. We have found numerically that $`r_s`$ has angular dependence, which is such that, for a defect located at $`z=0`$, $`r_{sx}=r_{sy}=0.54d`$, and $`r_{sz}=0.38d`$. We then approximate the screening radius due to this mechanism as

$$r_{s1}=(r_{sx}r_{sy}r_{sz})^{1/3}=0.48d.$$
(16)
The second screening mechanism is similar to that of Thomas-Fermi. It can be described as follows: When the potential fluctuates due to charge fluctuation on a given defect, the energies $`E`$ and $`E^{}`$ of the neighboring defects become shifted. In response, those neighboring defects change their occupation numbers thus screening the potential of the defect, which originally caused the fluctuation. The equilibrium description of this mechanism SE  results in the screening radius

$$r_{_{\text{eq}}}=\sqrt{\frac{ϵ}{4\pi e^2n_D\nu (T)}},$$
(17)
where

$$\nu (T)=\frac{1}{4k_BT}\left\{_{\mathrm{}}^{\mathrm{}}\stackrel{~}{p}(E)𝑑E+_{\mathrm{}}^{\mathrm{}}\stackrel{~}{p}^{}(E^{})𝑑E^{}\right\};$$
(18)
$$\stackrel{~}{p}(E)=p(E)\text{sech}^2\left(\frac{E\mu }{2k_BT}c\right),$$
(19)
and

$$\stackrel{~}{p}^{}(E^{})=p^{}(E^{})\text{sech}^2\left(\frac{E^{}\mu }{2k_BT}+c\right).$$
(20)
The value of $`\nu (T)`$ depends on $`T`$ very weakly. The product $`n_D\nu (0)`$ should be recognized as the density of the defect states at the chemical potential.

The equilibrium description of screening is applicable, when the sources of the Coulomb potential are static, and thus all the defects can contribute to screening. However, if a given defect fluctuates on time scale $`\tau `$, its potential can only be screened by other defects fluctuating on timescales not slower than $`\tau `$. To take this observation into account, we simply multiply $`n_D`$ in formula (17) by $`b(\tau )`$, the fraction of defects that fluctuate fast enough to take part in the screening of a defect characterized by the fluctuation time $`\tau `$. We estimate that fraction as

$`b(\tau )=`$ $`{\displaystyle \frac{1}{4k_BT\nu (T)}}\{{\displaystyle _{_+(\tau )}}\stackrel{~}{p}(E)P(E_B)dEdE_B`$ (21)
$`+{\displaystyle _{_{}(\tau )}}\stackrel{~}{p}^{}(E^{})P(E_B)dE^{}dE_B\},`$
where $`_+(\tau )`$ is the integration region limited by condition $`\tau _+(E,E_B)<\tau `$, and $`_{}(\tau )`$ the integration region limited by $`\tau _{}(E^{},E_B)<\tau `$. The screening radius due to the second mechanism thus becomes

$$r_{s2}(\tau )=\sqrt{\frac{ϵ}{4\pi e^2n_D\nu (T)b(\tau )}}.$$
(22)
For the slower fraction of defects, $`r_{s2}(\tau )`$ approaches $`r_{_{\text{eq}}}0.18\mu `$m, which is shorter than $`r_{s1}=0.43\mu `$m. At the same time, for the faster ones, $`r_{s2}`$ is very large, i.e. they are predominantly screened by the first mechanism. We, therefore, approximate the combined effect of the two mechanisms by using the following expression for the screening radius:

$$r_s^1(\tau )=r_{s1}^1+r_{s2}^1(\tau ).$$
(23)
The typical value of $`r_s`$ given by Eq.(23) is $`r_s^{}=0.2\mu `$m.

## VI evaluation of the noise spectrum

Now we are in position to evaluate $`\varphi (t,𝐫)\varphi (0,𝐫^{})`$ with $`\varphi (t,𝐫)`$ given by Eq.(4).

Since an electron emitted by one deep defect is most likely absorbed not by another defect but by the contact layer, it is appropriate to assume that different defects fluctuate independently, which implies that, for $`ij`$, $`\mathrm{\Delta }q_i(t)\mathrm{\Delta }q_j(0)=0`$. Keeping also in mind that

$$\mathrm{\Delta }q_i(t)\mathrm{\Delta }q_i(0)=\mathrm{\Delta }q_i^2\mathrm{exp}(t/\tau _i),$$
(24)
we take the disorder average of $`\varphi (t,𝐫)\varphi (0,𝐫^{})`$, i.e. we replace the discrete summation by the integrations over the relevant probability distributions. This gives

$`\varphi (t,𝐫)\varphi (0,𝐫^{})={\displaystyle \frac{n_D}{ϵ^2}}{\displaystyle 𝑑E_BP(E_B)}`$ $`\{{\displaystyle }dEp(E)F(𝐫,𝐫^{},\tau _+(E,E_B))\mathrm{\Delta }q(E)^2_+\mathrm{exp}({\displaystyle \frac{t}{\tau _+(E,E_B)}})`$ (25)
$`+{\displaystyle }dE^{}p^{}(E^{})F(𝐫,𝐫^{},\tau _{}(E^{},E_B))\mathrm{\Delta }q(E^{})^2_{}\mathrm{exp}({\displaystyle \frac{t}{\tau _{}(E^{},E_B)}})\},`$
where

$`F(𝐫,𝐫^{},\tau )`$ $`{\displaystyle \frac{\text{exp}\left(\frac{|𝐫𝐚|+|𝐫^{}𝐚|}{r_s(\tau )}\right)}{|𝐫𝐚||𝐫^{}𝐚|}d^3𝐚}`$ (26)
$`2\pi r_s(\tau )\mathrm{exp}\left({\displaystyle \frac{|𝐫𝐫^{}|}{r_s(\tau )}}\right).`$
The approximation in Eq.(26) can be justified by observing that the main spatial dependence of the integral for $`F(𝐫,𝐫^{},\tau )`$ has form $`\mathrm{exp}\left(\frac{|𝐫𝐫^{}|}{r_s(\tau )}\right)`$, and then the natural choice for the prefactor is $`F(0,0,\tau )=2\pi r_s(\tau )`$ (obtained by direct integration).

From Eq.(25), the value of $`|e\varphi (t,𝐫)|`$ can be estimated as: $`e\sqrt{\varphi ^2(0,0)}e^2ϵ^1\sqrt{2\pi n_D\nu (T)k_BTr_s^{}}3.5\text{meV}`$. Since $`k_BT30`$meV, the assumption $`|e\varphi (t,𝐫)|k_BT`$ made earlier was adequate.

Finally, we evaluate $`C_V(t)`$ by substituting Eq.(25) into Eq.(7), and then take the Fourier transform (1) to obtain

$`{\displaystyle \frac{S_V(f)}{V^2}}={\displaystyle \frac{2\pi ^2e^4n_D}{ϵ^2(k_BT)^2Az_T^2}}`$
$`\times {\displaystyle }dE_BP(E_B)\{{\displaystyle }D(f,\tau _+(E,E_B))\stackrel{~}{p}(E)dE`$
$`+{\displaystyle }D(f,\tau _{}(E^{},E_B))\stackrel{~}{p}^{}(E^{})dE^{}\},`$ (27)
where

$$D(f,\tau )=\frac{\tau r_s^4(\tau )\left[4z_T3r_s(\tau )+\text{e}^{\frac{2z_T}{r_s(\tau )}}(2z_T+3r_s(\tau ))\right]}{1+4\pi ^2f^2\tau ^2}.$$
(28)
Recalling the approximations associated with (i) the assumption of the constant resistivity layer; (ii) the assumption of the constant profile of defect concentration; (iii) the treatment of the screening radius; and (iv) the evaluation of $`F(𝐫,𝐫^{},\tau )`$ (Eq.(26)), we estimate that the integrated noise intensity obtained from Eq.(27) entails a factor of two theoretical uncertainty for the noise mechanism considered.

The complex appearance of Eq.(27) is due to the fact that it includes separate terms for the two transitions $`D^+D^0`$ and $`D^0D^{}`$. Such a separation is necessary only because the expressions (13, 14) for the activation times $`\tau _+`$ and $`\tau _{}`$ are slightly different from each other. This difference, however, vanishes if one assumes that the width of the distributions $`\stackrel{~}{p}(E)`$ and $`\stackrel{~}{p}^{}(E^{})`$ is much smaller than the width of $`P(E_B)`$, i.e. $`2k_BT\mathrm{\Delta }E_B`$.

The complicated form of the numerator in the right-hand side of Eq.(28) is the analytic result of the two integrations (7) over the volume of the “constant resistivity” layer. This expression reflects the competition between two length scales: the half-thickness of that layer ($`z_T`$) and the screening radius ($`r_s`$). In the “thin layer limit” $`z_Tr_s`$, the expression in the square brackets in Eq.(28) can be approximated by $`\frac{2z_T^2}{r_s}`$, whereas in the “bulk limit” $`z_Tr_s`$ that expression approaches $`4z_T`$.

Now we simplify Eq.(27) by taking limits $`z_Tr_s`$, $`2k_BT\mathrm{\Delta }E_B`$ and substituting $`\nu (T)=\frac{1}{2\mathrm{\Delta }E}`$ (just an approximate numerical fact), which gives

$$\frac{S_V(f)}{V^2}=\frac{8\pi ^2e^4n_D}{ϵ^2k_BT\mathrm{\Delta }EA}\frac{\tau _0(E_B)r_s^3\left(\tau _0(E_B)\right)P(E_B)dE_B}{1+4\pi ^2f^2\tau _0^2(E_B)},$$
(29)
where

$`{\displaystyle \frac{1}{\tau _0(E_B)}}`$ $``$ $`{\displaystyle \frac{1}{\tau _+(\mu ,E_B))}}{\displaystyle \frac{1}{\tau _{}(\mu ,E_B))}}`$ (30)
$`=`$ $`1.5\omega _0(E_B)\mathrm{exp}\left({\displaystyle \frac{E_B\mu }{k_BT}}\right).`$
Our actual system is characterized only by the weaker inequalities $`z_T<r_s`$, $`2k_BT<\mathrm{\Delta }E_B`$. Therefore, for comparison with experiments we shall still use the original formula (27). At the same time, for the qualitative analysis, we will focus on Eq.(29), but all the conclusions will be fully applicable to Eq.(27).

A remarkable feature of Eq.(29), which can be traced back to Eq.(7), is that, even though the resistance noise is caused by the fluctuations in the number of conduction electrons, the resulting noise spectrum is independent of their equilibrium concentration. Furthermore, that spectrum is only weakly dependent on all the defect parameters involved (see Ref. numbers ). In particular, the dependence of the prefactor on the concentration of defects $`n_D`$ is balanced by the $`n_D^{1/2}`$ dependence of $`r_{s2}`$, which then enters $`r_s^3`$ via Eq.(23).

Like in the BDAE model, the $`1/f`$-like spectral shape generated by formula (29) is the result of the broad distribution of the activation energies $`P(E_B)`$, but, at the same time, formula (29) has also two new features, namely: (i) the $`1/T`$ dependence of the prefactor and (ii) the energy-dependent weight $`r_s^3\left(\tau (E_B)\right)`$ multiplying $`P(E_B)`$. The second feature is particularly important for extracting the correct distribution of $`P(E_B)`$ from experimental data.

For comparison with experiment in Section VII we will need the integral of $`\frac{S_V(f)}{V^2}`$ over all frequencies, which is equal to $`C_V(0)`$. The expression for $`C_V(0)`$ can be obtained by substituting $`\frac{1}{4}r_s^4(\tau )\left[4z_T3r_s(\tau )+\text{e}^{\frac{2z_T}{r_s(\tau )}}(2z_T+3r_s(\tau ))\right]`$ instead of $`D(f,\tau )`$ into Eq.(27). The mathematical structure of this expression is similar to Eq.(27) but otherwise not very illuminating to be written explicitly one more time. Instead, we give an estimate for the noise integral corresponding to the thin film limit (29). We do it with one further simplification. Namely, we replace $`r_s(\tau (E_B))`$ in Eq.(29) by the typical value $`r_s^{}=0.2\mu `$m and then obtain:

$$_0^{\mathrm{}}\frac{S_V(f)}{V^2}𝑑f=\frac{2\pi ^2e^4n_Dr_{s}^{}{}_{}{}^{3}}{ϵ^2k_BT\mathrm{\Delta }EA}.$$
(31)
Although the bulk limit $`z_Tr_s`$ appears to be incompatible with the assumption of independent defect fluctuations (to be explained in Section VIII), it is still instructive to give an expression for the noise integral in this case too. Combining the limit $`z_Tr_s`$ with all other approximations used to derive Eq.(31), we obtain

$$_0^{\mathrm{}}\frac{S_V(f)}{V^2}𝑑f=\frac{4\pi ^2e^4n_Dr_{s}^{}{}_{}{}^{4}}{ϵ^2k_BT\mathrm{\Delta }EAz_T}.$$
(32)
Up to numerical prefactors, both approximations (31) and (32) can be summarized as follows: the noise intensity is proportional to $`\left(\frac{e}{k_BT}\right)^2`$, multiplied by the mean squared amplitude of potential fluctuations $`\frac{e^2n_Dk_BTr_s^{}}{ϵ^2\mathrm{\Delta }E}`$, further multiplied by the screening volume of a typical fluctuation $`r_{s}^{}{}_{}{}^{3}`$, and, finally, divided by the volume of the space, where the defects contributing to the potential fluctuations are located. In the thin layer limit that volume is roughly $`2r_sA`$, whereas in the bulk limit it is $`2z_TA`$. The combination $`n_Dk_BT/\mathrm{\Delta }E`$ appearing in the estimate of the potential fluctuation should be identified with the concentration of “thermally active” defects, i.e. those defects that fall in the thermal energy window around the chemical potential.

Now we give various estimates of the prefactor in front of the approximate $`1/f`$ dependence of the spectrum. For the thin film limit, we start from Eq.(29) and then, assuming $`r_s(\tau (E_B))=r_s^{}`$, $`P(E_B)=\frac{1}{2\mathrm{\Delta }E_B}`$ and $`\omega _0(E_B)=\omega _0(E_{B0})`$, obtain

$$\frac{S_V(f)}{V^2}=\frac{\pi ^2e^4n_Dr_{s}^{}{}_{}{}^{3}}{ϵ^2A\mathrm{\Delta }E\mathrm{\Delta }E_B}\frac{1}{f}.$$
(33)
For the bulk limit, the analogous expression is

$$\frac{S_V(f)}{V^2}=\frac{2\pi ^2e^4n_Dr_{s}^{}{}_{}{}^{4}}{ϵ^2Az_T\mathrm{\Delta }E\mathrm{\Delta }E_B}\frac{1}{f}.$$
(34)
Equation(34) admits yet another remarkable simplification, if one substitutes $`r_{\text{eq}}`$ given by Eq.(17) instead of $`r_s^{}`$. (The equilibrium self-screening mechanism described by Eq.(17) should indeed be a proper description for the slower fraction of the defects, which means the lower frequency part of the spectrum.) In this case

$$\frac{S_V(f)}{V^2}=\frac{\mathrm{\Delta }E}{2Az_Tn_D\mathrm{\Delta }E_B}\frac{1}{f}\frac{1}{N_Df},$$
(35)
where $`N_D=2Az_Tn_D`$ is the total number of defects in the constant resistivity layer. \[In the above approximation, we assumed $`\mathrm{\Delta }E\mathrm{\Delta }E_B`$.\] Although this result looks like a standard statistical factor, its origin is, in fact, far more complex. It can be traced back to an accidental interplay between the resistivity fluctuations and the screening mechanism. Given the limitations of the “bulk limit”, it is unlikely that formula (35) represents a clean limit in any realistic system. At the same time, for $`r_sz_T`$, this formula still gives a reasonable estimate.

It is, finally, interesting to observe that, for the numbers characterizing our system, the Hooge’s formulaHooge  $`\frac{S_V(f)}{V^2}=\frac{\alpha }{N_Cf}`$ with coefficient $`\alpha 10^3`$, would produce an estimate close to that of Eq.(35), (Here $`N_C`$ is the total number of conduction electrons.) This simply reflects the fact, that in our case the concentration of conduction electrons is about thousand times smaller than the concentration of defects. The exact ratio of these concentrations is strongly temperature dependent, as is the value of $`\alpha `$ required to fit the noise intensity in our experiments.

## VII Comparison with experiment

Theoretical spectra computed from Eq.(27) are compared with the experimental ones in Fig. 2a (see also Ref. supplement ). The experimental spectra were obtained by subtracting the zero-current noise from the total noise observed with $`V=50`$ meV. They were consistent with the spectra reported in Ref. Verleg , though the film was newly grown. All the theoretical spectra were obtained with $`P(E_B)`$ shown in Fig. 2b and characterized by $`E_{B0}=\mu +0.90`$eV and $`\mathrm{\Delta }E_B=0.09`$eV. The uncertainties of the fit for $`E_{B0}`$ and $`\mathrm{\Delta }E_B`$ are $`0.05`$ eV and $`0.02`$ eV, respectively. Figure 2c shows $`r_s(E_B)`$ obtained by averaging $`r_s\left(\tau _+(E,E_B)\right)`$ and $`r_s\left(\tau _{}(E^{},E_B)\right)`$ over $`E`$ and $`E^{}`$ at $`T=434`$K. Finally, Fig. 2d illustrates the role of the $`r_s^3`$ weight by showing $`P(E_B)r_s^3(E_B)`$. If $`P(E_B)`$ were to be obtained by fitting the experimental spectra to the BDAE model Verleg , the result would look like Fig. 2d, i.e. the maximum would be shifted by about 0.1 eV in comparison with Fig. 2b.

With the above value of $`\mathrm{\Delta }E_B`$, the estimates (33,34,35) for the prefactor in front of $`1/f`$ would give respectively: $`910^{12}`$, $`310^{11}`$ and $`110^{11}`$ — all in reasonable agreement with experiments. \[For the estimates (34) and (35) we used $`z_T=0.13\mu \text{m}`$.\]

In Fig. 3 we present another test of our theory, which is independent of the choice of $`E_{B0}`$ and $`\mathrm{\Delta }E_B`$. Namely, we compare the theoretical and experimental values of the integrated noise intensities for the four temperatures indicated in Fig. 2a.

The experimental evaluation of the integrated intensity of an $`1/f`$-like noise is a task notorious for its ambiguity. Fortunately, in our case, the power law extrapolations of all four noise spectra were convergent, which allowed us to use the following procedure:

First, we obtained the lower ends of the error bars by integrating the experimental noise spectra only in the frequency range of the actual experimental observations. Then, the upper ends were obtained by making power law extrapolations of the spectra beyond the frequency range of observation (up to $`10^6`$ Hz for small frequencies and $`10^8`$ Hz for large frequencies), and then adding the integrals over the extrapolated tails to the lower end values of the error bars. Finally, the “experimental” points indicated in Fig. 3 were chosen as the middle points of the above error bars.

The theoretical points presented in Fig. 3 were obtained for the full spectrum (27) as described in Section VI. One can also check that the estimate (31) gives numbers, which are not much different, e.g.: $`5.610^{11}`$ for $`T=340`$ K, and $`4.410^{11}`$ for $`T=434`$ K.

Given the large uncertainties of the extrapolation and the theoretical uncertainties indicated earlier, the factor of two agreement between the theoretical and experimental points in Fig. 3 should be considered as a successful consistency test of the noise mechanism proposed.

## VIII Possible generalizations

Now we discuss briefly various modifications required in order to generalize our treatment to $`1/f`$ noise experiments in other settings involving a-Si:H.

The application of the present theory to the films with coplanar currentsKK ; GJK ; Johanson ; BK ; Kasap  and also to the bulk samples would encounter an essential complication related to the absence of contact layers. In this case, the charge fluctuations of different defects become correlated via emission and capture of the same electron. In such a process, the potential fluctuation is simply translated in space from one defect to another, which means that, unless the two defects have different screening radii, the total resistance does not fluctuate at all. Our analysis in Section V revealed one possibility for different defects to have different screening radii (as a function of their activation times). However, even with this possibility, the overall effect of the absence of contact layers should be a noticeable reduction of the noise intensity.

The contact layers also play a role in the screening of defects. Without them (and without the $`n`$-layers) the first screening mechanism considered in Section V is not operational. Therefore, the treatment of screening of the faster fluctuating defects becomes more difficult.

The theory has to be further modified to include the hole conduction, if, as a result of doping, the chemical potential moves in the middle or below the middle of the band gap. An interesting situation may arise, when hole resistivity is equal to electron resistivity. In this case one should expect a drastic reduction of the noise intensity, because the change of electron resistivity induced by potential fluctuations will be compensated by the change in the hole resistivity.

Analysing experiments, one should also keep in mind that creation and annealing of charged defects represents an alternative way to potential fluctuations. In particular, in a-Si:H such a process may be associated with the diffusion of hydrogen atomsStreet , which can inhibit or expose the dangling bonds. The parameters describing this process are not very well known. In our case, the mechanism involving the emission and the capture of conduction electrons by already existing defects appears to be sufficient, but, in general, creation and annealing of defects can lead to a comparable contribution to the integrated noise intensity.

Summarizing this Section, we would like to emphasise that all the above complications can still be treated on the basis of formula (7), but different ingredients for the evaluation of the potential fluctuations will be required.

## IX Conclusions

In conclusion, we have developed a microscopic theory of $`1/f`$ noise in $`nin`$ sandwich structures of a-Si:H and found a very good agreement between this theory and our experiments. The noise mechanism proposed should be quite general (in particular, in semiconductors): it merely requires the presence of charge traps in regions of poor screening. The full calculation of the actual noise spectrum necessarily involves a fair amount of material specific details, which, to some extent, mask the generality of our treatment. However, one should remember, that our calculation is based on formula (7), which appears to have quite a broad range of applicability.

From a different perspective, our analysis introduces a new method of characterizing the deep defects in a-Si:H. In particular, our finding, that $`E_{B0}E_{c0}=0.27`$eV, which is substantially greater than 0.1 eV — the scale of potential fluctuations seen in the drift mobility experiments Street , may indicate that the defects are not distributed randomly in the background of the medium-range random potential but, instead, located in the regions of greater local strain, i.e. around the peaks of that potential.

Our result also limits the spectral intensities due to all other possible noise mechanisms to the difference between the experimental spectra and the spectra given by Eq.(27). Even with the freedom of varying $`E_{B0}`$ and $`\mathrm{\Delta }E_B`$, this is a very strong constraint.

###### Acknowledgements.

We thank R. E. I. Schropp for providing us with the sample, and A. Buchleitner, B. Farid and M. Weissman for their helpful comments on the manuscript. The work of B. V. F. was supported by the Foundation of Fundamental Research on Matter (FOM), which is sponsored by the Netherlands Organization for the Advancement of Pure Research (NWO).


text:
# 1 Introduction

## 1 Introduction

Most contemporary explanations of Gamma Ray Bursts (GRBs) attribute them to high entropy “fireballs” which convert their energy into a matter-dominated (baryonic) jet within which relativistic electrons and electromagnetic field are subsequently re-created (e.g. Piran 1999). In our opinion, there are a number of problems with this scenario. Besides the commonly recognized problems of low efficiency and baryon mass and magnetic-field fine-tuning, the creation of the plasma-dominated flows is problematic. It is usually envisaged that the energy is released electromagnetically (e.g., MacFadyen & Woosley 1999; Kim et al. 2002; Wheeler et al. 2000) and is dissipated right away into “lepto-photonic” plasma. Yet, there is no clear mechanism for transforming magnetic energy into a hot plasma on an outflow time-scale. Also, the baryonic jets are inferred to have Mach numbers $`300`$ and to maintain a ratio of bulk to internal energy $`>10^5`$, (or equivalently to maintain a flow that deviates locally by angles of less than $`10^{}`$), while dissipating heavily through internal shocks.

We explore an alternative model for GRBs in which the energy remains in electromagnetic form all the way from the its origin to the sites of $`\gamma `$-ray and afterglow emission.

### 1.1 The case for electromagnetic energy transport

Several arguments point to the importance of electromagnetic field in producing GRBs: (i) Magnetic fields naturally collimate outflows and have been invoked in all other jet outflows. (ii) Magnetic outflows are “clean” – instead of requiring a finely tuned baryon fraction, baryons are may be absent completely. (iii) There is no need to convert Poynting flux into fluid energy flux and then back again. (iv) Electromagnetic energy is “high quality” – it is in a low entropy form and can be efficiently dissipated through current instabilities during particle acceleration. (v) Electromagnetic dissipation is intrinsically intermittent (c.f. solar flares); the observed GRB variability need not be tied to the source and can consequently arise at much larger radii than in internal shock models, thereby relaxing the $`\gamma `$-ray opacity constraint. (vi) “Standard candle” – simple electromagnetic models offer an explanation for the surprisingly narrow reported burst energy distribution – approximately the same burst energy will be inferred independent of the viewing angle (c.f. Frail et al. 2001; Panaitescu et al. 2002; Lazzati et al. 2002). (vii) Ultrarelativistic electromagnetic outflows are formally subsonic which implies that they are more naturally established than hypersonic, fluid jets. (viii) There are astrophysical examples of collimated Poynting flux dominated outflows (pulsars, AGNs, black hole candidates).

### 1.2 Sources of Electromagnetic Outflow

We assume that the GRB source (sometimes referred to as “millisecond magnetar”, Usov 1992) is a strongly magnetized neutron star (e.g., Kouveliotou et al., 1998) or a stellar mass black hole-accretion torus with mass $`M110M_{}`$, size given by the effective light cylinder radius $`r_010`$ km, angular velocity $`\mathrm{\Omega }r_0/c10^4`$ rad s-1 and magnetic field strength $`B_010^{15}`$ G. Such super-strong magnetic fields may be generated by $`\alpha \omega `$ dynamos (e.g., Thompson & Murray, 2001), or just by differential rotation (e.g., Kluźniak & Ruderman, 1998). A rotational energy, $`E0.1Mr_0^2\mathrm{\Omega }^210^{52}`$ erg – a one parameter quantity, is available to power GRB bursts (intermediate mass black holes can supply more energy.) The field is chosen to produce a characteristic electromagnetic power $`L\mathrm{\Delta }\mathrm{\Omega }L_\mathrm{\Omega }B_0^2r_0^6\mathrm{\Omega }^4/c^310^{50}`$ erg s-1, comfortably larger than that inferred for GRBs and giving a source lifetime $`100`$ s matched to long bursts.

### 1.3 The Electrical Circuit

Independent of the source, we suppose that the combination of magnetic flux $`\mathrm{\Phi }B_0r_0^210^{27}`$ G cm2 and angular velocity combine to create a unipolar inductor with EMF $`0.1\mathrm{\Omega }\mathrm{\Phi }/c10^{22}`$ V and that the axisymmetric part of the electromagnetic field dominates at large radius. An equivalent and useful way to think about this is to say that there is a strong, quadrupolar current distribution outward along the axes and inward along the equator (or vice versa). Under electromagnetic conditions, the effective impedance is roughly that of free space – $``$ 100 $`\mathrm{\Omega }`$ is fine for estimates – so the current flowing around the circuit is $`I10^{20}`$ A and the power that is delivered to the “load”‘is $`LI10^{49}`$ erg s-1. (This interpretation of the load resistance is valid even if there is no dissipation and the electromagnetic Poynting flux just propagates away from the source with little reflection.) Our proposal differs from the conventional interpretation principally through the assumption that the current flows all the way out to the expanding blast wave, rather than completes close to the source (c.f. Fig. 2). In addition, we make the conjecture that the complex magnetic field geometry, that must be present within $`r_0`$, sorts itself out and becomes axisymmetric and primarily toroidal for $`r>>r_0`$, by analogy to what is observed in the quiet solar wind. This can occur completely without (e.g. Blandford 2002) or with partial (e.g. Lyutikov 2002b) dissipation of magnetic field.

### 1.4 Lepton Loading

The sources that are envisaged generally release energy though baryonic, leptonic and electromagnetic channels. The first is usually parameterized by the initial baryon rest mass fraction of the total energy density $`\eta `$, which we shall assume to be quite small. It is also useful to define a quantity $`\sigma `$ which is the ratio of the electromagnetic to the total matter energy density.

Somewhat paradoxically, it gets harder to convert electromagnetic energy directly to pair plasma the stronger the field becomes. The reason is that plasma tries to short out any electric field along magnetic field on a very short time scale, a few Langmuir periods. Even if there is not enough plasma to enforce $`𝐄𝐁=0`$ conditions the required amount of pairs will be created through vacuum breakdown. Hoverer, pair creation is not likely to drain all the potential EMF since the potential differences required to create pairs through a vacuum breakdown are typically $`1`$ GV and never more than $`1`$ TV, which is orders of magnitude smaller than the EMF required to account for ultra-relativistic outflows $`10`$ ZV in GRBs. When the field strength is as high as $`10^{15}`$ G, the minimum density of plasma needed to short out the electric field – the Goldreich-Julian density $`n_{\mathrm{GJ}}`$– is tiny in comparison with the equipartition density $`n_{\mathrm{eq}}`$. Put another way, $`\sigma _0<\sigma _{\mathrm{max}}n_{\mathrm{eq}}/n_{\mathrm{GJ}}/\mathrm{\Delta }V_{\mathrm{vac}}\omega _\mathrm{G}/\mathrm{\Omega }10^{16}`$.

The actual value of $`\sigma _0`$ is model-dependent. At one extreme, if there is an evacuated neutron star or black hole magnetosphere, only enough pairs need be created to supply the necessary space charge and current so that $`\sigma _0\sigma _{\mathrm{max}}`$. At the other, there might be a weakly magnetized hot torus radiating neutrinos with energy above $`5`$ MeV so that $`\sigma <<1`$ and the standard fireball model would apply. All intermediate values are possible (e.g., dissipation of large scale magnetic field through pair creation is likely to give $`\sigma _0/\mathrm{\Delta }V_{breakdown}10^{10}`$). If we suppose that the flow is initially electromagnetically dominated, with $`\sigma 1`$, then about the only way that it seems possible to create entropy directly and reduce $`\sigma `$ appreciably is for an electromagnetic turbulence cascade to develop, operating down to wavelengths, $`\lambda _{\mathrm{min}}`$, so small that electromagnetic energy can be dissipated in accelerating pairs, somewhat analogous to the viscous dissipation that terminates a fluid turbulence spectrum. If this really can operate then it is hard to see how more than a few percent of the electromagnetic energy will be dissipated on a source dynamical time scale, more specifically, $`\sigma \mathrm{ln}(r_0/\lambda _{\mathrm{min}})100`$, and the flow will remain electromagnetically dominated.

In fact, even if there is an appreciable pair content close to the source (but still $`\sigma 1`$), this will quickly become irrelevant because the pairs will annihilate when their temperature falls to $`20`$ keV, just like in the early universe. The Thomson optical depth of the plasma will then drop to a low enough value that the photons can escape and decouple from the flow. If this happens in GRBs then a weak $`\gamma `$-ray precursor is predicted (Lyutikov & Usov 2000; Fig. 1) and this may have been observed (Preece 2002).

### 1.5 Outflow Phases

Before discussing a few of the details, we should distinguish the principal stages through which the outflow passes in our GRB model.

1. Flow formation ($`r<r_0`$) There must be a quasi-steady source of electromagnetic power within $`r_010`$ km. For a long burst, this should last for $`10^6`$ dynamical times. The initial temperature of the pair plasma is $`T_03(\sigma _0/100)^{1/4}`$ MeV and the associated optical depth $`\tau _\mathrm{T}(B/B_\mathrm{L})(n/n_{\mathrm{GJ}})`$, where $`B_\mathrm{L}5\times 10^{15}`$ G is the Larmor field for which the Larmor radius of a mildly relativistic electron equals the classical electron radius. The flow is initially sufficiently optically thick to trap the radiation.

2. “Warm” acceleration ($`r_0<r<r_{\mathrm{thin}}`$). The plasma is accelerated by the electromagnetic field and the bulk Lorentz factor increases linearly with radius until the flow becomes optically thin at $`r_{\mathrm{thin}}10^8(\sigma _0/100)^{1/4}`$. A small admixture of baryons will extend $`r_{\mathrm{thin}}`$ but it cannot delay the escape of photons by much.

3. Electromagnetic bubble. Beyond the photosphere the radiation and $`e^\pm `$ decouple from magnetic field - loading drops by 7 orders of magnitude, leaving plasma strongly magnetized $`\sigma 10^9`$. At this moment the flow becomes a relativistically expanding magnetically dominated bubble. Several stages of the bubble expansion can be identified:

3.1. Magnetic acceleration ($`r_{\mathrm{thin}}<r<ct_\mathrm{s}`$). The radial Lorentz factor of the frame in which the electric vector vanishes can reach $`\mathrm{\Gamma }10^4`$, as long as it does not exceed $`\eta ^1`$. A filled electromagnetic bubble will be produced with polar and axial current that complete as a “Chapman-Ferraro” current along a relativistically expanding contact discontinuity (CD) that separates the bubble from the swept up circumstellar medium. The magnetic field will be primarily toroidal and the electric field should be primarily poloidal in the frame of the explosion. This phase will continue until $`ct_\mathrm{s}3\times 10^{12}`$ cm, when the source switches off.

3.2. “Coasting” electromagnetic shell ($`ct_\mathrm{s}<r<r_{\mathrm{sh}}(L_\mathrm{\Omega }t_\mathrm{s}^2/\rho c)^{1/4}`$). At this time, the bubble becomes a relativistically expanding shell of thickness $`ct_\mathrm{s}`$. The shell still contains toroidal magnetic field but the current now detaches from the source and completes along the shell’s inner surface. At this stage the CD is constantly re-energized by the fast-magnetosonic waves propagating from the central source so that the Lorentz factor of the CD is $`\mathrm{\Gamma }(L_\mathrm{\Omega }/\rho c^3)^{1/4}r^{1/2}`$ (in a constant density medium) or $`\mathrm{\Gamma }\mathrm{const}`$ (in a $`\rho r^2`$ wind). This stage is limited by one dynamical time scale $`t_{dyn}2t_\mathrm{s}\mathrm{\Gamma }^2`$.

3.3. Self-similar electromagnetic shell ($`r_{\mathrm{sh}}<r<r_{\mathrm{NR}}(L_\mathrm{\Omega }t_\mathrm{s}/\rho c^2)^{1/3}`$). After one dynamical time scale all the regions of the bubble come into a causal contact – most of the waves reflected from the CD have propagated throughout the bubble. As the expanding bubble performs a work on the surrounding medium its total energy decreases; the amount of energy that remains in the bubble at the self-similar stage needs to be calculated numerically. Most energy in the bubble is still concentrated in a thin shell with $`\mathrm{\Delta }RR/\mathrm{\Gamma }^2`$ near the surface of the bubble which is moving according to $`\mathrm{\Gamma }\sqrt{E_\mathrm{\Omega }/\rho c^2}r^{3/2}`$ (in a constant density medium), or $`\mathrm{\Gamma }r^{1/2}`$ (in a $`\rho r^2`$ wind). Interestingly, the structure of the magnetic bubble resembles at this stage the structure of the hydrodynamical relativistic blast wave wave (Blandford&McKee 1976). After the bubble came into a causal contact it starts to evolve in lateral direction trying to adjust magnetic pressures. This evolution may be accompanied by electro-magnetic instabilities which lead to particle acceleration and $`\gamma `$-rays production. The result of the lateral energy redistribution is a creation of an anisotropic expansion, being faster at a given time In addition, $`\gamma `$-rays may be produced throughout the CD surface due to development of current instabilities or inertial acceleration (Smolsky & Usov 1996).

4. Relativistic blast wave ($`r_{\mathrm{sh}}<r<r_{\mathrm{NR}}(L_\mathrm{\Omega }t_\mathrm{s}/\rho c^2)^{1/3}`$). As the bubble expands its energy is gradually transfered to the preceding forward shock wave. Most efficiently this transfer occurs at the end of the coasting phase (at the coasting phase $`dE/dtr`$). During the self-similar phase the energy in the magnetic bubble decreases slowly (logarithmically in times since at this stage $`dE/dt1/r`$), so that the relativistic blast wave stage is coexistent with the self-similar stage of the magnetic bubble. At the forward shock wave relativistic particles are accelerated producing the observed afterglow. Magnetic flux is either incorporated into the blast wave from the bubble or, perhaps, amplified at the relativistic shock front. This is the afterglow phase which usually becomes unobservably faint after the expansion speed becomes mildly relativistic at $`r10^{18}`$ cm.

5. Non-relativistic blast wave ($`r_{\mathrm{NR}}<r`$). Eventually the source will expand non-relativistically and become more spherical with time, resembling a normal supernova remnant.

## 2 Dynamics of magnetic explosions

### 2.1 Relativistic MHD vs force-free formalism

The conventional method for handling relativistic, magnetized flows is to use the relativistic extension of regular, non-relativistic magnetohydrodynamics (RMHD). However, there is a simpler extension which is appropriate when the plasma is sufficiently tenuous that its inertia can be ignored, though sufficiently dense that it can supply the space charge and current density ($`1<<\sigma <<10^{16}`$). Under these circumstances, we adopt the relativistic force-free (RFF) approximation

$$\rho 𝐄+𝐣\times 𝐁=0$$
(1)
This implies that $`𝐄𝐁`$ and its temporal derivative can be set to zero. In addition, we restrict our attention to the case $`E<B`$. Eq. (1) allows us to define an electromagnetic velocity $`𝐯=𝐄\times 𝐁/B^2`$ perpendicular to the magnetic field (it is not useful to define a component along the magnetic field).

Now, Maxwell’s equations can be written

$$\frac{𝐄}{t}=\times 𝐁4\pi 𝐣,\frac{𝐁}{t}=\times 𝐄,$$
(2)
where we can use $`(𝐄𝐁/t=0`$ to derive

$$𝐣=\left[(𝐄\times 𝐁)𝐄+(𝐁\times 𝐁𝐄\times 𝐄)𝐁\right]/B^2.$$
(3)
This RFF equation set represents a simple, evolutionary dynamical system (Uchida 1997). When one includes the constraints $`𝐄𝐁=𝐁=0`$, there are four, independent electromagnetic variables to evolve and four characteristics along which information is propagated. In the linear approximation, these correspond to forward and backward propagating fast and intermediate wave modes with phase speeds $`c`$ and $`c\widehat{𝐤}\widehat{𝐁}`$ respectively.

RFF dynamics can be developed in a manner that is quite analogous to regular hydrodynamics, with the anisotropic Maxwell stress tensor taking the place of the regular pressure and the electromagnetic energy density playing the role of inertia. There is an important difference, though, in that the existence of a luminal fast mode means that electromagnetic “flows” do not become truly “supersonic”.

The RMHD approach, which must be used in the presence of significant inertia is based upon the fluid velocity $`𝐮`$. The electric field is generally supposed to be related to the magnetic field through an infinite conductivity Ohm’s law, $`𝐄+𝐮\times 𝐁=0`$ and there are seven independent variables to evolve, with seven independent characteristics (two fast, two intermediate, two slow and one adiabatic) to track.

### 2.2 Suppression of lateral dynamics

Relativistic flows do not collimate easily on account of relativistic, kinematic effects. A simple way to see this is by analyzing the MHD force balance equation in the $`\theta `$ direction

$$_t[(w+b^2)\gamma ^2\beta _\theta ]+\frac{1}{r\mathrm{sin}\theta }_\theta [\mathrm{sin}\theta (p+b^2/2)]\mathrm{cot}\theta \frac{pb^2}{r}=0$$
(4)
where $`w`$ is the enthalpy and $`b=(B^2E^2)^{1/2}`$ is the magnetic field in the comoving frame. Eq. (4) shows that typically on a flow expansion time $`\beta _\theta 1/\gamma ^2\theta `$ \- the lateral dynamics is frozen-out for ultrarelativistic flow.

In the context of our model, the flow pattern is established in the vicinity of the source. In the subsequent motion, the flow at different latitudes drops out of causal contact until $`rr_{\mathrm{dyn}}`$, when connection is re-established and lateral motion becomes possible (the situation is analogous to the “Hello–Goodbye–Hello” kinematics familiar from the theory of inflation). As a consequence, steady state solutions based on the Grad-Shafranov equation are not likely to become valid on an expansion time-scale.

### 2.3 Radial motion

As a consequence of freezing of lateral dynamics we can separate out the radial motion, except, perhaps, close to the axis. To find the simplest RFF solution during the expanding shell phase, we suppose that the only non-zero field components are $`B_\varphi ,E_\theta `$. The system (3) can be solved by separation of variables:

$$B_\varphi =r^1[f_1(tr)+f_2(t+r)]g(\theta ),E_\theta =r^1[f_1(tr)f_2(t+r)]g(\theta )$$
(5)
We now impose two boundary conditions at the contact discontinuity, $`r=R`$, specifically pressure balance and velocity matching

$$\frac{B^2}{8\pi \mathrm{\Gamma }^2}=2\rho _{\mathrm{ext}}c^2\mathrm{\Gamma }(R)^2,E_\theta =0B(R)^2E_\theta (R)^2=B(R)^2/\mathrm{\Gamma }(R,\theta )^2$$
(6)
During the self-similar wave phase, the waves emitted by the source have all caught up with the CD – the magnetic bubble relaxes to a self-similar (in $`tr`$ coordinates) structure. Assuming that in the self-similar regime $`\mathrm{\Gamma }^2=\mathrm{\Gamma }_0^2t^m`$ (the coasting phase with the constant energy source is equivalent to self-similar solution with power supply, $`m=1`$) we find

$`f_1(tr)=\sqrt{16\pi \rho _{\mathrm{ext}}c^2}\mathrm{\Gamma }^2t^{2m/(1+m)}\left(2(m+1)\mathrm{\Gamma }^2\left(tr\right)\right)^{(1m)/(1+m)}`$
$`f_2(t+r)=\sqrt{16\pi \rho _{\mathrm{ext}}c^2}(t+r)/8`$ (7)
The reflected (in-going wave), $`f_2`$, is $`\mathrm{\Gamma }^2`$ times weaker than the outgoing - Doppler red-shifting during the reflection from the CD. Magnetic field:

$$B\frac{f_1(tr)}{r}g(\theta )=\sqrt{16\pi \rho _{\mathrm{ext}}c^2}\mathrm{\Gamma }^2\frac{t}{r}\chi ^{(1m)/(1+m)}g(\theta )$$
(8)
$`\chi =2(m+1)\mathrm{\Gamma }^2(1r/t)`$. The self-similar structure of the magnetic bubble looks similar to the structure of hydrodynamic blast wave (Blandford&McKee 1974)! This result is a bit surprising since the two systems are completely different (electro-magnetic bubble and hydrodynamical shock wave) and solutions come from different equations (Maxwell and Euler). An important property of the self-similar bubble is that for $`1<m<3`$ the magnetic field is concentrated near CD in a thin sheath $`R/\mathrm{\Gamma }^2`$.

In particular, for point explosion in a homogeneous medium $`\rho _{\mathrm{ext}}=`$const ($`m=3`$) we find

$`B{\displaystyle \frac{1}{t^2r\mathrm{sin}\theta \sqrt{\chi }}},E_\mathrm{\Omega }{\displaystyle \frac{E_0}{\mathrm{sin}^2\theta }}\mathrm{ln}{\displaystyle \frac{E_0}{\rho _{\mathrm{ex}}c^2r^3}}`$
$`\mathrm{\Gamma }^2={\displaystyle \frac{1}{\mathrm{sin}^2\theta }}\left({\displaystyle \frac{r_{nr}}{r}}\right)^3,r_{nr}=\left({\displaystyle \frac{3E_0}{2\rho _{\mathrm{ex}}c^2}}\right)^{1/3}10^{18}\mathrm{cm}`$ (9)
Interestingly, the energy contained within the magnetic bubble decreases logarithmically  with time, so that the forward blast wave fully decouples from the magnetic bubble only when the flow becomes trans-relativistic.

### 2.4  Lateral dynamics and collimation.

At large distances, $`r10^{16}`$ cm, the flow slows down and the lateral dynamics “un-freezes”. The flow then tries to adjust to a lateral force balance, which becomes $`g=1/\mathrm{sin}\theta `$, $`\mathrm{\Gamma }1/\mathrm{sin}^2\theta `$, $`L_\mathrm{\Omega }1/\mathrm{sin}^2\theta `$ when the toroidal field dominates over pressure. Thus, at $`r10^{16}cm`$ the flow relaxes to a universal lateral distribution of energy independent of the initial conditions.The energy flow is strongly peaked along the axis.

## 3 $`\gamma `$-ray emission

In the RFF limit the fast speed is the speed of light and so no fast shocks form. If we add a limited quantity of plasma and use RMHD, the shocks are weak and not likely to be efficient particle accelerators. There are no intermediate shocks in the RFF limit, though rotational discontinuities can be present.

We therefore propose that the $`\gamma `$-ray-emitting electrons are accelerated by current instabilities during the magnetic shell phase. Intense line and sheet currents are traditionally unstable as is well-documented by laboratory experiments. There are two possible locations of the emission, along the poles and in the body of the magnetic shell. We consider these in turn.

### 3.1 Polar Emission

We can continue the cosmological analogy and note that as the flow decelerates, angular scales $`\mathrm{\Gamma }^1`$ will “enter the horizon”. Specifically, electromagnetic instabilities on the axis can grow on progressively larger scales. The fastest growing modes are likely to be pinch, kink and helical modes. with progressively larger scale modes becoming unstable as the flow develops, also analogous to cosmology.

There are several ways in which these instabilities can develop. Let us outline one scenario that we find to be quite plausible. If we continue to make the RFF approximation and ignore the inertia and pressure of the plasma outflow then the solution outlined above will have an unbalanced stress near the axis where the current flows. This will cause the electromagnetic field to flow towards the axis. This inflow will be accompanied by development of current instabilities – pinches and kinks. We suspect that the nonlinear development of these instabilities will produce a turbulence spectrum down to a high $`k`$ inner scale where the particle acceleration is able to absorb the wave energy on the wave turnover time-scale. Simple estimates suggest that the slope of the wave turbulence spectrum is $`3/2`$ as opposed to the Kolmogorovian value of $`5/3`$ and that the energy of the emitted radiation can, indeed, lie in the $`\gamma `$-ray band, (after Doppler boosting). Numerical computations of this nonlinear development are underway (MacFadyen & Blandford, in preparation).

One merit of this explanation is that it allows the GRB to originate at a much larger radius, up to $`10^{16}`$ cm than in the standard, baryonic, intermediate shock model $`10^{13}`$ cm. This means that the bulk outflow Lorentz factor can be much smaller at the emission generation radius. (One objection that has been raised to this type of model is that the it may not be possible to achieve the large variability on $`10`$ ms time-scales that is observed. However, this turns out not to be a problem if the emission zone contains large amplitude waves moving with speed $`c`$, as is likely to occur in an intermittent turbulence spectrum.)

### 3.2 Shell Emission

The expanding shell will be preceded by a surface current emanating from the poles and flowing to the equator (or vice versa). It will be followed, at a distance of the order of the shell thickness by a reverse current. We strongly suspect (though have not demonstrated) that this electromagnetic configuration is strongly unstable and, after a few light crossing times in the frame of the source, the shell will break up and create a local turbulence spectrum which can accelerate particles as described above. This will occur at $`r10^{16}`$ cm. The variation, particle acceleration and emission are much as in the polar model. Alternatively, acceleration of external electrons due to inertia-induced electric fields at the contact discontinuity (Smolsky & Usov 1996) may lead to the production of $`\gamma `$-rays.

An attractive feature of the shell emission is that it predicts $`\gamma `$-ray emission in all directions and few orphan afterglows. The strength of the burst is expected to be related to the inferred energy per sterad in the afterglow although large variations are indeed possible. Burst seen at small angle (as may be inferred from achromatic breaks in the afterglow emission) should be seen to larger distances although a large burst to burst dispersion should be expected. It could be that short bursts are polar and long bursts are from the shell. Simulations of shell emission are also underway which will clarify some of these points.

## 4 Afterglow Emission

Under this model, the afterglow is produced during the relativistic blast wave phase when the energy contained in the electromagnetic bubble is transferred to the swept up circumstellar medium. However, if the electromagnetic dynamics is as conjectured above, then this will be reflected in the subsequent blast wave dynamics. The form of the blast wave would reflect both the form of the driver (magnetic bubble) and subsequent evolution of the shock. We have considered the expansion of non-spherical relativistic blast waves in a relativistic Kompaneets approximation (Kompaneets 1960, Shapiro 1980, Lyutikov 2002b). We find that only extremely strongly collimated shocks, with the opening angle $`\mathrm{\Delta }\theta 1/\mathrm{\Gamma }^2`$ show modification of profiles due to sideways expansion. Thus, the motion of the forward shock will be determined by the form of the electro-magnetic driver. In particular, the energy per sterad in the shock is $`L_\mathrm{\Omega }1/\theta ^2`$. What is inferred as a jet is a non-spherical relativistic outflow. This implies that, whatever the observer orientation, the inferred blast wave energetics is roughly the same, to within a logarithmic factor. The starting Lorentz factor for the afterglow is likely to be $`\mathrm{\Gamma }>>\theta ^1`$. The emission should exhibit an achromatic break when $`\mathrm{\Gamma }`$ falls to $`\theta ^1`$. Thus, all bursts are observable, although the fluency depends strongly upon angle of observation. This has important implications for the incidence of orphan afterglows and perhaps for the interpretation of X-ray flashes.

## 5 Conclusion

We have explored the “electromagnetic hypothesis” for ultra-relativistic outflows, namely that they are essentially electromagnetic phenomena which are driven by energy released by spinning black holes or neutron stars and that this electromagnetic behavior continues into the source region even when the flows become non-relativistic. The most striking implications of the electromagnetic hypothesis are that particle acceleration in the sources is due to electromagnetic turbulence rather than shocks and that the outflows are cold and electromagnetically-dominated, with very few baryons, at least until they become strongly dissipative.

The most important prediction that is likely to be tested in the coming years is the form of the external shock wave generated by the magnetic bubble. We predict $`L_\mathrm{\Omega }\theta ^2`$. Observational properties of such blast wave have been already discussed (Rossi at al. 2002). Current data are consistent with such an energy distribution. Further data are expected from optical polarization observations, which should exhibit a characteristic temporary evolution as well as correlation with intensity and the times of the jet achromatic breaks (Rossi at al. 2002).

## Acknowledgements

We are indebted to Andrew MacFadyen, Martin Rees and Elena Rossi for advice and encouragement. Support under NASA grant NAG 5-2837, 5-7007 is gratefully acknowledged.

Fig. 1. Temporal structure of GRBs (from Lyutikov & Usov 2000)

Fig. 2. Magnetic fields and currents in the bubble

## Questions

J. Trier Frederiksen: What are the assumptions of the medium that the magnetic bubble is pushing?

M. Lyutikov: There is no specific requirement on the external medium (except that it is tenuous): it can be hot or cold, homogeneous or a power-law in radius. The internal structure of the bubble will be somewhat different in these cases, but will remain qualitatively similar.

S. Shibata: What is the physics leading to ”delayed reconnection” in the distance far from the central engine?

M. Lyutikov: Relativistic radial motion leads to an effective freezing of lateral motion and of current instabilities, which becomes important only when the flow slows down.

E. Berger: What are the relevant length scales ($`10^{16}`$ cm for GRB, $`10^{17}`$ for afterglow, $`10^{18}`$ cm for non-relativistic phase) and what determines these scales?

M. Lyutikov: $`\gamma `$-rays are emitted at $`ct\mathrm{\Gamma }^2ct_s(Et_s/\rho c)^{1/4}10^{16}`$ cm at the end of the coasting phase when lateral dynamics of the bubble un-freezes. Starting from this radius, most of the energy is transferred to the forward shock which becomes non-relativistic at $`r(E/\rho c^2)^{1/3}10^{18}`$ cm.