text:
# Entanglement cost of mixed states

## Abstract

We compute the entanglement cost of several families of bipartite mixed states, including arbitrary mixtures of two Bell states. This is achieved by developing a technique that allows us to ascertain the additivity of the entanglement of formation for any state supported on specific subspaces. As a side result, the proof of the irreversibility in asymptotic local manipulations of entanglement is extended to two-qubit systems.

Developing a theory of entanglement is considered a priority in the field of quantum information, where quantum correlations are a precious resource for information processing Ni99 . In particular, the quest for proper entanglement measures has received much attention in recent years Ho01 . From the identification and study of properties of such measures a gain of insight into the nature of entanglement is expected. In turn, their computation for particular states provide us with an account of the resources present in those states.

Two measures of entanglement stand out due to their physical meaning. Both of them refer to the possibility of transforming entangled states of a bipartite system by means of local operations and classical communication (LOCC). The distillable entanglement Be96 ; Ra99  $`E_d(\rho )`$ quantifies how much pure-state entanglement can be extracted from $`\rho `$. More specifically, it gives the ratio $`M/N`$ in the large $`N`$ limit, where $`M`$ is the number of ebits \[i.e. entangled bits, or maximally entangled states $`(|00+|11)/\sqrt{2}`$ of a two-qubit system\] that can be distilled from the state $`\rho ^N`$ using LOCC. The entanglement cost Be96 ; Ha00  $`E_c(\rho )`$ quantifies, instead, the amount of pure–state entanglement needed to create $`\rho `$. It is defined in the limit of large $`N`$ as the ratio $`M/N`$, where $`M`$ is the number of ebits required to prepare $`\rho ^N`$ using LOCC.

The outputs produced so far by entanglement theory concerning these two entanglement measures include, among others, the following remarkable results:

($`i`$) All forms of bipartite pure-state entanglement are equivalent in the asymptotic limit Be96b , in the sense that for large $`N`$ and any bipartite pure state $`|\mathrm{\Psi }`$, $`|\mathrm{\Psi }^N`$ can be reversibly converted into ebits. Thus, for pure states $`E_d(|\mathrm{\Psi })=E_c(|\mathrm{\Psi })`$, with the so–called entropy of entanglement $`E(|\mathrm{\Psi })`$ denoting the resulting unique measure.

($`ii`$) Two forms of bipartite entanglement, namely free and bound entanglement Ho98 , have been identified for mixed states. The first form corresponds to mixed states that can be distilled, i.e. $`E_d>0`$. Bound entangled states were defined as those that cannot be distilled into pure-state entanglement, i.e. $`E_d=0`$, in spite of the fact that they cannot be produced \[in the non-asymptotic regime\] by just mixing product (i.e. unentangled) pure states.

($`iii`$) Contrary to the pure–state case, the asymptotic manipulation of some entangled mixed states is irreversible Vi01 . This follows from the gap observed between the distillable entanglement and the entanglement cost, $`E_d<E_c`$, for some mixed states. This phenomenon occurs both for bound entangled states and for distillable states.

Notice that, as far as mixed states are concerned, the above results are qualitative footnote2 . In particular, the entanglement cost $`E_c`$ has not been computed for any mixed state. This problem is related to the one of the additivity of the entanglement of formation $`E_f(\rho )`$ Be96 ; Wo01 , an auxiliary measure that quantifies how much pure–state entanglement —as given by $`E`$— is required to create a single copy of the mixed state $`\rho `$. In particular, it is not known whether $`E_f(\rho ^N)=NE_f(\rho )`$, which would imply that $`E_c=E_f`$.

In this paper we compute the value of the entanglement cost $`E_c`$ for all mixed states $`\rho _V`$ supported on some specific subspaces $`VH_AH_B`$. This is achieved by showing that the entanglement of formation $`E_f`$ is additive for the tensor product $`\rho _V\sigma `$,

$$E_f(\rho _V\sigma )=E_f(\rho _V)+E_f(\sigma ),$$
(1)
where $`\sigma `$ is an arbitrary bipartite state, which by iteration implies $`E_c(\rho _V)=E_f(\rho _V)`$. We also present a technique that allows us to evaluate $`E_f`$ for some classes of mixed states.

Our considerations include, in a two-qubit system, a mixture $`\rho _p`$ of two Bell states, $`|\mathrm{\Phi }^\pm (|00\pm |11)/\sqrt{2}`$,

$$\rho _p(1p)|\mathrm{\Phi }^+\mathrm{\Phi }^+|+p|\mathrm{\Phi }^{}\mathrm{\Phi }^{}|,p[0,\frac{1}{2}],$$
(2)
for which we obtain

$$E_c(\rho _p)=H_2(\frac{1}{2}+\sqrt{p(1p)}),$$
(3)
$`H_2(x)`$ being the Shannon entropy $`S(x,1x)`$. The distillable entanglement of $`\rho _p`$ reads distillable

$$E_d(\rho _p)=1H_2(p),$$
(4)
and thus $`E_d(\rho _p)<E_c(\rho _p)`$ for all $`p(0,1/2)`$. That is, even the process of preparing the elementary mixture $`\rho _p`$ is irreversible, in that not all the pure–state entanglement employed can be subsequently recovered by asymptotic LOCC. This constitutes a new, remarkably simple instance of the irreversibility that takes place in the asymptotic manipulation of entanglement.

Mathematically, the entanglement of formation of a mixed state $`\rho (H_AH_B)`$ can be expressed as Be96 ; Wo98

$$E_f(\rho )=\underset{dD_\rho }{inf}\underset{k}{}p_kE(|\psi _k),$$
(5)
where the entropy of entanglement $`E(|\psi )`$ is given by $`S(\rho _A)`$, $`\rho _A\mathrm{tr}_B|\psi \psi |`$, and the minimization is performed over the set $`D_\rho `$ of all pure–state realizations $`d\{p_k,|\psi _k\}`$ of $`\rho `$, $`\rho =_kp_k|\psi _k\psi _k|`$. The entanglement cost, in turn, corresponds to Be96 ; Ha00

$$E_c(\rho )=\underset{N\mathrm{}}{lim}\frac{E_f(\rho ^N)}{N}.$$
(6)
Our first goal is to show the additivity of $`E_f`$ for some mixed states $`\rho _V`$, as expressed in Eq. (1), which implies that $`E_c(\rho _V)=E_f(\rho _V)`$. For concreteness, we start by discussing a simple example. After that a theorem announces the result in its full generality.

Example 1. Let us consider two qubits $`A`$ and $`B`$, with Hilbert spaces $`H_A=H_B=𝒞^2`$, and a subspace $`VH_AH_B`$ spanned by the vectors noteTensor

$`|0_V`$ $``$ $`|0_A|0_B,`$
$`|1_V`$ $``$ $`|1_A|1_B.`$ (7)
Notice that, in particular, the mixture $`\rho _p`$ of Eq. (2) is supported on $`V`$, $`\rho _p(V)`$. For any vector $`|\varphi _VV`$, $`|\varphi _V=c_0|0_V+c_1|1_V`$, we can define a vector $`|\varphi _AH_A`$ as $`|\varphi _Ac_0|0_A+c_1|1_A`$. Then the operation of tracing out qubit $`B`$

$$|\varphi _V\varphi |tr_B(|\varphi _V\varphi |)=\underset{\alpha =0,1}{}|c_\alpha |^2|\alpha _A\alpha |,$$
(8)
which is a trace-preserving, completely positive (TPCP) map from $`(V)`$ to $`(H_A)`$, can also be described as a TPCP map $``$ from $`(H_A)`$ to $`(H_A)`$,

$$tr_B(|\varphi _V\varphi |)=(|\varphi _A\varphi |),$$
(9)
given by

$$(X)=\underset{\alpha =0,1}{}\mathrm{tr}(|\alpha _A\alpha |X)|\alpha _A\alpha |.$$
(10)
The relevant feature of subspace $`V`$ is that $``$ is an entanglement–breaking map Sh02 , as Eq. (10) makes manifest entbreak . The theorem below establishes that this property alone guarantees that Eq. (1) holds for any $`\rho _V(V)`$, and thus $`E_c(\rho _V)=E_f(\rho _V)`$. Then we can use the closed expression for $`E_f`$ in two-qubit systems Wo98  to evaluate $`E_f(\rho _V)`$, and thereby obtain, for instance, the value of $`E_c(\rho _p)`$ displayed in Eq. (3).

More generally, the theorem considers four arbitrary quantum systems, denoted $`A`$, $`B`$, $`a`$ and $`b`$, and refers to the entanglement between $`Aa`$ and $`Bb`$. $`V`$ is a subspace of $`H_{AB}`$ such that, for any $`|\mathrm{\Psi }_{Vab}VH_{ab}`$, tracing out subsystem $`B`$ destroys all existing entanglement between $`AB`$ with $`ab`$,

$$\mathrm{tr}_B(|\mathrm{\Psi }_{Vab}\mathrm{\Psi }|)=\underset{l}{}q_l|\mu _l_A\mu _l||\nu _l_{ab}\nu _l|,$$
(11)
that is, such that the map $`(V)(H_A)`$ given by $`\rho _V\mathrm{tr}_B\rho _V`$ is entanglement–breaking Sh02 .

Theorem: Let $`\rho _V(V)`$ and $`\sigma _{ab}(H_{ab})`$. Then

$$E_f(\rho _V\sigma _{ab})=E_f(\rho _V)+E_f(\sigma _{ab}).$$
(12)
Proof: Notice that $`E_f(\rho _V\sigma _{ab})E_f(\rho _V)+E_f(\sigma _{ab})`$, because from optimal pure–state decompositions of $`\rho _V`$ and of $`\sigma _{ab}`$ a (possibly non-optimal) decomposition for $`\rho _V\sigma _{ab}`$ can be constructed with average $`E`$ given by $`E_f(\rho _V)+E_f(\sigma _{ab})`$. In what follows we will show that

$$E_f(\rho _V\sigma _{ab})E_f(\rho _V)+E_f(\sigma _{ab}).$$
(13)
Let us consider a decomposition $`\{p_k,|\mathrm{\Psi }_k\}`$ of $`\rho _V\sigma _{ab}`$,

$$\rho _V\sigma _{ab}=\underset{k}{}p_k|\mathrm{\Psi }_k\mathrm{\Psi }_k|,$$
(14)
such that it is optimal, that is,

$$E_f(\rho _V\sigma _{ab})=\underset{k}{}p_kE(|\mathrm{\Psi }_k).$$
(15)
We recall that all $`|\mathrm{\Psi }_k`$ must belong to $`VH_{ab}`$. Next we argue that in order to prove Eq. (13) —and therefore the theorem— it is sufficient to show that for any pure state $`|\mathrm{\Psi }_{Vab}VH_{ab}`$,

$$E(|\mathrm{\Psi })E_f(\varrho _V)+E_f(\pi _{ab}),$$
(16)
where

$`\varrho _V`$ $``$ $`\mathrm{tr}_{ab}(|\mathrm{\Psi }_{Vab}\mathrm{\Psi }|),`$ (17)
$`\pi _{ab}`$ $``$ $`\mathrm{tr}_{AB}(|\mathrm{\Psi }_{Vab}\mathrm{\Psi }|).`$ (18)
Indeed, denoting by $`\varrho _V^k`$ and $`\pi _{ab}^k`$ the reduced density matrices of systems $`AB`$ and $`ab`$ for each $`|\mathrm{\Psi }_k`$, we would have

$`E_f(\rho _V\sigma _{ab})`$ $`=`$ $`{\displaystyle \underset{k}{}}p_kE(|\mathrm{\Psi }_k)`$ (19)
$``$ $`{\displaystyle \underset{k}{}}p_kE_f(\varrho _V^k)+{\displaystyle \underset{k}{}}p_kE_f(\pi _{ab}^k)`$
$``$ $`E_f({\displaystyle \underset{k}{}}p_k\varrho _V^k)+E_f({\displaystyle \underset{k}{}}p_k\pi _{ab}^k)`$
$`=`$ $`E_f(\rho _V)+E_f(\sigma _{ab}).`$
where the first inequality assumes Eq. (16), the second inequality uses that $`E_f`$ is a convex function and the last step follows from the fact that $`\rho _V=_kp_k\varrho _V^k`$ and $`\sigma _{ab}=_kp_k\pi _{ab}^k`$. Let us then move to justify Eq. (16). $`E(|\mathrm{\Psi }_{Vab})`$ is given by the von Neumann entropy of the reduced density matrix $`\xi _{Aa}\mathrm{tr}_{Bb}(|\mathrm{\Psi }_{Vab}\mathrm{\Psi }|)`$. Because of Eq. (11) we have

$`\xi _{Aa}`$ $`=`$ $`{\displaystyle \underset{l}{}}q_l|\mu _l_A\mu _l|\mathrm{tr}_b(|\nu _l_{ab}\nu _l|)`$ (20)
$``$ $`{\displaystyle \underset{l}{}}q_l|\mu _l_A\mu _l|\eta _a^l.`$
Define $`\gamma _A\mathrm{tr}_a(\xi _{Aa})=_lq_l|\mu _l_A\mu _l|`$. Then we have

$`E(|\mathrm{\Psi }_{Vab})=S(\xi _{Aa})`$ $``$ $`S(\gamma _A)+{\displaystyle \underset{l}{}}q_lS(\eta _a^l)`$ (21)
$``$ $`E_f(\varrho _V)+{\displaystyle \underset{l}{}}q_lE(|\nu _l_{ab})`$
$``$ $`E_f(\varrho _V)+E_f(\pi _{ab}).`$
The first inequality follows from the strong subadditivity of the entropy notestrong , as shown in Sh02 . In the second inequality we have used that $`S(\gamma _A=\mathrm{tr}_B(\varrho _V))E_f(\varrho _V)`$ notebound , also that $`S(\eta _a^l=\mathrm{tr}_b|\nu _l_{ab}\nu _l|)=E(|\nu _l_{ab})`$. Finally, the last inequality follows from the fact that $`\{q_l,|\nu _l_{ab}\}`$ is a (possibly non–optimal) realization of $`\pi _{ab}`$, $`\pi _{ab}=_lq_l|\nu _l_{ab}\nu _l|`$.

Thus, as illustrated in example 1, we can use this theorem to relate the asymptotic entanglement cost $`E_c`$ of the mixed states supported on some subspace $`VH_AH_B`$ to their entanglement of formation $`E_f`$. All that is needed is to identify subspaces $`VH_{AB}`$ that fulfill the above requirements.

Recall that, other than for two–qubit mixed states, the value of $`E_f`$ is only known in very few cases Wo98 . In this sense, another class of subspaces $`V^{}H_AH_B`$ of particular interest are those such that all their vectors are related by local unitary transformations. Since their reduced density matrices have the same spectrum, all these states are equally entangled. Let $`E(V^{})`$ denote their entropy of entanglement. Then, because a mixed state $`\rho _V^{}`$ supported on $`V^{}`$ is necessarily a mixture of vectors of $`V^{}`$, we conclude that $`E_f(\rho _V^{})=E(V^{})`$.

Example 2. Let us consider a two-qutrit system, $`H_A=H_B=𝒞^3`$, and the antisymmetric subspace $`V^{}H_AH_B`$, spanned by the vectors

$`|0_V^{}`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|1_A|2_B|2_A|1_B),`$
$`|1_V^{}`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|2_A|0_B|0_A|2_B),`$
$`|2_V^{}`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|0_A|1_B|1_A|0_B).`$ (22)
Notice that $`\mathrm{tr}_B(|\alpha _V^{}\beta |)=(\delta _{\alpha ,\beta }I_A|\beta _A\alpha |)/2`$. Therefore, for any vector $`|\varphi _V^{}=_\alpha c_\alpha |\alpha _V^{}`$, its reduced density matrix $`\rho _A^\varphi \mathrm{tr}_B|\varphi _V^{}\varphi |`$ can be expressed as a linear combination of the identity operator $`I_A`$ and the transposition of a projector onto the state $`|\varphi _A_\alpha c_\alpha |\alpha _A`$,

$$\rho _A^\varphi =\frac{1}{2}I_A\frac{1}{2}(|\varphi _A\varphi |)^T$$
(23)
(equivalently, any two antisymmetric states are related by local unitary transformations). It is then straightforward to obtain $`E_f(\rho _V^{})=1`$ ebit for any state $`\rho _V^{}(V^{})`$. Unfortunately, the subspace $`V^{}`$ does not meet the requirements of the theorem, and we can not compute $`E_c`$ for this family of mixed states.

Finally, in some cases one can combine the two results discussed in this paper to determine both $`E_f`$ and $`E_c`$, as the following examples illustrate.

Example 3. Let us consider a qubit-qutrit system, $`H_A=𝒞^2`$ and $`H_B=𝒞^3`$, and the subspace $`V^{\prime \prime }`$ spanned by

$`|0_{V^{\prime \prime }}{\displaystyle \frac{1}{\sqrt{3}}}(|0_A|2_B\sqrt{2}|1_A|0_B),`$
$`|1_{V^{\prime \prime }}{\displaystyle \frac{1}{\sqrt{3}}}(|1_A|2_B\sqrt{2}|0_A|1_B).`$ (24)
In this case $`\mathrm{tr}_B(|\alpha _{V^{\prime \prime }}\beta |)=(2\delta _{\alpha ,\beta }I_A|\alpha _A\beta |)/3`$, and therefore, for any vector $`|\varphi _{V^{\prime \prime }}=_\alpha c_\alpha |\alpha _{V^{\prime \prime }}`$, we find

$$\mathrm{tr}_B(|\varphi _{V^{\prime \prime }}\varphi |)=\frac{2}{3}I_A\frac{1}{3}|\varphi _A\varphi |,$$
(25)
where $`|\varphi _A_\alpha c_\alpha |\alpha _A`$. It follows that entanglement is constant in $`V^{\prime \prime }`$, $`E_f(\rho _{V^{\prime \prime }})=E(|\varphi _{V^{\prime \prime }})=H_2(1/3)`$. In addition, by noticing that the TPCP map $`^{\prime \prime }`$,

$$^{\prime \prime }(|\varphi _A\varphi |)\frac{2}{3}I_A\frac{1}{3}|\varphi _A\varphi |,$$
(26)
is entanglement–breaking, since it can be expanded as entbreak

$$^{\prime \prime }(X)=_{H_A}d|\varphi _A\mathrm{tr}(|\varphi _A\varphi |X)(I_A|\varphi _A\varphi |),$$
(27)
it follows, because of the theorem, that also $`E_c(\rho _{V^{\prime \prime }})=H_2(1/3)`$.

So far, we have calculated the entanglement cost for rank 2 density operators. Obviously, we can use our methods to determine this quantity for higher rank operators, as the following example shows.

Example 4. Let us consider $`H_A=𝒞^3`$, $`H_B=𝒞^6`$, and the subspace $`V^{\prime \prime \prime }`$ spanned by

$`|0_{V^{\prime \prime \prime }}`$ $``$ $`{\displaystyle \frac{1}{2}}(|1_A|2_B+|2_A|1_B+\sqrt{2}|0_A|3_B),`$
$`|1_{V^{\prime \prime \prime }}`$ $``$ $`{\displaystyle \frac{1}{2}}(|2_A|0_B+|0_A|2_B+\sqrt{2}|1_A|4_B),`$
$`|2_{V^{\prime \prime \prime }}`$ $``$ $`{\displaystyle \frac{1}{2}}(|0_A|1_B+|1_A|0_B+\sqrt{2}|0_A|5_B).`$ (28)
Since $`\mathrm{tr}_B(|\alpha _{V^{\prime \prime \prime }}\beta |)=(\delta _{\alpha ,\beta }I_A+|\beta _A\alpha |)/4`$, for any vector $`|\varphi _{V^{\prime \prime \prime }}=_\alpha c_\alpha |\alpha _{V^{\prime \prime \prime }}`$ we have

$$\mathrm{tr}_B(|\varphi _{V^{\prime \prime \prime }}\varphi |)=\frac{1}{4}I_A+\frac{1}{4}(|\varphi _A\varphi |)^T,$$
(29)
where $`|\varphi _A_\alpha c_\alpha |\alpha _A`$. It follows that the spectrum of the reduced density matrix for system $`A`$ is $`\{1/2,1/4,1/4\}`$, that is, the same for any pure state $`|\varphi _{V^{\prime \prime \prime }}`$. Therefore the entanglement is also constant in $`V^{\prime \prime \prime }`$, $`E_f(\rho _{V^{\prime \prime \prime }})=E(|\varphi _{V^{\prime \prime \prime }})=1.5`$ ebits. Finally, the TPCP map

$$^{\prime \prime \prime }(|\varphi _A\varphi |)\frac{1}{4}I_A+\frac{1}{4}(|\varphi _A\varphi |)^T,$$
(30)
is entanglement–breaking alternative , since when applied to the maximally entangled state $`(_{i=0}^2|i_A|i_C)/\sqrt{3}`$, where $`C`$ denotes an auxiliary system, the resulting state $`P_+/6`$, proportional to the projector $`P_+`$ onto the symmetric subspace of $`H_AH_C`$, is known to be separable Du00 . Consequently, the theorem implies that $`E_c(\rho _{V^{\prime \prime \prime }})=E_f(\rho _{V^{\prime \prime \prime }})=1.5`$ ebits.

Summarizing, we have shown that the entanglement of formation $`E_f`$ is additive for mixed states supported on a subspace such that tracing out one of the parties corresponds to an entanglement–breaking channel Sh02 . This has allowed us to evaluate the entanglement cost $`E_c`$ of several families of mixed states. A series of examples have been selected to illustrate these results.

Whether the entanglement of formation is additive for general mixed states remains an open question, which certainly deserves further investigation.

This work was supported by European Community, under project EQUIP (contract IST-1999-11053), grants HPMF-CT-1999-00200 (G.V.) and HPMF-CT-2001-01209 (W.D.) \[Marie Curie fellowships\] and through the ESF; by the Institute for Quantum Information GmbH, and by the United States of America through the NSF, under Grant. No. EIA-0086038.


format: text

text:
# Pair production of charged Higgs bosons at future linear 𝒆⁺𝒆⁻ colliders

## I Introduction

The exploration of the origin of mass and of the mechanism of electro-weak symmetry breaking represents one of the main issues on the physics agenda of future colliders. If the Higgs mechanism is indeed responsible for providing matter and force particles with their masses, the detailed exploration of the Higgs sector will most likely require the combination of data from the Lhc and lepton colliders. Today there are compelling indications that New Physics is required beyond the Standard Model (SM). In several SM extensions, a second complex scalar Higgs doublet is added yielding five physical Higgs bosons, two of which charged $`H^\pm `$, instead of just one. This is the case in the minimal supersymmetric SM extension (MSSM) and also in non-supersymmetric models (2HDM). The extra Higgs doublet introduces an additional fundamental parameter, describing the ratio of the vacuum-expectation values of the two doublets: $`\mathrm{tan}\beta =v_2/v_1`$.

The negative results of charged Higgs boson searches in the Lep-2 and Tevatron data and the indirect limits from the rate of the $`bs\gamma `$ process have set lower bounds on its mass, $`M_{H^\pm }`$, but the bulk of the parameter space remains still unexplored. At the Lhc, charged Higgs bosons can be observed but, according to present studies, only in the limited region of the parameters defined by $`\mathrm{tan}\beta >1030`$, depending on its mass and detectable decay modes.

In this paper, we discuss the potential of a high energy, high luminosity $`e^+e^{}`$ linear collider (LC) in the discovery of the charged Higgs boson and in the study of its properties. Since $`M_{H^\pm }`$ is not constrained in the models and can vary from the present 78.5 GeV Lep-2 limit lep  up to, and beyond, 1 TeV, different collider energies need to be considered. We concentrate on a TeV-class LC able to deliver $`e^+e^{}`$ collisions at $`\sqrt{s}`$ = 0.8 TeV, such as Tesla tesla  or an X-band LC nlc , and a multi-TeV LC operating at $`\sqrt{s}`$ = 3 TeV, such as Clic clic , and discuss the experimental issues specific to the different LC parameters.

Two mass values have been chosen for the charged Higgs boson: $`M_{H^\pm }`$ = 300 GeV with moderate $`\mathrm{tan}\beta `$, representative of SUSY scenarios with a relatively light Higgs sector, and $`M_{H^\pm }`$ = 880 GeV with $`\mathrm{tan}\beta `$ = 35, as point J from a set of SUSY benchmarks recently proposed msugra .

## II Production cross section for charged Higgs bosons in $`e^+e^{}`$ collisions

In $`e^+e^{}`$ collisions, charged Higgs bosons are pair produced via an intermediate photon or $`Z^0`$ boson. At tree level, the production cross section for $`e^+e^{}H^+H^{}`$ only depends on $`M_{H^\pm }`$ and on the centre-of-mass energy $`\sqrt{s}`$:

$$\sigma =\frac{e^4}{48\pi s}\left(1\frac{4m_H^2}{s}\right)^{3/2}\left(1+\frac{2c_V^{}c_V}{1m_Z^2/s}+\frac{c_{V}^{}{}_{}{}^{2}(c_V^2+c_A^2)}{(1m_Z^2/s)^2}\right)$$

However, in calculating the effective cross sections at high energy $`e^+e^{}`$ collisions, initial state radiation (ISR) and beam-beam effects must be taken into account. Because of the very small transverse dimensions of the beams at the interaction point, the electrons and positrons undergo significant radiation, in the field of the incoming beam, before collision (beamstrahlung). The average energy loss is of the order of 5% for $`\sqrt{s}`$=0.8 TeV at Tesla but reaches 31% for $`\sqrt{s}`$=3 TeV at Clic.

When these effects are taken into account, the tree level cross section is folded with the luminosity spectrum. For the low values of $`M_{H^\pm }`$ the effective cross section is enhanced compared to the tree level cross section, while close to the kinematical threshold only a smaller fraction of the energy spectrum is available and the effective cross section is reduced as exemplified in Figure 1. Corrections to the $`H^+H^{}`$ production cross section have been computed at one loop and found to become sizeable in some region of the parameter space hhcorr . A precise determination of the cross section can thus provide additional information on model parameters.

When considering only decays to SM particles, the dominant $`H^{}`$ decay modes are those involving the heaviest quark or lepton pairs accessible, i.e. $`\overline{t}b`$ and $`\tau ^{}\overline{\nu }_\tau `$111In the following the notation of a particle implies also its anti-particle, as shown in Figure 2 for two values of $`\mathrm{tan}\beta `$. The dependence of the dominant branching fractions on $`\mathrm{tan}\beta `$ can be exploited for determining this fundamental parameter of the theory. In this study we consider the fully hadronic $`e^+e^{}H^+H^{}(t\overline{b})(\overline{t}b)`$ channel, which leads to two $`W`$ bosons and four $`b`$ quarks in the final state and provides with a direct reconstruction of $`M_{H^\pm }`$.

## III Event Reconstruction

The Simdet parametrised detector simulation ref-simdet  has been used to account for the anticipated experimental resolutions of a LC detector. The $`b`$-tagging performance, which is of major importance in this analysis, has been included, assuming an efficiency $`ϵ_b`$ = 0.90. The characteristic signal final state, including four $`b`$ jets and the intermediate $`W`$ and $`t`$ mass constraints, guarantees an efficient rejection of most multi-fermion background processes. The cross section of the irreducible $`t\overline{b}\overline{t}b`$ background has been estimated using the Comphep program comphep  at 0.8 TeV and 3 TeV and found to be 5.5 fb and 1.5 fb respectively.

Hadronic jets have been reconstructed using the CAMJET clustering algorithm, for the 0.8 TeV analysis, and the Lund algorithm for 3 TeV, tuning the cut-off to maximise the number of eight jet events for the signal. In order to distinguish the production of charged Higgs bosons from the associated background processes and to accurately measure its mass, it is important to obtain a clean mass distribution of the multi-jet final states. Therefore only hadronic $`W`$ boson decays have been considered by removing events with significant transverse missing energy or tagged leptons.

The mass reconstruction has been performed as follows. In a first step, the presence of two $`W`$ bosons decaying hadronically has been tested. The assignment of the four non-$`b`$ tagged jets to their correct $`WW`$ pair has been done by choosing the combination which minimises the difference between the di-jet masses and the $`W`$ mass. In a second step, each $`t`$ quark has been reconstructed from a $`W`$ candidate paired with one of the four $`b`$ tagged jets, taking the combination giving the invariant mass closest to the top mass. Each of the charged Higgs bosons has then been reconstructed from a $`t`$ candidate paired with one of the two remaining $`b`$ tagged jets and selecting the configuration minimising the mass difference of the two $`tb`$ pairs. Events with no jet combination compatible with the $`W`$ or $`t`$ masses within the observed resolution have been discarded. An event kinematical fit has been finally applied to improve the mass resolution. The fit uses energy and momentum conservation, the $`W`$ and $`t`$ mass constraints and imposes $`H`$ boson equal mass.

### III.1 $`M_H`$ = 300 GeV at a TeV-class LC

At $`\sqrt{s}`$ = 0.8 TeV, the reconstructed mass peak has a width of 6.2 GeV with a natural $`H^\pm `$ width of 4.3 GeV, for $`M_{H^\pm }`$ = 300 GeV. This corresponds to a mass resolution of 1.5%. The signal reconstruction efficiency has been found to be 0.022, corresponding to 149 signal with 24 background events, in a signal mass window of $`\pm 2.5\sigma `$ around the mass peak. The good mass resolution achieved with the kinematical fit gives a signal significance $`S/\sqrt{B}`$=30 (see Figure 3). The product $`\sigma (e^+e^{}H^+H^{})\times \mathrm{BR}(H^{}\overline{t}b)`$ can be measured to an accuracy of 8.8%.

The maximum reach in mass for a $`5\sigma `$ discovery has been studied by analysing signal samples generated at different Higgs mass. The largest observable Higgs mass corresponds to about 350 GeV, corresponding to the detection of 47 signal events, with the same number from background for 1 ab-1 of data at $`\sqrt{s}`$ = 0.8 TeV.

At moderate $`\mathrm{tan}\beta `$ values, as well as in scenarios where the couplings to fermions are suppressed, the $`H^{}W^{}h^0`$ decay, leading to the same $`W^{}b\overline{b}`$ decay final state, may become sizeable. In order to study this process, a similar mass reconstruction procedure can be applied except for the top quark constraint being replaced by the $`h^0`$ one. By comparing the $`\chi ^2`$ of the fit for the two hypotheses, the $`H^{}\overline{t}b`$ and $`H^{}W^{}h^0`$ decays can be efficiently distinguished.

### III.2 $`M_H`$ = 880 GeV at a multi-TeV LC

At a multi-TeV collider the decay reconstruction is complicated by the larger beamstrahlung, as mentioned above, and by the $`\gamma \gamma \mathrm{hadrons}`$ background bkg  overlayed to each $`e^+e^{}`$ collision.

The reconstruction procedure for the charged Higgs boson follows closely that described above, but special care has been taken to make the result robust in presence of the $`\gamma \gamma `$ background, while the broader luminosity spectrum affects the kinematic fit, as discussed below. The additional hadrons generated by these $`\gamma \gamma `$ collisions, mostly affecting the forward regions, may either be merged into the jets coming from the $`H^\pm `$ decay or result in extra jets being reconstructed. It is therefore important to minimise the impact from this background source on the event reconstruction. In this analysis only the four leading non-$`b`$ jets have been considered, together with the $`b`$-tagged jets. In addition the resulting invariant masses of the intermediate states are biased towards higher values due to the additional hadron contribution. This needs to be taken into account by calibrating the mass response with simulation including the $`\gamma \gamma `$ background.

Due to the significant loss of energy of the colliding $`e^+`$ and $`e^{}`$, energy and momentum conservation constraints cannot be applied on the reconstructed system for the nominal $`\sqrt{s}`$. Instead the kinematical fit allows for an extra particle to be radiated but imposes its transverse momentum to be zero. The kinematic fit allows to improve, by a factor of two, the resolution on the mass reconstruction for $`H`$ candidates. Furthermore, the use of the kinematical fit, allows to correct for the effect of the $`\gamma \gamma `$ background resulting in a mass peak of the charged Higgs boson at the position of the generated mass (see Figure 4).

The signal event rate has been estimated for $`M_H=880\text{GeV}`$, first disentangling the $`\gamma \gamma `$ background and then assuming that the detector integrates the background from 15 bunch crossings (BX). In the first case with an integrated luminosity is 3 ab-1 and assuming BR($`Htb`$) = 1.0, the r.m.s. of the mass peak, after the kinematical fit, is 33 GeV with a $`H^\pm `$ natural width of 21 GeV. 47 $`e^+e^{}H^+H^{}(t\overline{b})(\overline{t}b)`$ signal events are expected to be reconstructed within the $`\pm 2.5\sigma `$ signal mass region, corresponding to a signal efficiency of $`0.02`$.

When the hadronic background, integrated over 15 BX is included, the peaks of the intermediate $`W`$ and $`t`$ states are shifted towards higher masses and get significantly wider. To account for this effect, the mass windows for the reconstructed $`W`$ boson and $`t`$ quark need to be shifted and made 50% wider. The fitted signal width becomes 42 GeV and there are 51 $`e^+e^{}H^+H^{}(t\overline{b})(\overline{t}b)`$ signal events with a reconstructed mass in the signal window, corresponding to an efficiency comparable to the previous case. The $`tbtb`$ background processes contribute 24 events for 3 ab-1. By scaling these results, the estimated maximum mass reach for a $`5\sigma `$ discovery with 3 ab-1 at 3 TeV is $``$1.0 TeV.

Beyond discovery the properties of the charged Higgs boson needs to be studied. These include the mass, the production rate and the decay branching fractions. The mass can be measured to an accuracy of 1% for 15 integrated BX and the product $`\sigma (e^+e^{}H^+H^{})\times \mathrm{BR}(H^{}\overline{t}b)`$ can be measured to an accuracy of $``$ 20%.

## IV Conclusions

The full reconstruction of the charged Higgs boson $`H^+H^{}t\overline{b}\overline{t}b`$ multi-jet final states has been studied for two scenarios: $`M_{H^\pm }`$ = 300 GeV at $`\sqrt{s}`$ = 800 GeV and $`M_{H^\pm }`$ = 880 GeV with $`\sqrt{s}`$ = 3 TeV. The analyses are based on the reconstruction of hadronic events with eight jets, four of which $`b`$-tagged. A kinematical fit taking into account the $`W`$ and $`t`$ intermediate mass constraints and energy-momentum conservation has been applied to reject the multi-fermion background and optimise the mass resolution for signal events. The Higgs mass can be reconstructed with an accuracy of $`1\%`$ and the product $`\sigma (e^+e^{}H^+H^{})\times \mathrm{BR}(H^{}\overline{t}b)`$ can be measured to an accuracy of 8.8% for $`M_{H^\pm }`$ = 300 GeV at 0.8 TeV and 15% for $`M_{H^\pm }`$ = 880 GeV at 3 TeV.

For further studies, a more efficient event selection needs to be envisaged. An interesting possibility is to tag the $`e^+e^{}H^+H^{}`$ event by reconstructing only one $`Htb`$ decay. This allows to increase the event selection efficiency and to study the $`H`$ decay in the hemisphere opposite to the tagged one in an unbiased way.

###### Acknowledgements.

We are grateful to E. Boos for his help with the Comphep package. The research activity of A.F. has been supported by a Marie Curie Fellowship of the European Community programme Improving Human Research Potential and the Socio-economic Knowledge Base under contract number HPMF-CT-2000-00865.


text:
# 1 (a) Snapshot configuration of the actual simulated time horizon for the one-dimensional one site per PE topology with 𝐿=10,000 PEs. The time horizon propagates “upwards” in the figure, and 𝑤 indicates the average width of the time horizon. (b) A short segment (50 sites) of the same configuration. The arrows on the left indicate the vertical scales in simulated time units (stu).

Going through Rough Times: from Non-Equilibrium Surface Growth to Algorithmic Scalability

G. Korniss,1 M.A. Novotny,2 P.A. Rikvold,3 H. Guclu,1 and Z. Toroczkai4
1Department of Physics, Applied Physics, and Astronomy, Rensselaer Polytechnic Institute,
110 8th Street, Troy, NY 12180-3590 USA
2Department of Physics and Astronomy and Engineering Research Center,
Mississippi State University, Mississippi State, MS 39762-5167 USA
3School of Computational Science and Information Technology, Department of Physics, and Center for Materials Research and Technology,
Florida State University, Tallahassee, FL 32306 USA
4Theoretical Division and Center for Nonlinear Studies, MS-B258
Los Alamos National Laboratory, Los Alamos, NM 87545 USA

ABSTRACT

Efficient and faithful parallel simulation of large asynchronous systems is a challenging computational problem. It requires using the concept of local simulated times and a synchronization scheme. We study the scalability of massively parallel algorithms for discrete-event simulations which employ conservative synchronization to enforce causality. We do this by looking at the simulated time horizon as a complex evolving system, and we identify its universal characteristics. We find that the time horizon for the conservative parallel discrete-event simulation scheme exhibits Kardar-Parisi-Zhang-like kinetic roughening. This implies that the algorithm is asymptotically scalable in the sense that the average progress rate of the simulation approaches a non-zero constant. It also implies, however, that there are diverging memory requirements associated with such schemes.

INTRODUCTION

Faithful and efficient simulation of complex systems with many interacting degrees of freedom is an important and challenging computational task. In a large class of systems the dynamic evolution is inherently stochastic and changes in the local configuration of the system occur randomly in space and time. The modeling of these systems can follow a “bottom-up” approach, starting with the definition of the “microscopic” dynamics. Then a master equation can be constructed with the corresponding transition probabilities. In most cases, for systems with many interacting degrees of freedom, even the numerical solution of the master equation (typically involving the numerical diagonalization of huge matrices) becomes insurmountable. This is when simulation becomes an invaluable tool for complex system modeling.

In the physics, chemistry, and biology communities these types of simulations are most commonly referred to as dynamic or kinetic Monte Carlo simulations. In computer science they are called discrete-event simulations. The updates (governed by the microscopic dynamics) in the (typically local) configuration of the system are considered discrete events. The basic notion of discrete-event simulation is that time is continuous and the discrete events occur instantaneously. Between events, the state (configuration) of the system remains unchanged. If the events occur at random instants of time, the dynamics can be referred to as asynchronous. Examples of such systems include the Ising model with the Glauber or Metropolis dynamics (the discrete events are the spin-flip attempts), cellular communication networks (call arrivals), spatial epidemic models (infections), financial markets (buy/sell orders), or internet traffic (packet transmission/reception). In this paper we focus on systems with general but short-ranged interactions on regular lattices and assume that the event dynamics can be described as a superposition of a large number of independent Poisson processes running in parallel (Poisson asynchrony).

When the size of the system becomes large, parallelization may be needed to obtain results within an acceptable time frame. Massively parallel simulation for complex systems with asynchronous dynamics, i.e., parallel discrete-event simulations (PDES), is a standard technique among computer scientists. It is somewhat surprising that despite PDES having a long history as far as applications and scalability are concerned , very few of the PDES techniques have filtered through, e.g., to the physics community. Even the simplest random-site update Monte Carlo schemes  where update attempts converge to Poisson arrivals in the large system-size limit, were long believed to be inherently serial. In this regard, Lubachevsky’s work  was rather illuminating, by illustrating how to apply the PDES scheme to the Ising model on a regular lattice with Glauber dynamics .

In a PDES scheme each processing element (PE) carries a subsystem of the full system via simple spatial decomposition. The difficulty of PDES is that the discrete events are not synchronized by a global clock since the dynamic is asynchronous. To put it simply, the paradoxical task is to (algorithmically) parallelize (physically) non-parallel dynamics. To achieve this, one must use the concept of local simulated times (or virtual times) and a synchronization scheme. The parallel algorithm must concurrently advance the local simulated time on each subsystem carried by a PE, without violating causality. In a “conservative” PDES scheme , only those PEs that are guaranteed not to violate causality are allowed to process their events and increment their local time. The rest of the PEs must “idle.” In an “optimistic” approach , the PEs do not have to idle, but since causality is not guaranteed at every update, the simulated history on certain PEs can become corrupted. This requires a complex “rollback” protocol to correct erroneous computations. Both simulation approaches lead to an evolving and fluctuating time horizon during algorithmic execution.

For systems which can be modeled on regular lattices with short-ranged interactions, the conservative scheme can be highly efficient. Recently, it was implemented for modeling magnetization switching  and the dynamic phase transition  in highly anisotropic thin-film ferromagnets. For example, the nearest-neighbor interaction implies that in order to ensure causality, PEs need to exchange their local simulated (or virtual) times only with “neighboring” PEs in the virtual PE topology. Communication times between PEs and possible idling due to the conservative synchronization protocol can be greatly suppressed by each PE carrying a large block of sites (spins) , yielding encouraging efficiencies and utilizations (fraction of non-idling PEs). Since there is a finite number of PEs in any architecture, one can typically be satisfied with these “experimental” observations. However, as the number of PEs available to simulating complex systems increases to tens of thousands, scalability questions become fundamental. Further, questions, such as how the utilization behaves in the asymptotic limit when the number of PEs goes to infinity, truly lie at the heart of any PDES scheme.

In this paper we address fundamental scalability questions for the general conservative parallel simulations for systems on regular lattices with short-range interactions. The way we tackle the problem, in some sense, goes opposite to the usual flow of a typical scientific modeling and simulation process. There, one develops advanced and sophisticated computational algorithms to study and understand systems in Nature. Here, by knowing how certain natural systems behave, we try to understand how advanced algorithms, PDES in particular, work: whether or not they are scalable, and how they can be optimized. Our approach is facilitated by a mapping  between non-equilibrium surface growth and the progress of the simulation (the evolution of the simulated time horizon). At the end, of course, based on the knowledge gained after answering fundamental scalability questions, one hopes to close this “loop” by devising and optimizing PDES schemes that can be used to investigate challenging problems in natural, artificial, or social systems. Analogies of a similar kind, e.g., that between phase transitions and computational complexity  also have turned out be highly fruitful to gain more insight in traditionally difficult problems in computer science. Also, exploiting analogies between the evolution of the time horizon and that of known physical systems appears to be rather helpful in understanding the performance for optimistic schemes as well. There is some evidence  that the time horizon in rollback-based schemes can exhibit self-organized criticality and power-law spatio-temporal correlations, which can be crucial to extract the scalability properties.

Based on our mapping between non-equilibrium surface growth and the progress of the simulation, it is clear that the tools and frameworks of modern statistical physics, in particular those of non-equilibrium interface/surface growth , can be extremely helpful in analyzing and understanding the asymptotic scalability properties of PDES schemes. To this end, one must look at the simulation scheme itself as an evolving interacting system of individual PEs where the synchronization rules among the PEs constitute the effective interaction. The evolution of this simulated time horizon, in particular its average progression rate and statistical spread, will determine the scalability properties of the corresponding PDES scheme.

In the next section we give an overview of the basic conservative PDES scheme . Then we characterize the morphological properties of the evolving random surface associated with the simulated time horizon. These findings yield direct implications for the scalability of the conservative algorithm for PDES.

THE BASIC CONSERVATIVE APPROACH

We consider a $`d`$-dimensional hypercubic regular lattice topology, where the underlying physical system has only nearest-neighbor interactions. However, our results hold for any short-range regular interaction pattern. In this paper we consider the case of simple Poisson asynchrony. Update attempts at each site are identical and independent Poisson processes (thus, the random simulated time increments between two successive update attempts are exponentially distributed) and are also independent of the state of the underlying physical system. The consequence of the former is that the evolution of the simulated time horizon completely decouples from the behavior and evolution of the underlying physical system. For simplicity, we discuss in detail the “worst-case” scenario, in which each PE carries one site (e.g., one spin). In this basic conservative scheme, each PE generates its own local simulated time for the next update attempt. (The actual update probabilities depend on the underlying systems, e.g., through energetics, but they do not affect the evolution of the time horizon.)

The set of local simulated times for $`N`$ PEs, $`\{\tau _i(t)\}_{i=1}^N`$, constitutes the simulated time horizon. Here $`t`$ is the discrete number of parallel steps simultaneously performed on each PE directly related to real/wall-clock time, or if the architecture operates in an asynchronous execution mode, $`t`$ is simply the continuous real time. On a regular $`d`$-dimensional hypercubic lattice $`N`$$`=`$$`L^d`$, where $`L`$ is the linear size of the lattice. In physics applications one typically specifies the initial configuration (i.e., at $`\tau `$=0) of the underlying physical system. This translates to $`\tau _i(0)`$$`=`$$`0`$ for every site for the initial condition of the parallel simulation. Then at each parallel update, only those PEs for which the local simulated time is not greater than the local simulated times of their nearest neighbors, can increment their local time by an exponentially distributed random amount, $`\eta _i(t)`$. The other PEs must idle. Without loss of generality we take independent, identically distributed (iid) exponential variables of mean one in simulated time units (stu), $`\eta _i(t)`$$`=`$$`1`$. Due to the continuous nature of the random simulated times, for $`t>0`$ the probability of equal-time updates for any two sites is of measure zero. The comparison with nearest neighbor simulated times and, if necessary, idling enforces causality. Also, at worst, the PE with the global minimum simulated time can make progress, so the algorithm is free from deadlock. For this basic conservative scheme, the theoretical efficiency or utilization (ignoring communication overheads) is simply the (average) fraction of non-idling PEs. This corresponds to the density of local minima of the simulated stochastic time horizon, which determines the average progress rate of the simulation. It can be written as

$$u(t)_N=\underset{i=1}{\overset{N}{}}\underset{jD_i^{\mathrm{nn}}}{}\mathrm{\Theta }\left(\tau _j(t)\tau _i(t)\right),$$
(1)
where $`D_i^{\mathrm{nn}}`$ is the set of nearest neighbors of site $`i`$, $`\mathrm{\Theta }()`$ is the Heaviside step function, and the $`\mathrm{}`$ denotes an ensemble average, i.e., an average over many independent simulations.

Another important aspect of the simulation is the width of the distribution of the local simulated times. This property can have serious effects on the “measurement part” of the algorithm, e.g. when one attempts to collect and compute simple statistics for the full underlying physical system “on the fly.” Therefore, one must determine the statistical spread (width) of the time horizon as was pointed out in Ref. . This quantity can be characterized by

$$w^2(t)_N=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left[\tau _i(t)\overline{\tau }(t)\right]^2,$$
(2)
where $`\overline{\tau }(t)`$$`=`$$`(1/N)_{i=1}^N\tau _i(t)`$ is the mean progress (“height”) of the time horizon. The behavior of the width, $`w^2(t)`$, alone typically captures and identifies the universality class of the non-equilibrium growth process .

EVOLUTION OF THE SIMULATED TIME HORIZON

The conservative synchronization protocol together with the virtual communication topology of the PEs (mimicking the interaction topology of the underlying physical system) fully specify the “microscopic” dynamics of the growth process associated with the evolution of the time horizon. The rules of the conservative synchronization can be compactly summarized as

$$\tau _i(t+1)=\tau _i(t)+\underset{jD_i^{\mathrm{nn}}}{}\mathrm{\Theta }\left(\tau _j(t)\tau _i(t)\right)\eta _i(t).$$
(3)
The above equation simply reflects that the local simulated time at site (PE) $`i`$ is incremented by an exponentially distributed random amount $`\eta _i(t)`$, provided that site $`i`$ is a local minimum of the time horizon. This stochastic time evolution equation is exact for the basic conservative scheme and can be easily simulated on a serial (!) workstation. Snapshots of the evolving and fluctuating time horizon obtained from direct simulations are shown in Fig. 1.

We discuss in detail the one-dimensional case ($`N`$$`=`$$`L`$) with periodic boundary conditions (i.e., ring topology). Replacing the Heaviside step function with a limiting smooth representation, one can perform standard coarse graining on Eq. (3) , yielding

$$\frac{\widehat{\tau }}{t}=\frac{^2\widehat{\tau }}{x^2}\lambda \left(\frac{\widehat{\tau }}{x}\right)^2+\widehat{\eta }(x,t).$$
(4)
Here $`\widehat{\tau }`$ is the coarse-grained “height” fluctuation, and the temporal and spatial derivatives are just the naive continuum interpretations of the differences, e.g., $`\widehat{\tau }/x`$$`=`$$`\tau _i\tau _{i1}`$. Similarly, $`\widehat{\eta }`$ is the coarse-grained noise, delta correlated in space and time. Higher order terms in Eq. (4) are dropped since they are irrelevant in a Renormalization Group (RG) sense. Equation (4) is the Kardar-Parisi-Zhang (KPZ) equation , which has turned out to be of central importance and have a wealth of applications for numerous artificial and natural growth processes over the past two decades, including molecular beam epitaxy, electrochemical deposition, fluid flow in porous media, and growth of bacterial colonies . Thus, we expect that the simulated time horizon exhibits kinetic roughening, the main feature of KPZ growth. Indeed, our direct simulation for the time horizon evolution, Eq. (3), confirms the coarse-graining approach \[Fig. 2(a) and (b)\]. There is a system-size dependent characteristic time scale, the crossover time, $`t_\times `$$``$$`L^z`$. For very early times, mostly microscopic details of the dynamics influence the width. For intermediate times, while $`t`$$``$$`t_\times `$ still, the width grows as a power law $`w^2(t)_L`$$``$$`t^{2\beta }`$, where $`\beta `$ is the growth exponent. For late times, $`t`$$``$$`t_\times `$, the width saturates for any finite system size. In this regime the surface reaches a steady-state evolution, and the fluctuations about the mean are stationary. The saturation or steady-state value of the width, however, scales as a power law with the system size, $`w^2(\mathrm{})_L`$$``$$`L^{2\alpha }`$, where $`\alpha `$ is the roughness exponent. The time horizon through its progress exhibits exactly the above scaling behavior with $`\beta `$$``$$`1/3`$ and $`\alpha `$$``$$`1/2`$, consistent with the exact one-dimensional KPZ exponents  \[Fig. 2(a)\]. This type of temporal and system-size scaling is consistent with the dynamic scaling hypothesis  and can be expressed through the Family-Vicsek scaling relation

$$w^2(t)_L=L^{2\alpha }f(t/L^z)$$
(5)
together with the important scaling law, $`\alpha `$$`=`$$`\beta z`$. Note that the scaling function $`f(x)`$ depends on $`t`$ and the linear system size $`L`$ only through the specific combination $`t/L^z`$, reflecting the importance of the crossover time $`t_\times `$. For small values of its argument $`f(x)`$ behaves as a power law, while for large arguments it approaches a constant

$$f(x)\{\begin{array}{cc}x^{2\beta }\hfill & \text{if }x1\hfill \\ \text{const.}\hfill & \text{if }x1\hfill \end{array}.$$
(6)
The existence of the above scaling function implies that if one plots the rescaled variables $`w^2(t)_L/L^{2\alpha }`$ vs $`t/L^z`$, then curves for different system sizes collapse for intermediate and late times. We have confirmed this data collapse for the simulated time horizon \[Fig. 2(b)\].

Direct simulation results for the average rate of progress of the time horizon (equivalent to the utilization in the PDES algorithm) for various system sizes are shown in Fig. 2(c). The utilization $`u(t)_L`$ decreases monotonically with time towards a long-time asymptotic limit well separated from zero, $`u(\mathrm{})_{\mathrm{}}`$$``$$`0.25`$. The fact that it cannot vanish in the infinite system-size limit can be argued based on an important universal feature of KPZ-like surfaces: The steady-state KPZ surface in one dimension is governed by the Edwards-Wilkinson Hamiltonian , i.e., it is essentially a random-walk profile . At coarse-grained length scales the local slopes become independent, yielding a non-zero average density of local minima, i.e., a non-zero average rate of progress of the simulation in the $`L`$$``$$`\mathrm{}`$ limit in the steady-state.

In higher dimensions we observe the same qualitative behavior as for $`d`$$`=`$$`1`$ . The surface roughens and saturates for any finite system, as seen in Fig. 3(a) and (c). Simultaneously, the density of local minima decreases monotonically towards its asymptotic ($`t`$$``$$`\mathrm{}`$) finite-size value \[Fig. 2(b) and (d)\]. Again, the steady-state density of local minima appears to be well separated from zero. For $`d`$$`=`$$`2`$ $`u(\mathrm{})_{\mathrm{}}`$$``$$`0.12`$, and for $`d`$$`=`$$`3`$ $`u(\mathrm{})_{\mathrm{}}`$$``$$`0.075`$. The $`u(\mathrm{})_{\mathrm{}}`$$``$$`𝒪(1/K)`$ behavior appears to be rather general , where $`K`$$`=`$$`2d`$ is the number of nearest neighbors on a regular lattice. Similar to the $`d`$$`=`$$`1`$ case, corrections to scaling are very strong, both for the surface width and the density of local minima. While for $`d`$$`=`$$`1`$ we were able to simulate large systems ($`L`$$``$$`10^3`$) to obtain the KPZ scaling exponents and the steady-state finite-size behavior of $`u(\mathrm{})_L`$, in higher dimensions the relatively small system sizes prevented us from extracting the scaling behavior of the width and the finite-size effects of the density of local minima. We conjecture that the simulated time horizon exhibits KPZ-like evolution in higher dimensions as well.

Next, we investigate in detail the steady-state scaling properties of the simulated time horizon for $`d`$$`=`$$`1`$, which directly translates to the asymptotic scalability properties of the corresponding PDES.

SCALING AND SCALABILITY

First, we provide further numerical evidence that the time horizon belongs to the KPZ universality class, in particular, in one dimension in the steady state, it is governed by the Edwards-Wilkinson (EW) Hamiltonian ,

$$_{\mathrm{EW}}\frac{1}{2}𝑑x(\widehat{\tau }/x)^2.$$
(7)
Then for small wave-vectors $`k`$, the average steady-state structure factor of the surface should behave as

$$S(k)=\tau _k\tau _k/L1/k^2,$$
(8)
where $`\tau _k`$$`=`$$`_{j=1}^Le^{ikj}\tau _j`$ is the spatial Fourier transform of the time horizon. This expectation is confirmed by direct simulations as shown in Fig. 4(a).

To further probe the universal properties of the surface in the steady state, we construct the width distribution $`P(w^2)`$, i.e., not just the average but the full probability density of the quantity in the brackets in Eq. (2) . This is a very powerful universal signature of rough surfaces, and in one dimension it can be tested against analytic results. The EW class is characterized by a universal scaling function, $`\mathrm{\Phi }(x)`$, such that $`P(w^2)=w^2^1\mathrm{\Phi }(w^2/w^2)`$ , where

$$\mathrm{\Phi }(x)=\frac{\pi ^2}{3}\underset{n=1}{\overset{\mathrm{}}{}}(1)^{n1}n^2e^{\frac{\pi ^2}{6}n^2x}.$$
(9)
Systems with $`L`$$``$$`10^3`$ show convincing data collapse \[Fig. 4(b)\] onto this exact scaling function.

Based on our results that the simulated time horizon exhibits KPZ-like kinetic roughening, we now address the implications for the scalability of conservative PDES schemes. We already argued that KPZ universal surfaces evolve toward a steady state where the coarse-grained local slopes, $`\widehat{\tau }/x`$, become independent \[see Eq. (7) in particular\]. This ensures that there is a finite density of local minima, so that the PDES algorithm progresses at a non-vanishing rate in the limit of an infinite number of PEs. Just as important for the utilization are the finite-size effects. Using exact calculations based on the Edwards-Wilkinson Hamiltonian in one dimension  and scaling arguments in higher dimensions , one can obtain the universal finite-size effects for the growth rate of generic KPZ-like processes

$$u(\mathrm{})_Lu(\mathrm{})_{\mathrm{}}+\frac{\text{const.}}{L^{2(1\alpha )}}.$$
(10)
The above equation can be used to estimate the utilization (average rate of progress of the simulation) and to extrapolate to the value for the infinite number of PEs, $`u(\mathrm{})_{\mathrm{}}`$. Equation (10) for the finite-size effects is in full agreement with the simulations \[Fig. 5(a)\] and yields $`u(\mathrm{})_{\mathrm{}}`$$`=`$$`0.246461(7)`$ for the one dimensional case ($`\alpha `$$`=`$$`1/2`$). While the actual asymptotic value of the density of local minima depends on “microscopic” measures, whether its asymptotic value vanishes or not, is fully governed by macroscopic characteristics and the corresponding universality class . We again emphasize that for the KPZ class this asymptotic value is non-zero.

The above findings for the density of local minima, which determines the average rate of progress of the simulation, imply that the “simulation part” of the conservative scheme is scalable. That is, if we run the simulation for long times, the average progress rate approaches a constant. However, the kinetic roughening exhibited by the time horizon has a disturbing implication: in the steady state the width (spread) of the simulated time horizon diverges with the number of PEs as

$$w^2(\mathrm{})_LL^{2\alpha }.$$
(11)
This scaling behavior for large $`L`$ is also confirmed by simulations \[Fig. 5(b)\], and it is contrary to the conclusions of Ref. . This property creates an additional difficulty for collecting statistics (e.g., to perform simple averages) “on the fly” during the course of the simulation. The diverging width means that the memory requirement per PE, for temporarily storing (buffering) data, diverges as we increase the number of PEs. In this sense we may call the “measurement part” of the bare conservative scheme asymptotically nonscalable. Thus, in an actual application, the programmer must implement some global synchronization or a moving “window” with respect to the global minimum of the time horizon . However such a “window” can have negative effects on the large $`L`$ utilization.

Along these lines of questioning, we are also interested in the extremal fluctuations of the time horizon. Namely, what is the typical size of the largest fluctuations above and below the mean, $`\mathrm{\Delta }_{\mathrm{max}}(t)`$$``$$`(\tau _{\mathrm{max}}(t)\overline{\tau }(t))`$ and $`\mathrm{\Delta }_{\mathrm{min}}(t)`$$``$$`(\overline{\tau }(t)\tau _{\mathrm{min}}(t))`$, where $`\tau _{\mathrm{max}}(t)`$ and $`\tau _{\mathrm{min}}(t)`$ are the global maximum and minimum simulated times among $`L`$ PEs, respectively. We found that in the steady state

$$\mathrm{\Delta }_{\mathrm{max}}^2(\mathrm{})_L\mathrm{\Delta }_{\mathrm{min}}^2(\mathrm{})_LL^{2\alpha },$$
(12)
i.e., they scale the same way as the width $`w^2(\mathrm{})_L`$ \[Fig. 5(b)\]. This should not come as a surprise, since this surface is highly correlated, dominated by long-wavelength fluctuations spreading over macroscopic length scales. This finding, again, is consistent with the extremal fluctuations found for general KPZ surfaces .

The above universal characteristics hold for the sensible and more efficient many sites per PE case and/or when the interaction and the corresponding communication pattern extends beyond nearest-neighbor PEs (but is still short ranged with a finite cutoff). For the case of many site per PE, however, the saturation occurs at a later time and the time horizon exhibits a substantially larger width .

SUMMARY AND OUTLOOK

We have studied the statistical properties of the basic conservative parallel scheme for regular lattice topologies. We found that the evolution of the simulated time horizon belongs to the well-known KPZ dynamic universality class of non-equilibrium surfaces. This type of growth is characterized by a non-zero density of local minima, implying a non-zero rate of propagation in the limit of an infinite number of PEs. We also determined the asymptotic finite-size corrections to this constant when the number of PEs is large but finite. Thus, the “simulation” part of the algorithm is scalable. Further, we showed that the spread (width) of the time horizon approaches a finite constant for a finite number of PEs, but this constant diverges as a power law in the limit of an infinite number of PEs. The same holds for the extremal fluctuations above and below the mean. This “macroscopic” roughness of the simulated time horizon means that the “measurement part” of the bare conservative scheme is not scalable. That is, there is an extra difficulty associated with taking statistical measurements “on the fly”. Intermittent data on each PE haved to be stored until all PEs reach the simulated time instant at which some statistics collection (e.g., simple averaging over the full physical application) is to be performed. The diverging spread of the time horizon implies a diverging storage need for this purpose on every PE. Thus, the programmer has to implement some global synchronization or windowing technique to limit the spread of the simulated time horizon in order not to exceed the memory constraint. By knowing exactly the finite-size dependence of the spread, for fixed $`L`$ one can determine the optimal time between global synchronizations or the optimal window size.

Our findings are universal in the sense that they hold for any short-range “interaction” topology for PEs on regular lattices. They are also valid in the case when each PE carries a block of sites. The asymptotic scaling behavior is again governed by the KPZ exponents, in such a way that for larger and larger blocks, there is a crossover from the almost “random deposition”  to KPZ-like growth at a later and later time.

We must mention that there have been earlier attempts to theoretically describe the scalability of the basic conservative PDES scheme. To obtain an analytically tractable scalability model, Greenberg et.al  introduced the $`K`$-random model. In this model at each update attempt, PEs compare their local simulated times to the local simulated times of $`K`$ randomly chosen PEs (rechosen at every update attempt). They showed that in the $`t`$$``$$`\mathrm{}`$, $`N`$$``$$`\mathrm{}`$ limit the average rate of progress of the simulation converges to a non-zero constant, $`1/(K+1)`$. They also showed that the evolution of the time horizon converges to a traveling-wave solution described by a finite width of the distribution of the local times. Finally, they suggested that convergence to a traveling-wave solution in the $`t`$$``$$`\mathrm{}`$, $`N`$$``$$`\mathrm{}`$ limit is universal and applicable for regular lattices as well. In obtaining this conclusion they made the assumption that replacing the “interaction” between nearest-neighbor PEs on a regular grid with the same interaction between randomly chosen PEs does not change the universality class of the time horizon. It does. Comparing the local time for each PE to $`K`$ randomly chosen others essentially turns the model in to a mean-field-like one where the time surface is short-range correlated and has a finite width in the limit of an infinite number of PEs. As we have shown in this paper, the time horizon of the conservative PDES for regular lattices and short-ranged interactions with finite cutoffs exhibits KPZ-like kinetic roughening.. This shows that the underlying communication topology of the PEs has crucial effects on the “universal” characteristics of the simulated time horizon.

Realizing the importance of the communication topology of the PEs, we are currently investigating how to turn the original conservative scheme on regular lattices into a fully scalable one, where both the “simulation” and the “measurement” parts are scalable, without the need for any global synchronization or windowing technique. Results of these studies will be published elsewhere.

ACKNOWLEDGMENTS

We thank B.D. Lubachevsky and Z. Rácz for useful discussions. We acknowledge the support of NSF through Grant Nos. DMR-9871455, DMR-0113049, and the support of the Research Corporation, Grant No. RI0761. Z.T. was supported by the U.S. Department of Energy under Contract No. W-7405-ENG-36.


text:
# Solitons in polarized double layer quantum Hall systems

Also at \]Department of Physics, Tbilisi State University, Tbilisi, Georgia

Also at \]Center for Semiconductor Physics in Nanostructures

Also at \]Center for Semiconductor Physics in Nanostructures

## Abstract

A new manifestation of interlayer coherence in strongly polarized double layer quantum Hall systems with total filling factor $`\nu =1`$ in the presence of a small or zero tunneling is theoretically predicted. It is shown that moving (for small tunneling) and spatially localized (for zero tunneling) stable pseudospin solitons develop which could be interpreted as mobile or static charge-density excitations. The possibility of their experimental observation is also discussed.

Much study has been devoted to the dynamical properties of charge-density excitations in double layer quantum Hall systems with total Landau filling factor $`\nu =1`$. An effective field-theoretical pseudospin model was developed in order to explain the experimentally observed spontaneous phase coherence even in the absence of tunneling moon ; yang ; brey  as well as for nonzero tunneling and in the presence of an applied in-plane magnetic field sheena . The main prediction stern ; balents ; fogler  of this model, the existence of the broken symmetry Goldstone mode fertig , has been confirmed by the observation spielman  of split off peaks in the tunneling conductance in the presence of in-plane magnetic field.

However, the study of double layer (pseudo) ferromagnets is usually restricted to the charge balanced state (see e.g. hanna  and references therein), although very recently the existence of a stable unbalanced state (the so called canted phase leo ; ramin1 ) in the presence of in-plane magnetic field was predicted. To date, static structures have dominated theoretical and experimental investigations. Only recently was the creation of a slowly moving pseudospin soliton investigated,kyr  but again this was done near the balanced state.

In this communication we study the dynamical problem of a strongly unbalanced double layer system in the presence of small or zero tunneling. Under such conditions the system is in highly nonequilibrium state but it is hard to reach the ground state if the interlayer tunneling amplitude is small, or even impossible if the amplitude is zero. On the other hand the nonlinear coupling induced by interlayer coherence is still present and thus the appearance of different exotic nonlinear creations might be expected. Indeed, as we will show below, stable, mobile or spatially localized pseudospin textures arise which could be described in terms of ac driven nonlinear Schrödinger (NLS) equation widely studied in literature (see e.g. Refs. baras ; baras2 ). This circumstance relates quantum Hall (pseudo)ferromagnets with completely different physical systems, such as yttrium iron garnet thin magnetic films kalinikos , long Josephson junctions lomdahl , relativistic electron-positron plasma berezhiani , etc., which are governed by the same NLS equation and where similar localized excitations are well studied both experimentally and theoretically.

The phenomelogical model for double layer quantum Hall (pseudo) ferromagnets, which is based upon microscopic considerations moon ; yang , is effectively described by the following Hamiltonian moon ; yang ; brey ; hanna :

$``$ $`=`$ $`{\displaystyle \frac{eV}{2}}n_z+{\displaystyle \frac{\rho _E}{2}}\left\{\left({\displaystyle \frac{n_x}{x}}\right)^2+\left({\displaystyle \frac{n_y}{x}}\right)^2\right\}+\beta n_z^2`$ (1)
$`\mathrm{\Delta }_{SAS}\left\{n_x\mathrm{cos}(Qx)+n_y\mathrm{sin}(Qx)\right\},`$
where $`\rho _E`$ is in-plane spin stiffness; $`\beta `$ is the a hard axis anisotropy; $`\mathrm{\Delta }_{SAS}`$ is a tunneling amplitude, $`QB_{}`$ is the “wave vector” associated with the tumbling of the electronic phase due to the in-plane magnetic field; $`V`$ describes the gate voltage; $`e`$ is an electron charge; $`\mathrm{}=1`$ and $`\widehat{n}(x,t)`$ is an order parameter unit vector. The component $`n_z(x,t)`$ is the local charge imbalance between the layers and can be expressed in terms of local filling factors of top and bottom layers, $`n_z=\nu _1\nu _2`$; the local filling factors $`\nu _1`$ and $`\nu _2`$ are proportional to the local electric densities $`N_1`$ and $`N_2`$ in the corresponding layers. For example, in the fully balanced case, $`\nu _1=\nu _2=1/2`$ (i.e. $`n_z=0`$) and $`N_1=N_2N`$. In typical experiments spielman  on double layer systems $`N=3.0\times 10^{10}cm^2`$. In the strongly unbalanced situation which is a subject of study in the present paper, one has $`n_z1`$ and thus the densities of top and bottom layers are close to $`2N`$ and $`0`$, respectively.

The time-space behavior of ordering vector could be described in terms of Landau-Lifshitz equationslichter :

$$\frac{\widehat{n}}{t}=\left(\widehat{n}\times \stackrel{}{H}_{eff}\right),$$
(2)
$$\stackrel{}{H}_{eff}=2\left\{\frac{}{\widehat{n}}\frac{}{x}\left[\frac{}{\frac{\widehat{n}}{x}}\right]\right\},$$

where $`\stackrel{}{H}_{eff}`$ is the effective (pseudo)magnetic field. We introduce the variables

$$n^\pm n_x\pm in_y$$
(3)
so that eq. (2) has the following form:

$`{\displaystyle \frac{n^+}{t}}`$ $`=`$ $`ieVn^++4i\beta n^zn^++2i\rho _En^z{\displaystyle \frac{^2n^+}{x^2}}`$ (4)
$`+2i\mathrm{\Delta }_{SAS}n_ze^{iQx}.`$
We will seek for the solution of this equation around $`n_z=1`$. This state corresponds to the highest energy value static solution. Any dynamical fluctuation around this state is the precession around $`n_z`$ which is allowed by the equations of motion. In the limit of small deviations from this state and small tunneling it is possible to apply the multiple scale approach taniuti ; ramaz :

$$n^+=\epsilon m^+(\xi ,\tau )e^{i(Qx\omega t)},$$

$$n_z=\sqrt{1|n^+|^2}=1\epsilon ^2\frac{|m^+|^2}{2},$$
(5)
$$\mathrm{\Delta }_{SAS}\epsilon ^3\stackrel{~}{\mathrm{\Delta }}_{SAS}$$

where $`m^\pm (\xi ,\tau )`$ is a slowly varying functions of the variables

$$\xi =\epsilon (xvt);\tau =\epsilon ^2t$$
(6)
and $`\epsilon `$ is a formal small parameter expressing the smallness or “slowness” of the object before which it appears, allowing for a multiple scale analysis of the problem. The range of acceptable $`\epsilon `$ will be determined in the course of the calculation. We are working in the regime:

$$|n^+|=|m^+|1.$$
(7)
To study the perturbative solution we substitute eq.(5) into the eq.(4), and collect terms of the same order of $`\epsilon `$. In the first order of $`\epsilon `$ the following equality is obtained:

$$\omega =eV+2(\rho _EQ^22\beta ),$$
(8)
while to second order we get the expression for the group velocity $`v`$:

$$v=4\rho _EQ.$$
(9)
Finally, to third order in $`\epsilon `$ we come to the following nonlinear equation:

$`i{\displaystyle \frac{m^+}{\tau }}+2\rho _E{\displaystyle \frac{^2m^+}{\xi ^2}}+(\rho _EQ^22\beta )|m^+|^2m^+=`$
$`2\stackrel{~}{\mathrm{\Delta }}_{SAS}e^{i\omega t},`$ (10)
which in the case of

$$\rho _EQ^2>2\beta ;\frac{eV+2\rho _EQ^24\beta }{\rho _EQ^22\beta }\epsilon ^2$$
(11)
reduces to the exactly solvable ac driven NLS equation considered in Refs. baras ; baras2 .

The physical meaning of the second expression from (11) together with definition (8) is that one has a near resonance driving of pseudospin mode with a wavenumber $`Q`$. The possibility of soliton, chaotic or other exotic solutions of the NLS equation (10) is studied in details in Refs. baras ; baras2  and we direct the reader there for more information. Here we only mention that the stable solutions of equation (10) could be interpreted as envelope solitons of order parameter $`\widehat{n}`$ moving with a group velocity proportional to applied in-plane magnetic field \[see expression (9)\] and characterized by carrier wave vector $`Q`$ which is also proportional to in-plane magnetic field. Note that according to the restriction of eq.(11), the in-plane magnetic field should be nonzero.

In Fig. 1 this excitation is presented as a moving bump of electric density through each layer. Note that envelope soliton nature only characterizes the transverse component of order parameter, while the physically measurable quantity - local charge imbalance - is just an ordinary soliton.

According to the general results of Refs. baras ; baras2  the stability of the mentioned soliton is restricted to the following condition:

$$h\frac{\mathrm{\Delta }_{SAS}\sqrt{\rho _EQ^22\beta }}{(eV+2\rho _EQ^24\beta )^{3/2}}<\sqrt{\frac{1}{27}}.$$
(12)
This quantity stands for the driving force in NLS equation baras  and therefore defines the soliton amplitude.

Here we also present the solutions for simple, but physically very important case of zero tunneling ($`\mathrm{\Delta }_{SAS}=0`$), which do not directly follow from the exact solutions for nonzero $`\mathrm{\Delta }_{SAS}`$ presented in Ref. baras . In the case of zero tunneling the in-plane magnetic field has no effect on the system. For simplicity we assume also $`V=0`$. Then the solution of the initial equation of motion (4) is sought for in the following simple form:

$$n^+=\epsilon m^+(x,t)e^{4i\beta t},n_z=1\epsilon ^2\frac{|m^+|^2}{2}.$$
(13)
Substituting this expression into Eq. (4) we come to the ordinary NLS equation:

$$i\frac{m^+}{t}+2\rho _E\frac{^2m^+}{x^2}2\beta |m^+|^2m^+=0,$$
(14)
which has a stable “dark soliton” solution abl . Since the physically measurable variable is $`n_z`$ (the local charge imbalance between the layers) we present here only an expression for its profile:

$$n_z=1\frac{D^2}{2}\left|\frac{\sqrt{1A^2}}{A}+i\mathrm{tanh}\left\{\frac{x}{\mathrm{\Lambda }}\right\}\right|^2,$$
(15)
where $`D`$ is a soliton amplitude; $`A`$ denotes the contrast of dark soliton (if $`A=1`$ one has a “black” dark soliton and “gray” dark for $`0<A<1`$) and soliton width $`\mathrm{\Lambda }`$ is defined as follows:

$$\mathrm{\Lambda }=\sqrt{\frac{2\rho _E}{\beta }}\frac{1}{D}.$$
(16)
The difference between gray dark and black dark solitons can be expressed the in terms of the charge densities in the top and bottom layers (see also Fig. 2): In a black dark soliton solution there exists some point along $`x`$ axis where the local electric density of top layer reaches the value $`2N`$, while the local density of bottom layer in the same point approaches zero. On the other hand, there always exists some deviation of the charge densities from $`2N`$ and $`0`$ for top and bottom layers, respectively, for the gray dark soliton solution.

Unlike the previous case of small tunneling where the soliton is characterized by nonzero carrier wave number $`Q`$, in the case of zero tunneling that is not the case and consequently transverse components of ordering vector have the similar form to the $`z`$ component and are not presented in Fig. 2.

We should emphasize again one more difference between two considered cases: In the case of nonzero (but small) tunneling the envelope soliton amplitude is linked with the driving coefficient $`\mathrm{\Delta }_{SAS}`$, while in the absence of tunneling the dark soliton amplitude is arbitrary within the restriction (7).

The above considerations are only valid for quasi-one dimensional samples, i.e. quantum Hall bars. Thus the main requirement for stability of soliton solutions is that the transverse dimension of the double layer system should be less than soliton width. In this case the large wavelength modulations do not grow (which destabilize the soliton) and solitons remain stable zakharov ; kivshar . Future work will investigate their stability.

In conclusion, we have found stable (mobile or static) charge-density excitations in strongly unbalanced double layer systems. For their observation the double layer quantum Hall system should be prepared so as $`n_z1`$. Then any perturbation will cause the appearence of either standing (zero tunneling) or moving (for small tunneling and applied in-plane magnetic field) charge density train of “bumps”, which could be easily detected experimentally. The perturbation could be induced by temporary application of local electric field perpedicular to the layers.

The work at the University of Oklahoma was supported by the NSF under grant No. EPS-9720651 and a grant from the Oklahoma State Regents for higher education. R. Kh. was supported by the NSF-NATO visiting scientist fellowship award No DGE-0075191. R.A. was supported by NSF MRSEC Grant DMR-0080054.


text:
# Mesoscopic Rectifiers Based on Ballistic Transport: Playing off Classical against Quantum Mechanics

## Abstract

Recent experiments on symmetry-broken mesoscopic semiconductor structures song1 ; song2  have exhibited an amazing rectifying effect in the transverse current-voltage characteristics with promising prospects for future applications. We present a simple microscopic model, which takes into account the energy dependence of current-carrying modes and explains the rectifying effect by an interplay of fully quantized and quasi-classical transport channels in the system. It also suggests the design of a ballistic rectifier with an optimized rectifying signal and predicts voltage oscillations which may provide an experimental test for the mechanism considered here.

The decreasing size of nanofabricated structures opens up new possibilities for mesoscopic semiconductor devices by exploiting ballistic transport in the 2-dimensional electron gas and quantization of confined electrons. Details of the geometry are crucial for the functioning of such devices. In the geometry of Fig. 1a, which was used in an experiment by Song et al. song1 , the symmetry was broken on purpose by introducing a triangle as shown. When current is injected at the source (S) and drawn out at the drain (D), a majority of charge carriers is deflected towards the lower voltage probe (L). Naively thinking and in the spirit of the Hall effect one might expect that a voltage difference builds up between the upper (U) and lower (L) voltage probe. Since the sample is symmetric with respect to the exchange of source and drain, the same voltage difference would arise on reversing the current and the sample would work as a rectifier.

On second thought, however, our physical understanding of mesoscopic systems leads us not to expect any voltage drop from top to bottom at all. As will be argued in detail below, symmetry considerations and the Landauer-Büttiker formalism LBF  in its most common linear form do not allow for such a voltage drop. It thus came as a surprise, when the experiment by Song et al. song1  revealed a rectifying effect. An interpretation of the experimental result in Refs. song1 ; song2  was based on the assumption of dissipation within the sample leading to self-consistent electric fields and a current dependence of the transmission probabilities through the sample. A phenomenological ansatz was made, as this current dependence would be exceedingly complicated to calculate in a microscopic model.

In the present letter we present a microscopic explanation of the rectifying effect which explicitely takes into account the energy dependence of the number of transverse modes for a system of two voltage probes. It does not require the existence of dissipation in the leads, instead the rectifying effect here originates in the interplay between purely quantized and quasi-classical transport in different channels of the system, which exhibit different energy dependences. A calculation of the tranverse current-voltage characteristics involving this mechanism shows good agreement with the experimental observations. This mechanism also leads to a prediction for the design of a ballistic rectifier with an optimized rectifying signal. For strong currents it predicts a reversal and even oscillations of the transverse voltage, which may provide a test for the explanation presented here.

Transport in mesoscopic systems like the one in Fig. 1a is typically described by the linear Landauer-Büttiker-formalism LBF

$$I_i=\frac{2e}{h}\left[(M_iR_i)\mu _i\underset{ji}{}\overline{T}_{ij}\mu _j\right].$$
(1)
Here $`I_i`$ is the net current in lead $`i`$ connecting the sample to a reservoir (contact) with chemical potential $`\mu _i`$. The leads are assumed to be ideal quantum leads with $`M_i`$ modes. $`R_i=\overline{T}_{ii}`$ is the reflection coefficient, which describes back-scattering from the sample into lead $`i`$ and $`\overline{T}_{ij}`$ are the transmission coefficients form lead $`j`$ into lead $`i`$. Transport across the sample is assumed to be purely elastic and dissipation and equilibration only take place in the reservoirs.

A prominent result of the Landauer-Büttiker formalism is the *reciprocity relation*

$$R_{ij,kl}(B)=R_{kl,ij}(B).$$
(2)
Here $`R_{ij,kl}=V_{kl}/I_{ij}`$ is the resistance obtained by dividing the voltage $`V_{kl}`$ measured between contact $`k`$ and $`l`$ by the current $`I_{ij}`$ flowing from contact $`j`$ to $`i`$. For the system of Fig. 1a at zero magnetic field we thus have

$$R_{UL,SD}=R_{SD,UL}.$$
(3)
Because of the symmetry of the system there can be no voltage build-up between $`S`$ and $`D`$ if the current is flowing from $`L`$ to $`U`$. Thus $`R_{UL,SD}=0`$ and hence by means of Eq. (3) we would expect $`V_{UL}`$ to vanish identically – in contrast to the experimental findings. (We would of course obtain the same results by solving the Eqs. (1) directly.) To overcome this apparent contradiction Song et al. song1 ; song2  suggested to include a phenomenological current dependence of the transmission coefficients due to dissipation inside the sample. In our treatment the necessary nonlinearity arises in the transport equations, if we allow for varying numbers of modes in the leads. We will show below that one needs not give up the conceptually attractive assumption of purely elastic transport inside the sample, that was so very successful in describing a wide range of experiments on transport in mesoscopic system (for reviews see e.g. mesotrans1 ; mesotrans2 ; mesotrans3 ).

To achieve this aim we use the Landauer-Büttiker formalism in a more general form datta . The current per unit energy (a quantity, which we will call *current density* in the following for simplicity) injected into the sample from reservoir $`i`$ through lead $`i`$ at energy $`E`$ is

$$i_i^+(E,\mu _i,T)=\frac{2e}{h}M_i(E)f_i(E).$$
(4)
Here $`f_i(E)=f(E,\mu _i,T)`$ is the Fermi-distribution in reservoir $`i`$ at temperature $`T`$. The outgoing current density in lead $`i`$ is

$$i_i^{}(E,\{\mu _l\},T)=\frac{2e}{h}\left[R_i(E)f_i(E)+\underset{ij}{}\overline{T}_{ij}(E)f_j(E)\right].$$
(5)
If we assume, that the transmission probabilities $`T_{ij}`$ are independent of energy and the mode number, we can write $`\overline{T}_{ij}(E)=T_{ij}M_j(E)`$ and Eq. (5) becomes

$$i_i^{}(E,\{\mu _l\},T)=\underset{j}{}T_{ij}i_j^+(E,\mu _j,T).$$
(6)
The incoming and outgoing currents are, respectively,

$$I_i^\pm =_{\mu _0}^{\mathrm{}}i_i^\pm (E)𝑑E,$$
(7)
where $`\mu _0`$ is an auxiliary quantity that is small enough so that $`f(\mu _0,\mu _i,T)1`$ holds for all $`i`$, but is otherwise arbitrary and will not show up in any measurable quantity. The balance equation for the current source with net current $`I`$ then is $`I_S^+I_S^{}=I`$ and for the drain $`I_D^+I_D^{}=I`$. A voltage probe is characterized by zero net current, i.e. $`I_i^+=I_i^{}`$.

The number of modes $`M_j(E)`$ can be written as

$$M_j(E)=\underset{n}{}\mathrm{\Theta }(E\epsilon _{j,n}),$$
(8)
where the $`\epsilon _{j,n}`$ are the energy eigenvalues of the transverse modes in lead $`j`$ and $`\mathrm{\Theta }(x)`$ denotes the Heaviside step function. For leads with a hard wall (i.e. box-like) cross section of width $`W_j`$ and for electrons of effective mass $`m^{}`$ we have

$$\epsilon _{j,n}=\frac{(\mathrm{}\pi n)^2}{2m^{}W_j^2}.$$
(9)
In this case the number of modes can also be expressed as

$$M(E)=\text{Int}\left[\frac{W_j}{\lambda (E)/2}\right]=\text{Int}\left[\frac{W_j\sqrt{2m^{}E}}{\mathrm{}\pi }\right],$$
(10)
where $`\text{Int}[]`$ denotes the integer part and $`\lambda (E)=h/\sqrt{2m^{}E}`$ is the de Broglie-wavelength of the electron at energy $`E`$. The incoming currents can be expressed as

$$I_j^+(\mu _j,T)=\frac{2e}{h}\underset{n}{}\left\{kT\mathrm{ln}\left[e^{(\epsilon _{j,n}\mu _j)/kT}+1\right]+\mu _j\epsilon _{j,n}\right\},$$
(11)
and the outgoing currents accordingly by integrating Eq. (6). These formulas were used in obtaining the numerical results presented below.

Let us now examine the very simple setup of Fig. 1b. A voltage probe reservoir (P) is connected via two identical ideal leads to source and drain. To simplify the calculations let us assume zero temperature, i.e. $`f(E,\mu _i,0)=\mathrm{\Theta }(\mu _iE)`$. Let us further assume $`\mu _S\mu _P\mu _D\mu _0`$, where $`\mu _S`$ and $`\mu _D`$ are given and $`\mu _P`$ is to be determined, and let us distinguish the cases of narrow and wide leads.

(a) If $`M(E)=M=`$constant in the range between $`\mu _0`$ and $`\mu _S`$ (this corresponds to narrow leads), the outgoing and incoming current densities from and to the probe reservoir are simply given by: $`i_P^+(E)=2(2e/h)M\mathrm{\Theta }(\mu _PE)`$ (the factor 2 arises because two leads are connected to the same reservoir) and $`i_P^{}(E)=(2e/h)\left[\mathrm{\Theta }(\mu _SE)+\mathrm{\Theta }(\mu _DE)\right]`$. These are trivially integrated and the current balance of the voltage probe $`I_P^+=I_P^{}`$ reads

$$\frac{4e}{h}M(\mu _P\mu _0)=\frac{2e}{h}M(\mu _S+\mu _D2\mu _0)$$
(12)
and thus the chemical potential $`\mu _P`$ of the voltage probe is independent of M and given by

$$\mu _P=(\mu _S+\mu _D)/2.$$
(13)
(b) If the leads are assumed to be wide compared to the Fermi wavelength $`\lambda _F`$, the number of modes will increase even under small changes in energy. In the case of a hard wall channel of width $`W>>\lambda _F`$ we can approximate Eq. (10) by a smooth function, i.e.

$$M(E)=Int\left[\frac{W}{\lambda (E)/2}\right]\frac{\sqrt{2m^{}}W}{\pi \mathrm{}}\sqrt{E}.$$
(14)
Since $`M(E)`$ grows as $`\sqrt{E}`$, the channels described in this approximation may be called *quasi-classical*, because the classical energy surface likewise increases as the square-root of the energy. The current densities are $`i_P^+(E)=2Q\sqrt{E}\mathrm{\Theta }(\mu _PE)`$ and $`i_P^{}=Q\sqrt{E}[\mathrm{\Theta }(\mu _SE)+\mathrm{\Theta }(\mu _DE)]`$ with $`Q=\sqrt{16m}eW/h^2`$. The current balance now reads

$$\frac{4Q}{3}\left(\mu _P^{3/2}\mu _0^{3/2}\right)=\frac{2Q}{3}\left(\mu _S^{3/2}+\mu _D^{3/2}2\mu _0^{3/2}\right),$$
(15)
and we find

$$\mu _P=\sqrt[3]{\left(\mu _S^{3/2}+\mu _D^{3/2}\right)/2}^2>(\mu _S+\mu _D)/2,$$
(16)
i.e. the potential $`\mu _P`$ deviates from the mean Eq. (13). As illustrated in Fig. 3 this is because the net current is transported from $`S`$ to $`P`$ by a larger number of modes in the energy interval $`[\mu _P,\mu _S]`$ than from $`P`$ to $`D`$ in the energy interval $`[\mu _D,\mu _P]`$ and thus $`\mu _P`$ raises above the mean to compensate for the additional current.

Equations (13) and (16) suggest how to create a rectifier based on ballistic transport: As sketched in Fig. 1c let us consider two separate pathways from source to drain, each via a voltage probe as in the above example (ignoring the third dotted pathway for the moment). If on both paths the number of modes is constant, the voltage probes will each be at the mean chemical potential between source and drain, i.e. there will be no voltage drop from top to bottom (this is what we showed based on the reciprocity relation in the beginning). The same holds true, if the current density grows identically with energy in both channels, e.g. for two quasi-classical channels. If on the other hand one path is narrow, i.e. has constant mode number $`M_1`$, whereas the other is wide with an increasing number of modes $`M_2(E)`$, we can observe a voltage drop between the probes, as indicated by the interval marked by the fat line in Fig. 3. Reverting the current yields the same voltage drop due to symmetry and thus we achieve the rectification of the signal.

In order to explain the experimental results of Ref. song1  we need to assume a third channel connecting the probes U and L in Fig. 1c (dotted lines) with $`M_3(E)`$ modes, since there is a direct connection between U and L in Fig. 1a. Now we consider finite temperatures again and choose $`M_1=1=`$ const. in the energy window, whereas for $`M_2(E)`$ and $`M_3(E)`$ we use the quasi-classical Eq. (14) with $`M_2(E_F)=20`$ and $`M_3(E_F)=15`$. These mode numbers were estimated by classical numerical calculations of the transmission coefficients $`\overline{T}_{US}`$ and $`\overline{T}_{LS}`$ in the geometry of Fig. 1a modeling the experiment. Curve A in Fig. 2 shows the resulting rectifying signal of this setup. If we give up the quasi-classical approximation for the wide channels and use the explicit sums instead, we obtain curve B, which nicely agrees with the experimental result (note that the experimental curve is not entirely symmetric due to an unintentional slight asymmetry of the sample about the vertical axis).

Based on the above mechanism we can now suggest the design of a ballistic rectifier with an optimized rectifying signal. If one manages to suppress the third channel between U and L, the voltage drop $`V_{UL}`$ can be enhanced. This is shown by curve C of Fig. 2 for the model of Fig. 1c without the dotted channel and using the quasi-classical approximation for the wide channel. Note that this curve can easily be calculated analytically for T=0. Using Eq. (13) and (16) one obtains $`\mu _U`$ and $`\mu _L`$, and hence the voltage difference, as a function of $`\mu _S`$, while leaving $`\mu _D`$ constant. The total current $`I`$ is given by the sum of the individual currents from the source S to U and to L, i.e. $`I=2Q/3(\mu _S^{3/2}\mu _L^{3/2})+2e/h(\mu _S\mu _U)`$, which is the current-voltage characteristics in analytic form. If again we give up the quasi-classical approximation for the wide leads, we find an even stronger rectifying signal as shown in curve D. As a realization of such a system which suppresses the third channel we suggest the geometry shown in Fig. 1d.

An interesting phenomenon arises in these structures, which combine narrow and wide leads, when the narrow channel opens up a new mode within the energy window: the voltage $`V_{UL}`$ undergoes a change in sign! This is demonstrated in Fig. 4, where the dotted line again shows curve D from Fig. 2. The solid line corresponds to slightly wider leads in the upper channel whereby $`\mu _S`$ becomes larger than $`\epsilon _2`$ of the narrow leads. This curve will eventually turn back to negative voltages as $`I`$ (following $`\mu _S`$) increases and shoot up again, when the third mode opens in the upper channel. This change in sign or even oscillations should be observable in an appropriate experiment, which would provide a test for the mechanism of ballistic rectifiers presented in this paper.

###### Acknowledgements.

The authors would like to thank Roland Ketzmerick for stimulating discussions. This research was supported in part by the the National Science Foundation under Grant No. PHY99-07949.