text:
# Dynamics of Competitive Evolution on a Smooth Landscape

## Abstract

We study competitive DNA sequence evolution directed by in vitro protein binding. The steady-state dynamics of this process is well described by a shape-preserving pulse which decelerates and eventually reaches equilibrium. We explain this dynamical behavior within a continuum mean-field framework. Analytical results obtained on the motion of the pulse agree with simulations. Furthermore, finite population correction to the mean-field results are found to be insignificant.

Competitive evolution such as breeding has been practiced for ages. With recent advances in molecular biology, this method is widely used to develop novel proteins and DNA sequences for a variety of applications REF:Arnold2001 . The basic idea of competitive molecular evolution is straightforward: in each generation, a number of molecules with certain desired characteristics are selected from the population; they are then diversified (via point mutation and/or recombination REF:Stemmer1994 ) and amplified back to the original population size. The “speed” of evolution as well as the final equilibrium distribution depend on a variety of factors such as the mutation rate, selection strength, molecule length, and population size. A systematic quantitative understanding of these dependencies is lacking thus far. Such understanding is not only of theoretical interest, but also helpful in improving the efficiency of the breeding processes. In this study, we develop a theoretical model for the simplest type of competitive evolution involving only point mutations on a smooth landscape. We achieve an understanding of this model with concepts and techniques developed in the study of front propagation REF:Collet1990 .

To make the discussion concrete, we focus on the in vitro evolution of DNA sequences due to competitive binding to proteins. An example of such a system is the recent experiment of Dubertret et. al. REF:Dubertret2001 , where DNA sequences are selected competitively according to their relative affinities for the lac-repressor protein. In this experiment, selection is accomplished by coating a beaker with lac-repressor molecules followed by subsequent washing, so that only the strongly-bound sequences remain. Mutation and amplification are then accomplished by multiple stages of polymerase chain reaction REF:Lodish2000 . While the experiment of Ref. REF:Dubertret2001  easily accomplished the goal of finding the best binding sequence starting from a pool of random sequences in a few generations, the shortness of the binding sequence \[$`20`$ base pairs (bp)\] makes it difficult to explore the interesting dynamics of the competitive evolution process. In our study we consider the evolution process of Ref. REF:Dubertret2001  applied to much longer sequences so that that the steady state dynamics can be examined. An example of such a system might be the histone-octamer, which is known to bind DNA sequences of $`147`$ bp REF:Widom1999 .

We consider a pool of $`N`$ DNA sequences of length $`L`$. Each sequence $`\stackrel{}{S}=(b_1,b_2,\mathrm{},b_L)`$ of nucleotides $`b_i`$ is subject to independent single-nucleotide mutations at a rate $`\nu _01`$ per nucleotide per generation. Selection is accomplished through protein-DNA binding. Let the binding energy of a sequence $`\stackrel{}{S}`$ to the protein be $`E_\stackrel{}{S}`$ and let the fraction of such sequences in the pool be $`n_\stackrel{}{S}`$. Assuming thermodynamic equilibrium for the binding process, the selection function is simply the binding probability, given by the Fermi function REF:VonHippel  $`P(E_\stackrel{}{S},\mu )=1/[1+\mathrm{exp}(E_\stackrel{}{S}\mu )]`$, where $`mu`$ is the chemical potential and all energies are expressed in units of $`k_\mathrm{B}T`$. Here $`\mu `$ serves as a (soft) selection threshold. and is determined by the fraction $`\varphi `$ of DNA sequences that remain bound to the proteins after selection, i.e., $`_\stackrel{}{S}P(E_\stackrel{}{S},\mu )n_\stackrel{}{S}=\varphi `$. It can be controlled by either the number of available proteins or, as in the experiment REF:Dubertret2001 , by the washing strength. The fraction $`\varphi `$ can be varied from $`\varphi 1`$ (weak competition) to $`\varphi 0`$ (strong competition). We define the evolution process iteratively whereby in each round, $`N`$ daughter sequences are chosen from the existing pool according to $`P(E_\stackrel{}{S},\mu )`$, and then point mutation are introduced with rate $`\nu _0`$ to generate the new sequence pool.

Finally we need to specify the binding energy $`E_\stackrel{}{S}`$. We assume that each nucleotide taking part in the binding contributes independently, and adopt a “two-state” model REF:VonHippel  which assigns an energy penalty $`ϵ`$ (of the order of a few $`k_\mathrm{B}T`$’s) for each nucleotide which does not match the one the protein prefers. This form of binding energy has been shown to work reasonably well for specific systems fields  and has been argued to hold for a wide class of regulatory proteins gmh , nonspecific . Given this energy model, a DNA sequence with $`r`$ mismatches has a binding probability $`P(r,r_0)=1/\left[1+e^{ϵ(rr_0)}\right]`$, where $`r_0\mu /ϵ`$ is the selection threshold in the “mismatch space” $`r`$REF:Gerland2001 .

The above evolution model is easily implemented on a computer. We fix three of the model parameters at $`N=5\times 10^5`$, $`L=170`$ and $`\nu _0=0.01`$ from here on, and vary only the selection strength through the choice of the selection stringency $`\varphi `$. A typical simulation result for strong selection ($`\varphi =0.1`$) is shown as the space-time plot of the mismatch distribution $`n(r,t)_\stackrel{}{S}n_\stackrel{}{S}\delta (E_\stackrel{}{S}rϵ)`$ in Fig. 1. We see that the distribution quickly forms a shape-preserving pulse (see the inset of Fig. 1), which moves, decelerates, and eventually reaches equilibrium in the neighborhood of the optimal sequence (at $`r=0`$). Basically, the selection eliminates weak binders in the population to improve the average binding energy, hence the selection threshold $`r_0`$ is decreased in the next round, while the change of $`r_0`$ further selects sequences with better binding energies. Along with new variety generated by mutation, a propagating pulse results.

We next investigate the dynamical behavior of the above evolution model analytically using a mean-field description. It will be convenient to describe the dynamics in the mismatch space $`r`$. Let us first consider the contribution from point mutation. For a sequence of length $`L`$, “alphabet size” $`𝒜`$ ($`𝒜=4`$ for nucleotides) and $`r`$ mismatches, there are $`L(𝒜1)`$ ways to mutate to a new sequence via a single point mutation. Among them, there are $`(Lr)(𝒜1)`$ ways to increase $`r`$ by one, and $`r`$ ways to reduce $`r`$ by one. Hence, a standard master equation can be written to describe the mutational dynamics of the distribution $`n(r,t)`$ in the mean field limit $`N1`$ REF:Peliti1997 , REF:Drossel2001 , REF:Gerland2001 . The effect of selection/amplification process can be phenomenologically modeled by an additive term proportional to $`\varphi ^1P(r,r_0)`$ for weak selection. Further taking the continuum limit in $`r`$ (valid in the limit of large $`L`$ and smooth population distribution), we arrive at the following mean-field description for $`n(r,t)`$ REF:Gerland2001 :

$`_tn=_r\left[_r(D(r)n)v(r)n\right]+U[r;n]n(r,t)`$ (1)
$`U[r;n]=\left[\varphi ^1P(r,r_0(t))1\right]/\tau ,`$ (2)
The first two terms on the right-hand side of (1) result from the (conservative) mutational processes, with

$$D(r)=\frac{\nu }{2}\left(1\frac{𝒜2}{𝒜1}\frac{r}{L}\right),v(r)=\nu \left(1\frac{𝒜}{𝒜1}\frac{r}{L}\right)$$
(3)
being the “diffusion coefficient” and “drift velocity” respectively REF:Gerland2001 , $`\nu \nu _0L`$. The $`r`$-dependences of $`v`$ and $`D`$ reflect the different phase space volume for the different mismatches. For example, the form of $`v(r)`$ ensures that the the distribution approaches the maximum entropy point with $`\overline{r}=\frac{𝒜1}{𝒜}L`$ mismatches by mutation alone. The third term in Eq. (1) represents the effect of the selection/amplification, controlled by the growth function $`U`$ defined in Eq. (2). (The factor $`\tau O(1)`$ denotes generation time.) Competition is explicitly manifested in the $`n`$ dependence of the growth function $`U`$, via the threshold $`r_0(t)`$ which is determined from the condition $`\varphi =𝑑rP(r,r_0(t))n(r,t)`$. In Eq. (2), an overall shift in $`U`$ by the constant $`1`$ has been included to ensure that the population size $`N`$ is conserved after selection/amplification, in accordance with the evolution process. This shift produces the desired competitive effect that individuals which bind better than the threshold $`r_0(t)`$ are reproduced and those not meeting the threshold decay away. In the actual analysis, we will approximate the Fermi function $`P(r,r_0)`$ by a step function $`\mathrm{\Theta }(r_0r)`$, which turns out not to affect the qualitative behavior.

We will see that the simplicity of the continuum mean-field equation provides an analytic understanding of generic features of the evolutionary dynamics, including the existence of the decelerating, shape-preserving population pulse; it also provides an analytical estimate of the smallness of the finite-$`N`$ correction. However that quantitative differences do exist between our simplified description and the breeding schemes employed in the simulations and experiments, due to the phenomenological nature of simplified description, and the continuum approximations used in both the mismatch space and time. (A quantitatively more accurate approach has recently been developed by Kloster and Tang REF:Kloster2 .)

We start with the simplest case of infinite sequence length (while keeping $`\nu `$ a finite constant), yielding constant coefficients $`D(r)=D`$ and $`v=\nu `$. Making the Ansatz in Eq. (1) that the distribution $`n(r,t)=n[y(r,t)]`$ where $`yrr_0(t)`$ and $`r_0(t)=ct`$ for some constant speed $`c`$, we obtain the ODE

$$Dn^{\prime \prime }(y)[(c+\nu )n(y)]^{}+u(y)n(y)=0,$$
(4)
where $`u(y)\left(\varphi ^1\mathrm{\Theta }(y)1\right)/\tau `$. A physically allowed solution of Eq. (4) exists for every $`cc_0\nu `$, where

$$c_0\sqrt{4D(\varphi ^11)/\tau }.$$

In fact, the smallest possible speed $`c_{\mathrm{min}}c_0\nu `$ is selected by the dynamics given a reasonably compact initial distribution. Here, velocity selection follows the familiar marginal stability mechanism REF:Collet1990 : The selected solution $`n^{}(y)`$ (with the velocity $`c_{\mathrm{min}}`$) is the one that decays most sharply at the pulse front (i.e., the $`r<r_0`$ end) among all the allowed solutions. Thus, as the front of the distribution broadens from the initial condition, it first reaches the asymptotic decay of $`n^{}`$. From then on, the distribution stops broadening and moves with the speed $`c_{\mathrm{min}}`$. Standard arguments show that this is equivalent to the condition that the front be marginally stable in the frame of reference moving with $`c_{\mathrm{min}}`$. Note that as $`c_0\sqrt{D}\sqrt{\nu }`$, $`c_{\mathrm{min}}<0`$ when the mutation rate $`\nu `$ is sufficiently large, indicating the worsening of the overall affinity of the sequences due to accumulation of deleterious mutations despite the presence of selection. These results apply also to the more general Fermi function, as it is only the asymptotic behavior of the growth term \[i.e., $`U(rr_0)`$\] that governs velocity selection.

In evolutionary dynamics, the population size $`N`$ often plays a very important role REF:Peliti1997 , REF:Drossel2001 . To see how $`N`$ enters, we note that the mean-field equation (1) has an inconsistency in that at the very front of the moving pulse, arbitrarily small $`n`$ gets the benefit of exponential amplification. But in reality, the number of individuals is discrete so that $`n`$ should always be greater than $`1/N`$. To deal with this problem, a cutoff procedure was proposed within the mean-field framework REF:Levine1996 , REF:Derrida1997 . Here we employ this procedure to estimate the effect of a finite population on the evolutionary velocity FNOTE:phasespace . Specifically, we modify the selection/amplification term $`u(y)`$ in Eq. (4) to $`u(y)\mathrm{\Theta }(nN^1)`$ (for $`y<0`$). A direct extension of the approach in Ref. REF:Derrida1997  leads to $`\delta c_0/c_0\frac{\pi ^2}{2}/\mathrm{ln}^2N`$ for the fractional change in $`c_0`$, which has the same scaling form as that for the Fisher equation REF:Derrida1997 . To test this result, we ran simulations (with a modified mutational scheme to achieve a constant drift) to measure the propagating speed of the pulse for different population sizes $`N`$. Our finding of $`\delta c/c0.06`$ between $`N=5\times 10^3`$ and $`5\times 10^8`$ is in line with the expectation and indicates that under typical experimental conditions, the fluctuation effect due to finite population is insignificant.

We next examine the more realistic situation of finite sequence length $`L`$. The important new effect is due to the $`r`$-dependence in the drift velocity $`v(r)`$ \[see Eq. (3)\], which, as the population approaches towards $`r=0`$, increasingly hinders its advance. This can already be appreciated if we assume a quasi steady-state dynamics and replace $`\nu `$ in the formula for $`c_{\mathrm{min}}`$ with $`v(r)=\nu \gamma r`$ \[where $`\gamma \frac{4}{3}\nu _0`$ according to (3)\]: We find a stable stationary position, $`\overline{r}^{\mathrm{EQ}}(\nu c_0)/\gamma `$ where $`c_{\mathrm{min}}=0`$. Here we identify this position naturally with the mean of the population $`\overline{r}^{\mathrm{EQ}}rn^{\mathrm{EQ}}(r)𝑑r`$.

To proceed with a more rigorous analysis, we neglect the $`r`$ dependence in $`D`$ which has only a small quantitative effect. Also, we assume that the equilibrium position $`\overline{r}^{\mathrm{EQ}}1`$ so that the boundary condition at $`r=0`$ can be safely ignored.

Returning to the mean-field equation (1), we use the same moving-pulse Ansatz as before except that we no longer fix a linear time dependence to the threshold $`r_0(t)`$. This Ansatz produces a linear ODE for $`r_0(t)`$:

$$\gamma r_0(t)+\dot{r}_0(t)=\gamma r_0^{\mathrm{EQ}}$$
(5)
where $`r_0^{\mathrm{EQ}}`$ is the equilibrium threshold, so that a static distribution can be achieved in the moving frame. The population mean $`\overline{r}`$ follows exactly the same motion (in fact, $`\overline{r}(t)r_0(t)`$ except when selection is very weak), i.e., a single exponential with time constant $`\gamma `$ (which depends only on the point mutation rate $`\nu _0`$). This is a generic result independent of the details of the fitness function, as long as a pulse solution exists for Eq. (1). The decay constant $`\gamma `$ obtained from simulation of the discrete model is in quantitative agreement with the expectation ($`\frac{4}{3}\nu _0`$) for weak selection ($`1>\varphi 0.25`$); an example is shown in Fig. 2(a). In fact, the same qualitative result (i.e., the existence of a shape-preserving pulse) holds for strong selection as shown already in Fig. 1 where $`\varphi =0.1`$.

The shape of the pulse, i.e. the equilibrium distribution $`n^{\mathrm{EQ}}`$, is governed by the same ODE as Eq. (4), except that the constant velocity $`c`$ is replaced by $`\gamma (y+r_0^{\mathrm{EQ}})`$. The resulting equation again has a continuum of physically allowed solutions, each having a different shape and corresponding to a different equilibrium position $`r_0^{\mathrm{EQ}}`$ (hence different $`\overline{r}^{\mathrm{EQ}}`$). Here we have an interesting generalization of velocity selection to the selection instead from a continuum of decelerating pulses. Again, starting from a compact initial distribution, the dynamics selects the solution REF:Solutions  whose front ($`y<0`$) decays most rapidly, (in this case a Gaussian falloff), whereas the other solutions have a power-law front.

The $`\overline{r}^{\mathrm{EQ}}`$ extracted from the selected solution agrees well with its heuristic approximation of $`(\nu c_0)/\gamma `$ when $`\gamma 1`$; see Fig. 2(b). The theory is quantitatively accurate FNOTE:Shape  when the selection is not too strong (e.g., $`\varphi ^1<2.5`$). For very strong selection, the equilibrium threshold position $`r_0^{\mathrm{EQ}}`$ approaches $`r=0`$ and the boundary condition there needs to be taken into account. When $`c_0`$ and $`\gamma `$ are expressed in terms of original parameters, $`\overline{r}^{\mathrm{EQ}}=(\nu c_0)/\gamma `$ suggests that for a population with sequence length $`L1`$, the population pulse could stall at $`\overline{r}^{\mathrm{EQ}}1`$. In order for the population to reach the optimal at $`r0`$, we need to increase the selection strength (i.e., lowering $`\varphi `$) or decrease the mutation rate so that $`\varphi ^11+\nu _0\tau L/2`$, to overcome the bigger entropic barrier associated with longer sequences.

As there have been extensive studies of evolution on various landscapes in the context of population genetics REF:Burger2000 , REF:Peliti1997 , REF:Drossel2001 , it is worth comparing the dynamical behavior of competitive evolution with that of more common evolutionary models. The traditional study of evolution focuses on fixed fitness landscapes, where every genotype (e.g., sequence $`\stackrel{}{S}`$) has a predetermined absolute fitness value (i.e., the reproductive rate of the sequence $`\stackrel{}{S}`$). Competitive evolution is different in that it is subject to a dynamic fitness landscape. That is, the fitness is measured relative to a dynamic selection threshold and progress towards the best binding sequence occurs via competition among the currently existing genotypes — being better is all-important, not being best. This aspect of competitive evolution leads to qualitatively different dynamical behavior. For comparison, we can consider the simplest and most widely studied fixed landscape, i.e., the smooth “Mt. Fuji” landscape REF:Burger2000 , REF:Peliti1997 , REF:Drossel2001 , REF:Woodcock1996 , where each nucleotide contributes independently and additively to the fitness of the sequence, thus forming a landscape on which fitness rises steadily toward a single peak. For infinite sequence length, the mean-field theory fails in that it produces an unphysical, run-away solution REF:Levine1996  due to the unlimited growth rate of $`n`$ at the high fitness states REF:Woodcock1996 , REF:Levine1996 , REF:Kessler1997andRidgeway1998 , and a finite population has a traveling speed that is essentially proportional to population size REF:Kessler1997andRidgeway1998 . For finite sequence length, the finite population dynamics is orders-of-magnitude slower in reaching equilibrium than the (now non-divergent) mean-field prediction REF:Kessler1997andRidgeway1998 . In contrast, finite population effects merely cause a small correction for the competitive evolutionary process.

To summarize, we investigated the dynamics of competitive evolution in the context of molecular evolution experiments. The major result concerns the existence and properties of a shape-invariant population pulse which propagates towards an eventual equilibrium configuration. Analytical results on the motion of the pulse obtained from the mean-field equation are in good agreement with simulations. Also, corrections due to finite population size are shown to be insignificant. An interesting aspect of our findings is the convergence of the evolution process to a solution far from optimal (i.e., $`\overline{r}^{\mathrm{EQ}}1`$), if the selection strength is not sufficiently strong or mutation rate not sufficiently low. In general, competitive evolution is rather different from the usual picture of climbing a fixed fitness landscape. This approach may be applicable more generally, e.g. to natural evolution in cases where competition for scarce resources is the primary driving force, as an organism only needs to be more efficient than its competitors to win the battle for evolutionary survival.

We thank M. Kloster for detail comments and discussions. We also acknowledge helpful discussions with D. Kessler, Q. Ouyang, L. Tang and J. Widom. This research is supported by the NSF through Grants DMR-0211308 and MCB-0083704.


format: text

text:
# Degree of entanglement for two qubits

## I Introduction

Quantum entanglement is the most surprising nonclassical property of composite quantum systems . As it is well-known, a qubit (or a spin-$`1/2`$ particle) is described by the $`2\times 2`$ density matrix $`\rho (𝐧)=(\mathrm{𝟏}+\stackrel{}{\sigma }𝐧)/2`$, $`|𝐧|1,`$ where $`\mathrm{𝟏}`$ is the unit matrix, $`\stackrel{}{\sigma }=(\sigma _1,\sigma _2,\sigma _3)`$ the Pauli matrices vector, and $`𝐧`$ the Bloch vector. $`|𝐧|=1`$ corresponds to a pure state, otherwise a mixed state. Whereas, an entangled pairs of two qubits is completely described by the following $`4\times 4`$ density matrix:

$$\rho _{AB}^{}=\frac{1}{4}(\mathrm{𝟏}\mathrm{𝟏}+\stackrel{}{\sigma }^A𝐮\mathrm{𝟏}+\mathrm{𝟏}\stackrel{}{\sigma }^B𝐯+\underset{i,j=1}{\overset{3}{}}\beta _{ij}\sigma _i^A\sigma _j^B),$$
(1)
from which one could obtain two reduced density matrices

$`\rho _A^{}=\mathrm{tr}_B(\rho _{AB}^{})={\displaystyle \frac{1}{2}}(\mathrm{𝟏}+\stackrel{}{\sigma }^A𝐮),`$ (2)
$`\rho _B^{}=\mathrm{tr}_A(\rho _{AB}^{})={\displaystyle \frac{1}{2}}(\mathrm{𝟏}+\stackrel{}{\sigma }^B𝐯),`$ (3)
for the two qubits A and B, where $`𝐮`$ and $`𝐯`$ are Bloch vectors for particles A and B, respectively; $`\beta _{ij}`$ are some real numbers.

It has been shown that entangled pairs are a more powerful resource than separable, i.e., disentangled, pairs in a number of applications, such as quantum cryptography , dense coding , teleportation  and investigations of quantum channels , communication protocols and computation  . The superior potentiality of entangled states has raised a natural question: “ How much are two particles entangled?”, since pairs with a high degree of entanglement should be a better resource than less entangled ones. Many measures of entanglement proposed in the past have relied on either the Schmidt decomposition  or decomposition in a magic basis . In an interesting paper, Abouraddy et al. devised a new measure of entanglement for pure bipartite states of two qubits, based on a decomposition of the state vector as a superposition of a maximally entangled state vector and an orthogonal factorizable one . Although there are many such decompositions, the weights of the two superposed states are remarkably unique. The square of the weight of the maximally entangled state vector (i.e., $`P_E=p^2`$) is then defined as the degree of entanglement for two qubits, such a measure is consistent with three measures of entanglement previously set forth: maximal violation of Bell’s inequality , concurrence  and two-particle visibility .

The purpose of this paper is to propose a new approach to the problem of defining the degree of entanglement for two qubits in a pure state. In Sec. II, a new measure is formulated to quantify the degree of entanglement. Some examples are given in Sec. III. Conclusion and discussion are made in the last section.

## II Formalism

Theorem: If $`\rho _{AB}^{}`$ is a pure state, then its degree of entanglement $`P_E`$ is equal to

$$P_E=(\mathrm{det}\widehat{\alpha })^{1/4}$$
(4)
where the matrix $`\widehat{\alpha }`$ is

$$\widehat{\alpha }=\left(\begin{array}{cccc}1\hfill & v_1\hfill & v_2\hfill & v_3\hfill \\ u_1\hfill & \beta _{11}\hfill & \beta _{12}\hfill & \beta _{13}\hfill \\ u_2\hfill & \beta _{21}\hfill & \beta _{22}\hfill & \beta _{23}\hfill \\ u_3\hfill & \beta _{31}\hfill & \beta _{32}\hfill & \beta _{33}\hfill \end{array}\right).$$
(5)
Proof: $`\rho _{AB}^{}`$ is a pure state implies that $`\rho _{AB}^2=\rho _{AB}^{}`$, from which one obtains the following constraints among $`u_i`$, $`v_i`$ and $`\beta _{ij}`$ $`(i,j=1,2,3)`$:

$$u_i=\beta _{i1}v_1+\beta _{i2}v_2+\beta _{i3}v_3,$$
(6)
$$v_i=\beta _{1i}u_1+\beta _{2i}u_2+\beta _{3i}u_3,$$
(7)
$$\underset{i,j}{}\beta _{ij}^2=3|𝐮|^2|𝐯|^2,$$
(8)
$$\beta _{ij}=u_iv_j(1)^{i+j}M_{ij},$$
(9)
where $`M_{ij}`$ is the algebraic complement of the matrix element $`\beta _{ij}`$ for the following $`\widehat{\beta }`$ matrix:

$$\widehat{\beta }=\left(\begin{array}{ccc}\beta _{11}\hfill & \beta _{12}\hfill & \beta _{13}\hfill \\ \beta _{21}\hfill & \beta _{22}\hfill & \beta _{23}\hfill \\ \beta _{31}\hfill & \beta _{32}\hfill & \beta _{33}\hfill \end{array}\right).$$
(10)
Eqs. (6) and (7) can be recast as $`\widehat{\beta }𝐯=𝐮`$, $`\widehat{\beta }^T𝐮=𝐯`$, where $`T`$ represents transpose and $`𝐮=(u_1,u_2,u_3)^T`$. An interesting result, i.e., $`|𝐮|=|𝐯|`$, will be obtained immediately from Eqs. (6) and (7) for the pure state $`\rho _{AB}^{}`$ . From Eq. (9) we have

$`\beta _{11}^2+\beta _{12}^2+\beta _{13}^2=\beta _{11}u_1v_1+\beta _{12}u_1v_2+\beta _{13}u_1v_3`$ (11)
$`[(1)^{1+1}\beta _{11}M_{11}+(1)^{1+2}\beta _{12}M_{12}+(1)^{1+3}\beta _{13}M_{13}].`$ (12)
Due to $`\mathrm{det}\widehat{\beta }=\beta _{11}M_{11}\beta _{12}M_{12}+\beta _{13}M_{13}`$ and Eq. (6), one obtains

$`\beta _{11}^2+\beta _{12}^2+\beta _{13}^2u_1^2=\mathrm{det}\widehat{\beta }.`$ (13)
Similarly,

$`\beta _{21}^2+\beta _{22}^2+\beta _{23}^2u_2^2=\mathrm{det}\widehat{\beta },`$ (14)
$`\beta _{31}^2+\beta _{32}^2+\beta _{33}^2u_3^2=\mathrm{det}\widehat{\beta }.`$ (15)
After combining Eqs. (8), (13), (15), and taking $`|𝐮|=|𝐯|`$ into account, one easily obtains $`\mathrm{det}\widehat{\beta }=1|𝐮|^2`$ . Consequently, we have

$`(\mathrm{det}\widehat{\alpha })^{1/4}=[(\mathrm{det}\widehat{\beta })(1|𝐮|^2)]^{1/4}=\sqrt{1|𝐮|^2}.`$ (16)
One can know from Ref.  that $`P_E=2\kappa _1\kappa _2`$, where $`\kappa _1`$ and $`\kappa _2`$ are the two coefficients in the Schmidt decomposition $`|\mathrm{\Psi }=\kappa _1|x_1,y_1+\kappa _2|x_2,y_2`$, $`\rho _{AB}^{}=|\mathrm{\Psi }\mathrm{\Psi }|`$, where $`\{|x_1,|x_2\}`$ and $`\{|y_1,|y_2\}`$ are orthogonal bases for the Hilbert spaces of particles $`A`$ and $`B`$, respectively. It is easy to prove that $`\kappa _1=\sqrt{(1+|𝐮|)/2}`$, $`\kappa _2=\sqrt{(1|𝐮|)/2}`$, which are square-roots of the two eigenvalues of the reduced matrix $`\rho _A^{}`$ or $`\rho _B^{}`$. Therefore we have $`P_E=(\mathrm{det}\widehat{\alpha })^{1/4}`$. This ends the proof.

## III Examples

Example 1. For the state $`|\mathrm{\Psi }=(|00+|01+|11)/\sqrt{3}`$, one obtains the density matrix

$`\rho _{AB}^{}=|\mathrm{\Psi }\mathrm{\Psi }|={\displaystyle \frac{1}{3}}\left(\begin{array}{cccc}1\hfill & 1\hfill & 0\hfill & 1\hfill \\ 1\hfill & 1\hfill & 0\hfill & 1\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 1\hfill & 1\hfill & 0\hfill & 1\hfill \end{array}\right)`$
with the Bloch vectors $`𝐮=(2/3,0,1/3)^T`$, $`𝐯=(2/3,0,1/3)^T`$, and the alpha matrix

$$\widehat{\alpha }=\frac{1}{3}\left(\begin{array}{cccc}3& 2& 0& 1\\ 2& 2& 0& 2\\ 0& 0& 2& 0\\ 1& 2& 0& 1\end{array}\right).$$
(18)
One can have $`P_E=2/3`$, which is consistent with the result in Ref. .

Example 2. For the state $`|\mathrm{\Psi }=[|00+2(|01+|11)]/3`$, the density matrix is

$`\rho _{AB}^{}=|\mathrm{\Psi }\mathrm{\Psi }|={\displaystyle \frac{1}{9}}\left(\begin{array}{cccc}1\hfill & 2\hfill & 0\hfill & 2\hfill \\ 2\hfill & 4\hfill & 0\hfill & 4\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 2\hfill & 4\hfill & 0\hfill & 4\hfill \end{array}\right)`$
with $`𝐮=(8/9,0,1/9)^T`$, $`𝐯=(4/9,0,7/9)^T`$, and the alpha matrix

$$\widehat{\alpha }=\frac{1}{9}\left(\begin{array}{cccc}9& 4& 0& 7\\ 8& 4& 0& 8\\ 0& 0& 4& 0\\ 1& 4& 0& 1\end{array}\right).$$
(20)
Hence the degree of entanglement is $`P_E=4/9`$.

Example 3. For the maximally entangled state $`|\mathrm{\Psi }=(|00+|11)/\sqrt{2}`$, one obtains the density matrix

$`\rho _{AB}^{}=|\mathrm{\Psi }\mathrm{\Psi }|={\displaystyle \frac{1}{2}}\left(\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 1\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 1\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right)`$
with the Bloch vectors $`𝐮=𝐯=(0,0,0)^T`$, and the alpha matrix

$$\widehat{\alpha }=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right).$$
(22)
Thus $`P_E=1`$ reaches the highest value.

Example 4. For the disentangled pure state $`\rho _{AB}^{}=\frac{1}{2}(\mathrm{𝟏}+\stackrel{}{\sigma }^A𝐮)\frac{1}{2}(\mathrm{𝟏}+\stackrel{}{\sigma }^B𝐯)`$, where $`|𝐮|=|𝐯|=1`$, we have the alpha matrix as

$$\widehat{\alpha }=\left(\begin{array}{cccc}1& v_1& v_2& v_3\\ u_1& u_1v_1& u_1v_2& u_1v_3\\ u_2& u_2v_1& u_2v_2& u_2v_3\\ u_3& u_3v_1& u_3v_2& u_3v_3\end{array}\right).$$
(23)
Obviously $`P_E=0`$ indicates that $`\rho _{AB}^{}`$ is disentangled.

## IV Conclusion and Discussion

In conclusion, we have presented a measure to quantify the degree of entanglement for two qubits in a pure state. We would like to make some discussion in the following:

(i) The similar idea developed in this paper could be generalized to quantify the degree of entanglement for two qu$`N`$its (i.e., $`N`$-state quantum systems, $`N=2`$ and $`N=3`$ correspond to a qubit and a qutrit, respectively)  in a pure state. For instance, the density matrix for two entangled qutrits could be written as

$$\rho _{AB}^{}=\frac{1}{9}(\mathrm{𝟏}\mathrm{𝟏}+\sqrt{3}\stackrel{}{\lambda }^A𝐮\mathrm{𝟏}+\sqrt{3}\mathrm{𝟏}\stackrel{}{\lambda }^B𝐯+\frac{3}{2}\underset{i,j=1}{\overset{8}{}}\beta _{ij}\lambda _i^A\lambda _j^B),$$
(24)
where $`\lambda _i`$ $`(i=1,2,\mathrm{},8)`$ are the eight hermitian generators of $`SU(3)`$ (namely the usual Gell-mann matrices). For the state of two entangled qutrits

$$|\mathrm{\Psi }=\frac{1}{3}(|00+|11+|22),$$
(25)
its corresponding density matrix is

$$\rho _{AB}^{}=\frac{1}{9}(\mathrm{𝟏}\mathrm{𝟏}+\frac{3}{2}\underset{i,j=1}{\overset{8}{}}\beta _{ij}\lambda _i^A\lambda _j^B),$$
(26)
with the non-zero coefficients $`\beta _{11}=\beta _{33}=\beta _{44}=\beta _{66}=\beta _{88}=1`$, $`\beta _{22}=\beta _{55}=\beta _{77}=1`$. The elements $`\beta _{ij}`$, 1, $`𝐮`$ and $`𝐯`$ form a $`9\times 9`$ matrix $`\widehat{\alpha }`$, it is easy to show that $`P_E=(\mathrm{det}\widehat{\alpha })^{1/4}=1`$, which indicates that the state $`|\mathrm{\Psi }`$ in Eq. (25) is just a maximally entangled state.

(ii) After making the parametrization $`𝐮=\widehat{𝐮}\mathrm{tanh}\varphi _𝐮`$, where $`\widehat{𝐮}=𝐮/|𝐮|`$, the density matrix of a qubit $`\rho (𝐮)=(\mathrm{𝟏}+\stackrel{}{\sigma }𝐮)/2`$ can be connected to the Lorentz boost matrix $`L(𝐮)=\mathrm{exp}(\phi _𝐮\stackrel{}{\sigma }\widehat{𝐮}/2)=\mathrm{𝟏}\mathrm{cosh}(\phi _𝐮/2)+\stackrel{}{\sigma }\widehat{𝐮}\mathrm{sinh}(\phi _𝐮/2)`$ as

$$\rho (𝐮)=\frac{L(𝐮)}{2\mathrm{cosh}\varphi _𝐮},\varphi _𝐮=\phi _𝐮/2.$$
(27)
Obviously, $`\rho (𝐮)`$ and $`L(𝐮)`$ are in one-to-one correspondence. For the former, the physical meaning of the vector $`𝐮`$ is the Bloch vector in quantum mechanics, while for the latter the relativistic velocity. Due to the rapidity $`\phi `$, i.e., the hyperbolic angle, special relativity can be formulated in terms of hyperbolic geometry. As a result, some physical quantities have been found to have geometric significance, such as the Thomas rotation angle corresponds to the defect of a hyperbolic triangle . After viewing the Bloch vector $`𝐮`$ as an analogous relativistic velocity, the Bures fidelity $`F(\rho _1,\rho _2)=[\mathrm{tr}\sqrt{\sqrt{\rho _1}\rho _2\sqrt{\rho _1}}]^2`$ was found to have a geometric interpretation in the framework of hyperbolic geometry . Similarly, with the aid of the parametrization $`𝐮=\widehat{𝐮}\mathrm{tanh}\varphi _𝐮`$, it is not difficult to find that the entanglement degree $`P_E=\sqrt{1|𝐮|^2}=1/\mathrm{cosh}\varphi _𝐮`$ for two qubits in a pure state is the reciprocal of the Lorentz factor  in the hyperbolic geometry. The extension of our approach to the mixed states of two entangled qubits will be discussed elsewhere.


text:
# Phenomenological model for the remanent magnetization of dilute quasi-one-dimensional antiferromagnets

## I Introduction

At low temperatures, quasi-one-dimensional magnetic systems exhibit a wealth of interesting behavior, such as dimensional crossover, Smith and Friedberg (1968); de Jonge et al. (1975); Wang (1997) random quantum paramagnetism, Nguyen et al. (1996) order-by-disorder phenomena, Oseroff et al. (1995); Azuma et al. (1997) and Griffiths phases, Fisher (1995); Young and Rieger (1996) which have motivated many experimental and theoretical investigations. In most of these systems, three-dimensional (3D) order is eventually induced by interchain interactions. Taking advantage of the many analytical results available for one-dimensional model systems, this situation has been modeled in a variety of ways. Most of the existent approaches are based on linear-chain approximations, Scalapino et al. (1975); Trudeau and Plumer (1995); Schulz (1996) which treat correlations along the chains in an exact way, while introducing interchain couplings via effective fields. These approximations have been quite successfully applied to pure systems, and have given rise to generalized Ginzburg-Landau theories, Scalapino et al. (1975); McKenzie (1995) which account for fluctuations. Also, they have been widely used to explain disorder effects, Imry et al. (1975); Hone et al. (1975); Schouten et al. (1980); Korenblit and Shender (1993); S.Eggert et al.  which are among the main topics of research on quasi-one-dimensional systems.

In the present work we consider a class of quasi-one-dimensional compounds, Paduan-Filho et al. (1998); Becerra et al. (2000) represented by (CH3NH3)MnCl3$``$2H2O (abbreviated MMC) and (CH3)2NH2MnCl3$``$2H2O (abbreviated DMMC), consisting of localized spin systems in which the Mn2+ ions (spin $`S=\frac{5}{2}`$) lie along the crystalline $`b`$ axis, forming chains, and are antiferromagnetically coupled to each other by an intrachain interaction $`J/k_B`$ around $`3`$ K. Magnetic susceptibility and specific heat measurements Simizu et al. (1984) indicate the onset of 3D long-range order at Néel temperatures $`T_N=4.12`$ K for MMC and $`T_N=6.36`$ K for DMMC, with the magnetic moments aligning along the $`a`$ axis. These temperatures are compatible with an interchain interaction $`\left|J_{}\right|\left|J\right|\times 10^2`$. The character of this interaction is not reported in the literature. However, owing to the behavior of the materials upon dilution with non-magnetic Cd2+ ions, it has been suggested (see below) that ferromagnetic interchain couplings are present. At temperatures above $`T10`$ K, susceptibility results are well described by a one-dimensional Heisenberg $`S=\frac{5}{2}`$ model, but at lower temperatures anisotropy effects (probably of dipolar origin) become relevant. Simizu et al. (1984) Calculations based on a classical anisotropic Heisenberg model, with parameters derived from experiments on DMMC, reinforce the importance of anisotropy. Schouten et al. (1980) In particular, the behavior of the correlation length along the chains is predicted to cross over from Heisenberg-like to Ising-like as the temperature is lowered.

Substitution of small amounts of non-magnetic Cd2+ for Mn2+ ions induces the appearance of a remanent magnetization Paduan-Filho et al. (1998); Becerra et al. (2000) below $`T_N`$, when samples are cooled in the presence of fields of a few Oe directed along the easy axis. This remanent magnetization is observed to vary linearly with temperature, except immediately below $`T_N`$, where demagnetization effects seem to be relevant. Paduan-Filho et al. (1998) Moreover, there is an excess parallel susceptibility, which is in general associated with the existence of uncompensated magnetic moments in odd-sized segments formed along the chains upon dilution. Dupas and Renard (1978) Apparently, the linear temperature dependence of the remanent magnetization is of universal character, as inferred from measurements Becerra et al. (2000) performed on DMMC doped with Cd2+ (non-magnetic) and Cu2+ ($`S=\frac{1}{2}`$). Experiments Carvalho et al. (2001) performed on similar compounds, CsMnCl3$``$2H2O (CMC) and CsMnBr3$``$2H2O (CMB) doped with Cu2+ (for which the signs of the interchain interactions are well known), revealed that a remanent magnetization appears in CMB, in which interchain couplings are ferromagnetic along one of the transversal directions and antiferromagnetic along the other; in contrast, no net magnetization is observable in CMC, where all interactions are antiferromagnetic. These experimental results, combined with the observation that some effective ferromagnetic coupling is expected in order to sustain a net remanent magnetization, have led to the idea that ferromagnetic interchain interactions should also be present in DMMC and MMC. Becerra et al. (2000) However, in the lack of experimental data, up to now no conclusive evidence on this point seems to be available.

In this paper we introduce and discuss a phenomenological model for the low-temperature magnetic behavior of those compounds. By virtue of the previously discussed anisotropy effects, we believe the qualitative aspects to be captured by a $`S=\frac{5}{2}`$ Ising model, which in the pure limit (and in the simplest case) is described by the Hamiltonian

$$=J\underset{𝒓}{}S_𝒓S_{𝒓+𝒃}+\underset{𝒓}{}\underset{𝜹}{}J_𝜹^{}S_𝒓S_{𝒓+𝜹},$$
(1)
where $`J>0`$, $`𝒓`$ is a lattice vector, $`𝒃`$ is the primitive vector along the crystalline $`b`$ direction, $`𝜹`$ is a vector connecting a site to its nearest neighbors in the $`ac`$ plane, $`J_𝜹^{}=J_{}>0`$ if $`𝜹`$ is parallel to the $`a`$ axis, and $`J_𝜹^{}=J_{}`$ if $`𝜹`$ is parallel to the $`c`$ axis. Our approach is based on a linear-chain approximation, which treats the intrachain couplings ($`J`$) exactly, while introducing the weak interchain interactions ($`J_{}J`$) via Curie-Weiss terms connecting all spins (in such a way that a staggered effective field results, combining both ferro- and antiferromagnetic interchain interactions in a cooperative manner). At very low temperatures the chains are antiferromagnetically ordered, with a characteristic two-sublattice structure. Upon dilution, a very long chain breaks into finite segments, and uncompensated magnetic moments appear at the ends of odd-sized segments. On phenomenological grounds, we assume these moments to correlate ferromagnetically, with their direction determined in the experiments by the cooling field. For each segment of spins, the partition function can be exactly calculated; the total free energy of the chain is obtained by summing the free energies of segments of all sizes, with proper weighting factors. This procedure is detailed in Sec. II.1. Then, in Sec. II.2, we include the Curie-Weiss terms and discuss the results of the approximation. We show that this approach reproduces the linear temperature dependence of the remanent magnetization and the existence of an excess susceptibility. The final section is devoted to a discussion and conclusions.

## II Phenomenological model

### II.1 Nearest-neighbor interactions

We initially consider an open segment of $`n`$ Ising spins with antiferromagnetic couplings and alternating fields, described by the Hamiltonian

$$_n=J\underset{j=1}{\overset{n1}{}}S_jS_{j+1}\underset{j=1}{\overset{n}{}}h_jS_jD\underset{j=1}{\overset{n}{}}S_j^2,$$
(2)
where $`J>0`$ and $`h_j=h_1`$ $`(h_2)`$ for odd (even) $`j`$; a crystal field $`D`$ is also introduced. The spin variables $`S_j`$ take the values $`\pm \frac{1}{2}`$, $`\pm \frac{3}{2}`$ and $`\pm \frac{5}{2}`$. The alternating fields are introduced to give room to a staggered effective field needed to describe long-range antiferromagnetic order in the presence of interchain interactions. According to the phenomenological hypothesis that there are uncompensated magnetic moments pointing in a preferred direction, determined by the cooling field, we assume that the spins at the ends of odd-sized segments are always under the action of a field $`h_1`$. When the field is removed, the moments would remain uncompensated due to pinning by the non-magnetic impurities. For even-sized segments, the particular choice of a field $`h_1`$ at $`j=1`$ is of no consequence, since in this case the partition functions are symmetric under inversion.

As we are considering finite values of $`n`$, we must treat separately the cases of odd- and even-sized segments. Using the transfer-matrix technique, we can write the partition functions as

$$𝒵_{n1}^{\mathrm{odd}}=𝒗_1\left|𝐓^{\frac{n2}{2}}\right|𝒗_1$$
(3)
and

$$𝒵_n^{\mathrm{even}}=𝒗_1\left|𝐓^{\frac{n2}{2}}𝐓_1\right|𝒗_2=𝒗_2\left|𝐓_2𝐓^{\frac{n2}{2}}\right|𝒗_1,$$
(4)
where $`n`$ is an even number, $`𝐓=𝐓_1𝐓_2`$, the elements of the $`6\times 6`$ matrices $`𝐓_1`$ and $`𝐓_2`$ are given by

$`T_1(S_i,S_j)`$ $`=e^{\beta JS_iS_j+\frac{1}{2}\beta \left(h_1S_i+h_2S_j\right)+\frac{1}{2}\beta D\left(S_i^2+S_j^2\right)},`$ (5)
$`T_2(S_i,S_j)`$ $`=T_1(S_j,S_i),`$ (6)
and the components of the vectors $`𝒗_1`$ and $`𝒗_2`$ are

$`v_1(S_j)`$ $`=`$ $`e^{\frac{1}{2}\beta \left(h_1S_j+DS_j^2\right)},`$ (7)
$`v_2(S_j)`$ $`=`$ $`e^{\frac{1}{2}\beta \left(h_2S_j+DS_j^2\right)}.`$ (8)
The free energies associated with odd- and even-sized segments are

$$F_{n1}^{\mathrm{odd}}=k_BT\mathrm{ln}𝒵_{n1}^{\mathrm{odd}}$$
(9)
and

$$F_n^{\mathrm{even}}=k_BT\mathrm{ln}𝒵_n^{\mathrm{even}}.$$
(10)
We now consider a very long chain, and assume that each of the $`N`$ sites is occupied by a spin with probability $`p`$. For $`0<p<1`$, the chain is composed of finite segments of spins separated by empty sites. In the $`N\mathrm{}`$ limit, the number of segments of size $`n`$ is $`NP(n)=N(1p)^2p^n`$. Assuming that each segment is described by the Hamiltonian in Eq. (2), the total free energy per spin is given by the infinite series

$$f_{\mathrm{nn}}(h_1,h_2,T)=\frac{1}{p}\underset{n\mathrm{even}}{}\left[P(n1)F_{n1}^{\mathrm{odd}}+P(n)F_n^{\mathrm{even}}\right].$$
(11)
For $`p<1`$, as $`nP(n)`$ becomes negligible for sufficiently large $`n`$, this infinite series can be truncated and readily evaluated numerically.

Let us denote by type 1 (type 2) those spins under the action of a field $`h_1`$ ($`h_2`$). The numbers $`N_1`$ and $`N_2`$ of spins of either type in a chain can be determined by noting that in a segment of size $`n`$ there are $`n/2`$ type-1 spins if $`n`$ is even and $`(n+1)/2`$ type-1 spins if $`n`$ is odd. Thus, the fractions of type-1 and type-2 spins are

$$\frac{N_1}{N}=\underset{n\mathrm{odd}}{}P(n)\frac{n+1}{2}+\underset{n\mathrm{even}}{}P(n)\frac{n}{2}=\frac{p}{1+p},$$
(12)
and

$$\frac{N_2}{N}=\underset{n\mathrm{odd}}{}P(n)\frac{n1}{2}+\underset{n\mathrm{even}}{}P(n)\frac{n}{2}=\frac{p^2}{1+p},$$
(13)
respectively. For $`p<1`$, the difference between the two fractions will obviously generate a non-zero magnetization at zero temperature.

### II.2 Linear-chain approximation

In order to mimic the weak interchain coupling in the real compounds, we now assume that, in addition to the nearest-neighbor couplings inside each segment, there are also ferromagnetic Curie-Weiss (CW) interactions connecting all spins in the chain. We further assume that the CW interactions between two type-1 or two type-2 spins have strength $`J_{\mathrm{cw}}/N`$, but that the CW interactions between a type-1 and a type-2 spin are weaker by a factor $`\gamma `$. This $`\gamma `$ factor is introduced to allow for off-plane interchain couplings; in the pure limit ($`p=1`$) the chains are expected to exhibit antiferromagnetic order, so that $`\gamma `$ must be smaller than unity. Upon dilution, we expect the antiferromagnetic arrangement to survive inside each segment, and in principle this could lead to a variation of $`\gamma `$ with the concentration $`p`$, since the magnetic arrangement in the planes perpendicular to the chains could be disturbed. In any case, our results suggest that $`\gamma `$ is very small, if not zero, in the compounds under consideration.

The contribution of the CW interactions to the total energy per spin is

$$\epsilon _{\mathrm{cw}}=pJ_{\mathrm{cw}}(m_1^2+2\gamma m_1m_2+m_2^2),$$
(14)
where $`m_1`$ ($`m_2`$) is the magnetization per magnetic ion due to spins of type 1 (type 2). Since $`\epsilon _{\mathrm{cw}}`$ depends only on the averages $`m_1`$ and $`m_2`$, and not on the detailed configuration of the spins, it is convenient to perform a change of variables. So, we now introduce the Helmholtz potential per spin $`a_{\mathrm{nn}}(m_1,m_2,T)`$, related to the nearest-neighbor interactions, and defined by the Legendre transform

$$a_{\mathrm{nn}}(m_1,m_2,T)=f_{\mathrm{nn}}(\stackrel{~}{h}_1,\stackrel{~}{h}_2,T)+m_1\stackrel{~}{h}_1+m_2\stackrel{~}{h}_2,$$
(15)
where $`\stackrel{~}{h}_1`$ and $`\stackrel{~}{h}_2`$ are effective fields and

$$m_1=\left(\frac{f_{\mathrm{nn}}}{\stackrel{~}{h}_1}\right)_{\stackrel{~}{h}_2,T}\mathrm{and}m_2=\left(\frac{f_{\mathrm{nn}}}{\stackrel{~}{h}_2}\right)_{\stackrel{~}{h}_1,T}.$$
(16)
For given values of $`m_1`$ and $`m_2`$ we can write the total Helmholtz potential,

$$a(m_1,m_2,T)=a_{\mathrm{nn}}(m_1,m_2,T)+\epsilon _{\mathrm{cw}},$$
(17)
from which we obtain the relation between the external magnetic fields $`h_1`$, $`h_2`$ and the effective fields,

$$h_1=\left(\frac{a}{m_1}\right)_{m_2,T}=\stackrel{~}{h}_12pJ_{\mathrm{cw}}(m_1+\gamma m_2),$$
(18)
and similarly

$$h_2=\left(\frac{a}{m_2}\right)_{m_1,T}=\stackrel{~}{h}_22pJ_{\mathrm{cw}}(\gamma m_1+m_2).$$
(19)
Comparing these last results (for $`\gamma =0`$) with the local field at a site $`𝒓`$ due to its $`q_{}`$ nearest-neighbors in adjacent chains, as given by the Hamiltonian in Eq. (1), we conclude that $`J_{\mathrm{cw}}`$ can be estimated as

$$J_{\mathrm{cw}}\frac{1}{2}pq_{}J_{}.$$
(20)
The thermodynamically stable magnetizations are those which minimize the free-energy functional

$`\mathrm{\Phi }(h_1,h_2,T;m_1,m_2)`$ $`=`$ $`a(m_1,m_2,T)m_1h_1m_2h_2`$ (21)
$`=`$ $`f_{\mathrm{nn}}(\stackrel{~}{h}_1,\stackrel{~}{h}_2,T)\epsilon _{\mathrm{cw}}.`$
For low temperatures and small ratios $`J_{\mathrm{cw}}/J`$, setting $`h_1=h_2=0,`$ the stable values of $`m_1`$ and $`m_2`$ have opposite signs. In the presence of dilution ($`p<1`$), since we have $`\left|m_1\right|\left|m_2\right|`$, the model predicts a remanent magnetization per lattice site, $`m_r`$, given by

$$m_r=p(m_1+m_2).$$
(22)
In the $`T0`$ limit, $`m_r`$ reaches a saturation value,

$$\underset{T0}{lim}m_r=\frac{N_1N_2}{N}S=\frac{p(1p)}{(1+p)}S,$$
(23)
where $`S=\frac{5}{2}`$. The zero-field differential susceptibility $`\chi _0`$ can be calculated by setting $`h_1=h_2=h`$ and taking the $`h0`$ limit,

$$\chi _0=\underset{h0}{lim}\frac{m_r}{h}.$$
(24)
The Néel temperature is obtained from the solution of the equation

$$\left\{\frac{^2\mathrm{\Phi }}{m_1^2}\frac{^2\mathrm{\Phi }}{m_2^2}\left(\frac{^2\mathrm{\Phi }}{m_1m_2}\right)^2\right\}_{m_1=m_2=0}=0,$$
(25)
in the absence of an external field.

In Fig. 1 we show the experimental data Becerra et al. (2000) for the temperature dependence of the remanent magnetization in DMMC doped with 4.5% Cd (this concentration was estimated from high-temperature fits to a Curie-Weiss law). We also show results of our calculations for the remanent magnetization with 4.5% dilution, $`J_{\mathrm{cw}}/J=1.5\times 10^2`$, $`\gamma =0`$, and $`D=0`$. We obtained the best fit for the linear portion of the curve by setting the theoretical Néel temperature to 1.14 times the experimental value (which amounts to fitting $`J`$). This is a reasonable procedure, since our calculations are of a mean-field nature, and thus we do not expect to obtain good quantitative results for the Néel temperature. Of course, the qualitative features of the calculations are not sensitive to small variations in the parameters; however, no strictly universal behavior (in the sense of data collapse) could be identified. We point out that setting the value of the crystal field to high positive values turns the system into a $`S=\frac{1}{2}`$ Ising model, and in this case the linear temperature dependence of $`m_r`$ could not be so well reproduced. Note that, in view of Eq. (20), the value of $`J_{\mathrm{cw}}/J`$ used in the fit is fully compatible with the estimated experimental value of $`J_{}/J`$ mentioned in the Introduction. The calculated ratio of the Néel temperatures of the diluted and pure systems is $`0.86`$, compared to the experimental estimate Becerra et al. (2000) of $`0.99`$ for the real material. From Eq. (23), the saturation value of $`m_r`$ for $`1\%`$ dilution is $`0.497\%`$ of the sublattice magnetization in the pure system, in excellent agreement with the experimental estimate Paduan-Filho et al. (1998) of $`0.5\%`$ for MMC with $`1\%`$ Cd.

In Fig. 2 we use the previous set of parameters to plot the calculated zero-field susceptibility $`\chi _0`$ both in the pure limit and for $`4.5\%`$ dilution. The broad maxima in the curves reflect short-range intrachain antiferromagnetic correlations, while the cusps (indicated in the figure by the arrows) correspond to the Néel temperatures ($`T_N`$). As it is clear in the inset, for the dilute case we observe other features at lower temperatures. The small maximum close to $`T=0`$ is due to isolated spins, whose sole energy scale is given by the very weak interchain couplings, while the neighboring shoulder is due to small odd-sized segments, whose end spins are uncompensated (even-sized segments give negligible contributions to $`\chi _0`$ at such low temperatures). This is in sharp contrast to the pure limit, in which the susceptibility vanishes exponentially for $`T<T_N`$.

It should be mentioned that the present approach is a generalization of that used by Slotte Slotte (1985) to investigate the dilute $`S=\frac{1}{2}`$ Ising chain with competing short- and long-range interactions. However, owing to the presence of competition, his approach does not contemplate the possibility of long-range antiferromagnetic order at finite temperatures, even in the pure limit.

## III Conclusions

We introduced a phenomenological model for the remanent magnetization ($`m_r`$) in a class of dilute quasi-one-dimensional antiferromagnets, composed of weakly interacting spin chains. The model assumes the existence of uncompensated spins at the ends of odd-sized segments formed along the chains upon dilution. These spins are supposed to remain ferromagnetically correlated after a cooling field is removed. By using a linear-chain approximation, in which interchain interactions are treated at a mean-field level, we were able to reproduce the linear temperature dependence of $`m_r`$ for a set of parameters compatible with experimental estimates.

Our linear-chain approximation is based on the assumption that, even upon dilution, each segment feels an staggered effective field. Of course, this assumption is subject to some restrictions. Depending on the impurity concentration $`1p`$, the existence of uncompensated moments pointing in a preferred direction could lead to a complete destabilization of the magnetic ordering perpendicular to the chains (this can be seen by considering the effect, in a particular chain, of two neighboring non-magnetic ions, which may invert the roles of the alternating sublattices). In this case, spins along the chains would feel the same interchain effective field, irrespective of their position. Actually, this would lead to a ferromagnetic transition (with a diverging susceptibility), and the long-range antiferromagnetic ordering would not be recovered even as $`p1`$. We have performed calculations near this limit, and have checked that the critical temperature depends linearly on $`1p`$, being thus too small compared with experimental findings. Moreover, it is not possible to reproduce the linear temperature dependence of the remanent magnetization. We conclude that, at least in the low impurity concentrations used here, for which the occurrence of two neighboring non-magnetic atoms in the same chain is a rare event, our approximation is reasonable.

There remains the topic of identifying the precise mechanism responsible for the persistence of ferromagnetic correlations between the uncompensated spins. Monte Carlo simulations based on the Hamiltonian in Eq. (2) could be used to verify whether it is sufficient or necessary to have both ferro- and antiferromagnetic interchain interactions present in order to give rise to a remanent magnetization in quasi-one-dimensional systems.

###### Acknowledgements.

We thank C. C. Becerra, A. Paduan-Filho and T. A. S. Haddad for fruitful discussions. This work was partially financed by the Brazilian agency CNPq and Fapesp.


text:
# References

Reply to Schlottmann and Zvyagin
X.-Y. Ge, M.D. Gould, J. Links and H.-Q. Zhou
Department of Mathematics,
The University of Queensland, 4072,
Australia

## Abstract

This is a reply to a comment by P. Schlottmann and A.A. Zvyagin.

In their comment , Schlottmann and Zvyagin raise several issues regarding the nature of integrable impurities in one-dimensional quantum lattice models, and claim to expose false statements in our recent work . In order to address these issues in a pedagogical manner, we feel that it is appropriate to discuss these questions in terms of models based on the $`gl(2|1)`$ invariant solution of the Yang-Baxter equation. However it is important from the outset to make it clear that the arguments we will present below are general and apply to other classes of models.

The first point that we would like to make is that it is claimed in  that there are two approaches to the algebraic Bethe ansatz. However, it appears to us that approach (i) described in item (2) of  is the co-ordinate Bethe ansatz and we are completely bemused why this should be referred to as an algebraic Bethe ansatz. In the co-ordinate Bethe ansatz approach one starts with a prescribed Hamiltonian and then solves the Schrödinger equation to get the two-particle scattering matrices and the particle-impurity scattering matrix. Together these scattering matrices form the monodromy matrix. It is impossible to infer the Hamiltonian from such a monodromy matrix. In this context we do not feel that item (1) of  answers our query about the existence of the impurity monodromy matrix in the algebraic approach (ii) of item (2) in . Hereafter we will focus our attention to this case.

The solution of the Yang-Baxter equation

$$R_{12}(uv)R_{13}(u)R_{23}(v)=R_{23}(v)R_{13}(u)R_{12}(uv)$$
(1)
associated with the Lie superalgebra $`gl(2|1)`$ is an operator $`R(u)\mathrm{End}(VV)`$, where $`V`$ is a three-dimensional $`_2`$-graded space with one bosonic and two fermionic degrees of freedom. Explicitly, this operator takes the form

$$R(u)=u.II+P$$
(2)
where $`P`$ is the $`_2`$-graded permutation operator. For the purposes of constructing integrable one-dimensional quantum systems on a closed lattice it is usual to introduce the the Yang-Baxter algebra with elements $`\widehat{A}_{ij}(u)`$, $`i,j=1,2,3`$. Putting these operators as the elements of a $`3\times 3`$ matrix $`A(u)`$, the algebraic relations amongst these operators are determined by the requirement

$$R_{12}(uv)A_1(u)A_2(v)=A_2(v)A_1(u)R_{12}(uv).$$
(3)
It is apparent from eq. (1) that the $`R`$-matrix (2) provides a representation of the Yang-Baxter algebra. The tensor product representation

$$T(u)=R_{0N}(u)\mathrm{}\mathrm{}R_{02}(u)R_{01}(u),$$
(4)
called the monodromy matrix, is also a representation of the Yang-Baxter algebra. The transfer matrix $`t(u)`$ is defined as the supertrace of the monodromy matrix, viz. $`t(u)=\mathrm{str}_0T(u)`$ and the Hamiltonian is the logarithmic derivative of the transfer matrix evaluated at $`u=0`$, viz. $`H=t^1(0)t^{}(0)`$, where the prime indicates the derivative with respect to the spectral parameter $`u`$. For the present case this yields the supersymmetric $`tJ`$ model as demonstrated in .

In order to incorporate an integrable impurity into the model, one looks for a different representation of $`A(u)`$, say $`L(u)`$, and then constructs the impurity monodromy matrix

$$T(u)=R_{0N}(u)\mathrm{}\mathrm{}R_{02}(u)L_{01}(u\theta )$$
(5)
with the transfer matrix and Hamiltonian defined as above. We agree with  that

* the position of the impurity Lax operator $`L(u\theta )`$ in the chain is inconsequential (here we put it in the first site for convenience)
* that the parameter $`\theta `$, which plays the role of a coupling parameter, can be chosen arbitrarily
* multiple impurity Lax operators can be placed in the monodromy matrix, each with an independent coupling $`\theta _i`$, and that in this case the position of the impurity Lax operators has no effect on the Bethe ansatz equations.

We also agree with  that when one extends this approach to the case of an open chain, which is achieved by the additional requirement that there is a solution to the reflection equations, this situation persists. Any impurity Lax operator $`L(u\theta )`$ for the closed chain is also an impurity Lax operator for the open chain, its position in the chain is immaterial regarding the Bethe ansatz solution, the coupling $`\theta `$ is arbitrary and the solution for the extension to multiple impurities with independent couplings is independent of the position of the impurities in the chain. The essential feature which seems to be completely missed by the authors of  is that the converse is not necessarily true! There exists classes of boundary impurities, which are obtained via solutions of the reflection equations, which have no analogue in the corresponding closed chain case. This is one of the primary results that we have tried to convey in  for the $`q`$-deformed $`tJ`$ model. The case of the supersymmetric $`tJ`$ model (with $`q=1`$) was studied in , which is relevant to our discussions here. Our results show the existence of boundary Kondo impurities for the supersymmetric $`tJ`$ model which have no analogue in the closed chain case. We will return to this point later.

Next we will look, in closer detail, at the types of impurity Lax operators $`L(u)`$ which exist for the present case of the supersymmetric $`tJ`$ model. We do not claim that the present classification is complete. One class takes the generic form

$$L(u)=u.II+\underset{i,j}{\overset{3}{}}(1)^{[i]}e_i^j\pi (E_j^i)$$
(6)
where $`e_i^j`$ is the matrix acting on $`V`$ with 1 in the $`(i,j)`$ position and zeroes elsewhere and $`\pi `$ is an arbitrary representation of the $`gl(2|1)`$ generators $`E_j^i`$ satisfying the commutation relations

$$[E_j^i,E_l^k]=\delta _j^kE_l^i(1)^{([i]+[j])([k]+[l])}\delta _l^iE_j^k.$$

Above, $`[i]`$ is the grading index so that $`[i]=0`$ if $`i`$ is a bosonic label and $`[i]=1`$ if $`i`$ is a fermionic label. Below we will adopt the convention $`[1]=[2]=1,[3]=0`$. Choosing the fundamental representation in eq. (6) for the $`gl(2|1)`$ generators with $`\theta =0`$ yields the $`R`$-matrix (2). The case of non-zero $`\theta `$ gives rise to an impurity Lax operator which was studied in . Other impurity Lax operators can be obtained by choosing different representations for $`gl(2|1)`$. For the case of the dual to the fundamental representation this yields the impurity model studied in . Choosing the four-dimensional representations yields the impurity model of . According to , these models belong to the class (ii) described in item (2), in which there are two levels of Bethe ansatz equations corresponding to charge and spin degrees of freedom. However, it is claimed in item (4) of  that “the effect of the impurity matrix in the charge sector (first level Bethe ansatz) changes the commutation relations in the spin sector (second level).” This statement is completely at odds with the explicit studies of the Bethe ansatz solutions conducted in  and is also not supported by the findings of  where the algebraic Bethe ansatz was performed for the $`gl(2|1)`$ invariant $`R`$-matrix in the most general context, with the only assumption being made the existence of a reference state on which the action of the monodromy matrix takes an upper triangular form. This raises a serious doubt over the validity of the Bethe ansatz solution presented in Appendix A of , which we have already discussed in the introduction of . Some concerns about the Bethe ansatz treatments of one of the authors of  have been raised by other researchers  and we believe that a similar lack of appreciation of the subtleties involved has caused the current misunderstanding regarding our work. Even worse, in (5) of  the authors attempt to counter our criticisms by claiming the operator $`\widehat{A}_{21}(u)`$ to be a “raising operator”. Appendix A of  claims that the operators $`\widehat{A}_{12}(u)`$, $`\widehat{A}_{13}(u)`$, $`\widehat{A}_{21}(u)`$ and $`\widehat{A}_{23}(u)`$ all vanish on the reference state. If $`\widehat{A}_{21}(u)`$ is a raising operator then $`\widehat{A}_{12}(u)`$, $`\widehat{A}_{13}(u)`$ and $`\widehat{A}_{23}(u)`$ are lowering operators. The fact that they each vanish on the reference state indicates that the reference state is a “lowest weight state”, rather than a “highest weight state”. Clearly, a proper treatment does not assume that the raising operator $`\widehat{A}_{21}(u)`$ vanishes on a lowest weight state. These assumptions were not made in the Bethe ansatz treatments of  ( is an exception where such an assumption is valid) and indeed  goes into a detailed discussion about this very issue in the general context.

We believe also that this assumption of  is inconsistent with the transfer matrix eigenvalue (A1) presented therein. Here we give our argument in detail. Let $`|\mathrm{\Psi }`$ denote the reference state for the algebraic Bethe ansatz procedure. It is assumed in  that the reference state is an eigenstate of the diagonal operators $`\widehat{A}_{ii}(u)`$, viz.

$$\widehat{A}_{ii}(u)|\mathrm{\Psi }=\lambda _i(u)|\mathrm{\Psi },$$
(7)
and

$$\widehat{A}_{21}(u)|\mathrm{\Psi }=\widehat{A}_{12}(u)|\mathrm{\Psi }=0.$$
(8)
The explicit commutation relations between the operators $`\widehat{A}_{ij}(u),i,j=1,2`$ can be conveniently expressed as

$$(uv)[\widehat{A}_{ij}(u),\widehat{A}_{kl}(v)]=\widehat{A}_{il}(u)\widehat{A}_{kj}(v)\widehat{A}_{il}(v)\widehat{A}_{kj}(u)$$

and in particular we have

$$(uv)[\widehat{A}_{12}(u),\widehat{A}_{21}(v)]=\widehat{A}_{11}(u)\widehat{A}_{22}(v)\widehat{A}_{11}(v)\widehat{A}_{22}(u).$$
(9)
Acting (9) on $`|\mathrm{\Psi }`$ immediately shows, through use of (7,8), that

$$\lambda _1(u)\lambda _2(v)=\lambda _1(v)\lambda _2(u)$$

for all values of the parameters $`u,v`$, which immediately implies that

$$\lambda _2(u)=c\lambda _1(u)$$

for some constant $`c`$ independent of $`u`$. The constant will be unity if there is spin reflection symmetry. The action of the transfer matrix $`t(u)`$ on $`|\mathrm{\Psi }`$ thus reads

$$t(u)|\mathrm{\Psi }=\left(\lambda _3(u)(1+c)\lambda _1(u)\right)|\mathrm{\Psi }$$

which is inconsistent with the result obtained by setting $`M=N=0`$ in (A1) of .

To further confuse the issue, item (5) of  mistakenly identifies the parameter $`l`$ introduced in the introduction of  as characterizing the symmetry of the impurity Lax operator. It is in fact a label for a class of atypical representations of the $`gl(2|1)`$ superalgebra. Applying such a representation to eq. (6) produces an impurity Lax operator which bears resemblence to that alluded to in Appendix B of . It should be clear in this instance why the operator $`\widehat{A}_{21}(u)`$ cannot vanish on the (lowest weight) reference state, since $`\widehat{A}_{21}(u)`$ is a spin raising operator, and the “higher spin” atypical representations of $`gl(2|1)`$, which are being used for the local impurity Hilbert space, do not admit a spin singlet. A description and discussion of the applicability of these representations is given in .

Item (6) of  provides a misrepresentation of our claims. We agree that it is “incorrect that impurities can only modify the first level Bethe ansatz equations.” The point that we make in the introduction of  is that if $`\widehat{A}_{12}(u)`$ and $`\widehat{A}_{21}(u)`$ both vanish on the reference state (as is assumed in ), then the second level Bethe ansatz equations will be unchanged (contradicting the final results presented in ).

We now return to the issue of impurity Lax operators. Another class can be obtained by taking the Lax operator associated with a subalgebra of the symmetry of the $`R`$-matrix. For example, in the present $`gl(2|1)`$ invariant case there is a natural $`gl(2)`$ subalgebra which has the Lax operator

$$L(u)=u.(e_1^1+e_2^2)I\underset{i,j}{\overset{2}{}}e_i^j\pi (F_j^i)$$
(10)
where the operators $`F_j^i`$ satisfy the $`gl(2)`$ commutation relations. This impurity operator can be used in (5) to construct an impurity model which would appear to describe a Kondo type impurity, except for one important facet. The operator (10) is singular, which in turn makes the transfer matrix singular, and does not yield a well defined Hamiltonian in terms of the logarithmic derivative of the transfer matrix. It is for this reason that we looked to the reflection equations for the study of Kondo impurities in the $`tJ`$ model in .

Finally we turn to the pertinent results of our works . In these papers we solve the reflection equations for the purpose of introducing Kondo impurities in the supersymmetric $`tJ`$ model and the $`q`$-deformed version, in a fashion which permits a well defined Hamiltonian. By this approach, the impurity interaction is purely of the Kondo type in the sense that there are no charge degrees of freedom, in spite of the protestations of  that this situation is not possible. The argument of  that charge degrees of freedom should be present for the impurity site, on the basis of symmetry, is wrong because it is possible to break symmetry at the boundary and maintain integrability, just as in the case of boundary scalar fields (cf.  for a discussion of this point in relation to integrable models with broken $`U_q(gl(n))`$ symmetry due to boundary interactions). In  the symmetry of the Hamiltonian is broken from $`gl(2|1)`$ to $`gl(2)`$ due to the inclusion of the boundary interactions.

Solutions of the reflection equation are specific to the construction of open chain models and there is no justification to expect that such integrable boundary impurities have an analogue in the bulk case. An attempt to do so for the model of  yields a singular open chain transfer matrix analogous to that discussed above for the periodic case, for which the Hamiltonian cannot be defined. We have argued vehemently in the conclusion of  that the boundary impurity $`K`$-matrices are operator valued and do not arise as the “dressing” of a scalar $`K`$-matrix. This indicates that these impurities are integrable only when situated on the boundary and do not have a bulk counterpart. (The situation really is no different to the case of integrable boundary scalar fields which also cannot be moved into the bulk. It is also true of model I in  studying a spin ladder model with defect in the rung coupling, where integrability holds only when the defect is on the boundary.) It is important to emphasize that our construction is entirely different from that discussed in  where the $`K`$-matrix is “the ordinary $`c`$-number reflection matrix of a free boundary sandwiched between two forward scattering impurity matrices.” This conclusion was confirmed by  who developed the “projection method” for the construction of such impurities in a general context. Application to the specific case of the $`tJ`$ model was studied in . Furthermore, the case of boundary impurities in the $`q`$-deformed $`tJ`$ model was independently studied in  and also confirms the validity of our approach.


text:
# 1 Introduction

## 1 Introduction

The recently reported experimental result on double beta decay ($`0\nu \beta \beta `$) shows the possible evidence for the non-zero Majorana mass of the electron neutrino in the range of $`m_{ee}=(0.050.86)`$ eV at 95% C.L. with the best value of $`m_{ee}=0.4`$ eV. There are certain important implications of this result. Firstly, it rules out all models predicting the Dirac neutrino masses, leaving only the option for the Majorana type of neutrinos. It also allows the lepton number violating processes such as leptogenesis in a natural way. Secondly, this result together with the earlier experimental data on the atmospheric and the solar neutrino oscillations, allows only the degenerate and the inverted hierarchical solutions for the three generation left-handed Majorana neutrinos.

In the above context, it is important to construct theoretical models which can predict the degenerate and inverted hierarchical patterns of the Majorana neutrino mass matrices within the framework of the grand unified theories(GUTs) with or without supersymmetry. In this short paper we attempt to generate the degenerate as well as the inverted hierarchical pattern of the left-handed Majorana neutrino mass matrices using the see-saw formula in a model independent way. This is, in fact, a continuation of our earlier work where the neutrino mixings are provided from the texture of the right-handed Majorana mass matrix $`M_{RR}`$, while keeping the Dirac neutrino mass matrix $`m_{LR}`$ in the diagonal form. We had taken the Dirac neutrino mass matrix $`m_{LR}`$ as either the charged lepton mass matrix ( $`m_{LR}`$ =$`tan\beta `$ $`m_l`$ referred to as case(i)) or the up-quark mass matrix ( $`m_{LR}=m_{up}`$ referred to as case(ii)). While referring to the earlier paper for details, the model successfully generated both the hierarchical and the inverted hierarchical (having opposite sign mass eigenvalues) neutrino mass matrices as a result of the proper choice of the parameters in texture of $`M_{RR}`$. In section 2 we present the generation of the degenerate as well as inverted hierarchical neutrino mass matrices using the seesaw formula, and their predictions on mass eigenvalues and mixing angles. Section 3 is devoted to summary and conclusion.

## 2 Neutrino mass matrices from See-saw formula

The left-handed Majorana neutrino mass matrix $`m_{LL}`$ is given by the celebrated see-saw formula,

$$m_{LL}=m_{LR}M_{RR}^1m_{LR}^T$$
(1)
where $`m_{LR}`$ is the Dirac neutrino mass matrix in the left-right (LR) convention. The leptonic (MNS) mixing matrix is now given by $`V_{MNS}=V_{\nu L}^{}`$ where $`m_{LL}^{diag}`$ = $`V_{\nu L}`$ $`m_{LL}`$ $`V_{\nu L}^T`$. Here both $`m_{LR}`$ and the charged lepton mass matrix $`m_l`$ are taken as diagonal, whereas the right-handed Majorana neutrino mass matrix $`M_{RR}`$ as non-diagonal. Using the see-saw formula (1) we generate both patterns of $`m_{LL}`$ viz., (I) nearly degenerate and (II) inverted hierarchical neutrino mass matrices with same sign mass eigenvalues. We conentrate here only on the cases which have bimaximal mixings with (+) entries in $`m_{LL}`$ listed in Table-I.

Table-I: Zeroth order neutrino mass matrices with texture zeros allowed by the result on $`0\nu \beta \beta `$ decay. The (+) entries in the matrices correspond to the LMA MSW solution with bimaximal mixings whereas (-) correspond to SMA MSW solution with single maximal mixing.

| Type | $`m_{LL}`$ | $`m_{LL}^{diag}`$ |
| --- | --- | --- |
| I(A) | $`\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& \pm 1\\ 0& \pm 1& 0\end{array}\right)m_0`$ | $`Diag(1,\pm 1,1)m_0`$ |
| I(B) | $`\left(\begin{array}{ccc}\pm 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)m_0`$ | $`Diag(\pm 1,1,1)m_0`$ |
| II(A) | $`\left(\begin{array}{ccc}1& 0& 0\\ 0& \pm \frac{1}{2}& \pm \frac{1}{2}\\ 0& \pm \frac{1}{2}& \pm \frac{1}{2}\end{array}\right)m_0`$ | $`Diag(1,\pm 1,0)m_0`$ |
| II(B) | $`\left(\begin{array}{ccc}0& 1& 1\\ 1& 0& 0\\ 1& 0& 0\end{array}\right)m_0`$ | $`Diag(1,1,0)m_0`$ |

I(A). Nearly degenerate mass matrix with opposite sign mass eigenvalues:

The most common texture for the nearly degenerate mass matrix $`m_{LL}`$ with opposite mass eigenvalues is represented by

$$m_{LL}=\left(\begin{array}{ccc}(12\delta _12\delta _2)& \delta _1& \delta _1\\ \delta _1& \delta _2& (1\delta _2)\\ \delta _1& (1\delta _2)& \delta _2\end{array}\right)m_0$$
(2)
where $`m_0`$ controls the overall magnitude of the masses of the neutrinos whereas $`\delta _1`$ and $`\delta _2`$ give the desired splittings for solar and atmospheric data. When $`\delta _1=\delta _2=0`$, Eq.(2) reduces to the zeroth order mass matrix of the Type I(A) in Table-I, with no splittings.

The diagonalisation of $`m_{LL}`$ in Eq.(2) leads to the following eigenvalues and mixing angles \[see Appendix\]:
$`m_{\nu _1}(12\delta _2(\sqrt{3}+1)\delta _1)m_0,`$
$`m_{\nu _2}(12\delta _2+(\sqrt{3}1)\delta _1)m_0,`$
$`m_{\nu _3}m_0,`$
$`\mathrm{sin}^22\theta _{12}=\frac{2}{3}`$, $`\mathrm{sin}^22\theta _{23}=1`$, $`\mathrm{sin}^22\theta _{13}=0.`$
For the choice of the values of the parameters $`m_0=0.4`$ eV, $`\delta _1=3.6\times 10^5`$, $`\delta _2=3.9\times 10^3`$, Eq.(2) leads to the following numerical predictions
Mixing angles:
$`sin^22\theta _{12}=0.67`$, $`sin^22\theta _{23}1.0`$, $`|V_{e3}|=1.5\times 10^{14}`$,

Mass eigenvalues:
$`m_{\nu _i}=(0.39684,0.396892,0.4)`$ eV, $`i=1,2,3`$; $`\mathrm{\Delta }m_{12}^2=4.13\times 10^5eV^2`$ and $`\mathrm{\Delta }m_{23}^2=2.48\times 10^3eV^2`$.

The prediction on solar mixing angle is consistent with the LMA MSW solution.

Case (i) where $`m_{LR}`$ =$`tan\beta `$ $`m_l`$

The degenerate mass matrix $`m_{LL}`$ in Eq.(2) is now generated through see-saw formula in Eq.(1) for the following choices of $`m_{LR}`$ and $`M_{RR}`$ respectively:

$$m_{LR}=tan\beta \left(\begin{array}{ccc}\lambda ^6& 0& 0\\ 0& \lambda ^2& 0\\ 0& 0& 1\end{array}\right)m_\tau ,$$
(3)
and

$$M_{RR}=\left(\begin{array}{ccc}(1+2\delta _1+2\delta _2)\lambda ^{12}& \delta _1\lambda ^8& \delta _1\lambda ^6\\ \delta _1\lambda ^8& \delta _2\lambda ^4& (1+\delta _2)\lambda ^2\\ \delta _1\lambda ^6& (1+\delta _2)\lambda ^2& \delta _2\end{array}\right)v_R$$
(4)
The inverse of $`M_{RR}`$ has a simple form for $`\delta _2>\delta _1`$ and $`\delta _1,\delta _2<<1`$,

$$M_{RR}^1=\left(\begin{array}{ccc}(12\delta _12\delta _2)\lambda ^{12}& \delta _1\lambda ^8& \delta _1\lambda ^6\\ \delta _1\lambda ^8& \delta _2\lambda ^4& (1\delta _2)\lambda ^2\\ \delta _1\lambda ^6& (1\delta _2)\lambda ^2& \delta _2\end{array}\right)(v_R^1)$$

The expression for $`m_0`$ in Eq.(2) is worked out as $`m_0=m_\tau ^2tan^2\beta /v_R`$. For input values of $`m_0=0.4eV`$, $`tan\beta =40`$, $`m_\tau =1.7`$ GeV, and $`\lambda =0.22`$, the see-saw scale is calculated as $`v_R10^{13}`$ GeV. This in turn gives the masses of the three right-handed Majorana neutrinos after the diagonalisation of $`M_{RR}`$:
$`|M_{RR}^{diag}|=(4.67\times 10^{11},1.296\times 10^5,5.06\times 10^4)`$ GeV.

Case(ii) where $`m_{LR}=m_{up}`$

The texture of $`m_{LL}`$ in Eq.(2) is again realised through the see-saw formula (1) with the following textures of $`m_{LR}`$ and $`M_{RR}`$:

$$m_{LR}=\left(\begin{array}{ccc}\lambda ^8& 0& 0\\ 0& \lambda ^4& 0\\ 0& 0& 1\end{array}\right)m_t,$$
(5)
and

$$M_{RR}=\left(\begin{array}{ccc}(1+2\delta _1+2\delta _2)\lambda ^{16}& \delta _1\lambda ^{12}& \delta _1\lambda ^8\\ \delta _1\lambda ^{12}& \delta _2\lambda ^8& (1+\delta _2)\lambda ^4\\ \delta _1\lambda ^8& (1+\delta _2)\lambda ^4& \delta _2\end{array}\right)v_R$$
(6)
We have $`m_0=m_t^2/v_R`$ in Eq.(2), and with the input values $`m_0=0.4`$ eV, $`m_t=200`$ GeV, we obtain $`v_R=10^{14}`$ GeV and the mass eigenvalues of the right-handed Majorana neutrinos: $`|M_{RR}^{diag}|=(1.105\times 10^{11},3.035\times 10^3,5.005\times 10^{11})`$ GeV.

I(B). Nearly degenerate mass matrix with same sign mass eigenvalues:

We propose another form of the nearly degenerate mass matrix,

$$m_{LL}=\left(\begin{array}{ccc}(12\delta _12\delta _2)& \delta _1& \delta _1\\ \delta _1& (1\delta _2)& \delta _2\\ \delta _1& \delta _2& (1\delta _2)\end{array}\right)m_0$$
(7)
The diagonalisation of $`m_{LL}`$ in Eq.(7) leads to \[see Appendix\]
$`m_{\nu _1}(12\delta _2(\sqrt{3}+1)\delta _1)m_0,`$
$`m_{\nu _2}(12\delta _2+(\sqrt{3}1)\delta _1)m_0,`$
$`m_{\nu _3}m_0`$,
$`\mathrm{sin}^22\theta _{12}=\frac{2}{3}`$, $`\mathrm{sin}^22\theta _{23}=1`$, $`\mathrm{sin}^22\theta _{13}=0`$.

The numerical solution leads to $`m_{\nu _i}=(0.39684,0.396892,0.4)`$ eV, $`i=1,2,3`$ for the same choices of the input values of $`\delta _{1,2}`$ and $`m_0`$ as in Eq.(2). Further, the predictions on the three mixing angles are the same as in Eq.(2). When $`\delta _1=\delta _2=0`$, it reduces to the type I(B) in the Table-I.

This form of the mass matrix(7) can be realised in the see-saw mechanism(1) using the following textures of $`m_{LR}`$ and $`M_{RR}`$:

Case(i) where $`m_{LR}`$ = $`\mathrm{tan}\beta `$ $`m_l`$

Here $`m_{LR}`$ is given in Eq.(3) and the right-handed neutrino mass matrix takes the form

$$M_{RR}=\left(\begin{array}{ccc}(1+2\delta _1+2\delta _2)\lambda ^{12}& \delta _1\lambda ^8& \delta _1\lambda ^6\\ \delta _1\lambda ^8& (1+\delta _2)\lambda ^4& \delta _2\lambda ^2\\ \delta _1\lambda ^6& \delta _2\lambda ^2& (1+\delta _2)\end{array}\right)v_R$$
(8)
The inverse of $`M_{RR}`$ has a simple form for $`\delta _2>\delta _1`$ and $`\delta _1,\delta _2<<1`$,

$$M_{RR}^1=\left(\begin{array}{ccc}(12\delta _12\delta _2)\lambda ^{12}& \delta _1\lambda ^8& \delta _1\lambda ^6\\ \delta _1\lambda ^8& (1\delta _2)\lambda ^4& \delta _2\lambda ^2\\ \delta _1\lambda ^6& \delta _2\lambda ^2& (1\delta _2)\end{array}\right)(v_R^1)$$

The expression for $`m_0`$ in Eq.(7) is again worked out as $`m_0=m_\tau ^2tan^2\beta /v_R`$.

Case(ii) where $`m_{LR}=m_{up}`$

The texture of $`m_{LR}`$ is given in Eq.(5), and the texture of $`M_{RR}`$ has the form

$$M_{RR}=\left(\begin{array}{ccc}(1+2\delta _1+2\delta _2)\lambda ^{16}& \delta _1\lambda ^{12}& \delta _1\lambda ^8\\ \delta _1\lambda ^{12}& (1+\delta _2)\lambda ^8& \delta _2\lambda ^4\\ \delta _1\lambda ^8& \delta _2\lambda ^4& (1+\delta _2)\end{array}\right)v_R$$
(9)
and $`m_0=m_t^2/v_R`$.

II(A). Inverted hierarchical mass matrix with same sign mass eigenvalues:

The most general form of the inverted hierarchical mass matrix $`m_{LL}`$ with the same sign mass eigenvalues can be expressed as

$$m_{LL}=\left(\begin{array}{ccc}(12ϵ)& ϵ& ϵ\\ ϵ& a& (a\eta )\\ ϵ& (a\eta )& a\end{array}\right)m_0^{}$$
(10)
where $`a=0.5`$ and $`m_0^{}`$ is the overall factor for the masses of the neutrinos. The parameters $`ϵ`$ and $`\eta `$ give the desired splittings for solar and atmospheric data.

The diagonalisation of $`m_{LL}`$ in Eq.(10) leads to the following eigenvalues and mixing angles \[see Appendix\]:
$`m_{\nu _1}(1(\sqrt{3}+1)ϵ\frac{\eta }{2}+\frac{\sqrt{\eta ϵ}}{6})m_0^{},`$
$`m_{\nu _2}(1+(\sqrt{3}1)ϵ\frac{\eta }{2}\frac{\sqrt{\eta ϵ}}{6})m_0^{},`$
$`m_{\nu _3}`$$``$ $`\eta `$ $`m_0^{},`$
and mixing angles:
$`\mathrm{sin}^22\theta _{12}=\frac{2}{3}`$, $`\mathrm{sin}^22\theta _{23}=1`$, $`\mathrm{sin}^22\theta _{13}=0.`$
When $`ϵ=\eta =0`$, Eq.(10) reduces to the zeroth order mass matrix of the type II(A) in Table-I, with no solar splitting. For solution of the LMA MSW solar data and atmospheric neutrino oscillation, we have the choice of the parameters $`m_0^{}=0.05`$ eV, $`ϵ=0.002`$ and $`\eta =0.0001`$ leading to the following predictions:

Mixing angles:
$`sin^22\theta _{12}=0.67`$, $`sin^22\theta _{23}1.0`$, $`|V_{e3}|=3.04\times 10^{13}`$,

Mass eigenvalues:
$`m_{\nu _i}=(0.05007,0.04973,0.000005)`$ eV, $`i=1,2,3`$; leading to $`\mathrm{\Delta }m_{12}^2=3.393\times 10^5eV^2`$ and $`\mathrm{\Delta }m_{23}^2=2.47\times 10^3eV^2`$.

Here the neutrino mass eigenvalues are of the same sign and this differs from the other inverted hierarchical mass matrix having the following form,

$$m_{LL}=\left(\begin{array}{ccc}ϵ& 1& 1\\ 1& \delta _1& \delta _2\\ 1& \delta _2& \delta _1\end{array}\right)m_0^{},ϵ,\delta _1,\delta _2<<1$$
(11)
which gives opposite sign mass eigenvalues. For $`\delta _1`$, $`\delta _2`$, $`ϵ`$ = 0, it leads to the type II(B) in Table-I. This structure has been successfully generated within this model, and the radiative correction has been found to be weak. We shall not address again this model here. Instead, we concentrate on the generation of the texture of $`m_{LL}`$ given in Eq.(10) from the see-saw mechanism (1).

Case(i) when $`m_{LR}`$ = $`tan\beta `$ $`m_l`$

The inverted hierarchical mass mtrix $`m_{LL}`$ in Eq.(10) can now be realised with the choice of $`m_{LR}`$ = $`tan\beta `$ $`m_l`$ given in Eq.(3) and $`M_{RR}`$ of the following form

$$M_{RR}=\left(\begin{array}{ccc}2a\eta (1+2ϵ)\lambda ^{12}& \eta ϵ\lambda ^8& \eta ϵ\lambda ^6\\ \eta ϵ\lambda ^8& a\lambda ^4& (a\eta )\lambda ^2\\ \eta ϵ\lambda ^6& (a\eta )\lambda ^2& a\end{array}\right)\frac{v_R}{2a\eta }$$
(12)
where $`M_{RR}^1`$ has a simple form for $`ϵ>\eta `$ and $`ϵ,\eta <<1`$,

$$M_{RR}^1=\left(\begin{array}{ccc}(12ϵ)\lambda ^{12}& ϵ\lambda ^8& ϵ\lambda ^6\\ ϵ\lambda ^8& a\lambda ^4& (a\eta )\lambda ^2\\ ϵ\lambda ^6& (a\eta )\lambda ^2& a\end{array}\right)(v_R^1)$$

The expression for $`m_0^{}`$ in Eq.(10) is given by $`m_0^{}=m_\tau ^2tan^2\beta /v_R`$. For input values of $`m_0^{}=0.05`$ eV, $`tan\beta =5`$, $`m_\tau =1.7`$ GeV, we obtain $`v_R=1.445\times 10^{12}`$ GeV which leads to $`|M_{RR}^{diag}|=(9.742\times 10^3,2.831\times 10^4,7.24\times 10^{15})`$ GeV.

Case(ii) when $`m_{LR}=m_{up}`$
The texture of $`m_{LL}`$ in Eq.(10) is again realised using $`m_{LR}=m_{up}`$ given in Eq.(5) and the texture of $`M_{RR}`$:

$$M_{RR}=\left(\begin{array}{ccc}2a\eta (1+2ϵ)\lambda ^{16}& \eta ϵ\lambda ^{12}& \eta ϵ\lambda ^8\\ \eta ϵ\lambda ^{12}& a\lambda ^8& (a\eta )\lambda ^4\\ \eta ϵ\lambda ^8& (a\eta )\lambda ^4& a\end{array}\right)\frac{v_R}{2a\eta }$$
(13)
which leads to $`m_0^{}=m_t^2/v_R`$ in Eq.(10). Using the input values $`m_0^{}=0.05eV`$, $`m_t=200`$GeV, we have $`v_R=8\times 10^{14}`$ GeV and $`|M_{RR}^{diag}|=(2.4\times 10^4,4\times 10^{18},2.4\times 10^9)`$ GeV where the mass of the heaviest right-handed Majorana neutrino lies above the GUT scale but below the Planck scale.

In the above calculation we have confined to the (+) entries of $`m_{LL}`$ in Table-I and these are consistent with the LMA MSW solution of the solar neutrino problem. The (-) entries in the texture of $`m_{LL}`$ in Table-I which generally lead to SMA MSW solar neutrino oscillation solution, can also be worked out in this model under the same framework with the proper choices of ($`\pm `$) signs before the elements of the $`M_{RR}`$.

A few comments on the stability condition under radiative corrections are in order. The nearly degenerate mass matrices $`m_{LL}`$ in Eqs.(2) and (7) are found to be unstable under radiative corection in minimal supersymmetric standard model. The inverted hierarchical mass matrix with the same mass eigenvalues given in Eq.(10) is also found to be unstable under radiative correction. However, the inverted hierarchical mass matrix given in Eq.(11) with opposite sign mass eigenvalues is stable under radiative corection .

## 3 Summary

In summary, we generate the textures of the nearly degenerate as well as the inverted hierarchical left-handed Majorana neutrino mass matrices from the see-saw formula using the diagonal form of the Dirac mass matrix and non-diagonal form of the right-handed Majorana neutrino mass matrix. The predictions on lepton mixing angles $`sin^22\theta _{12}0.67`$, $`sin^22\theta _{23}1.0`$ and $`|V_{e3}|0`$ are in excellent agreement with the experimental values. This is also true for the predictions of $`\mathrm{\Delta }m_{12}^2`$ and $`\mathrm{\Delta }m_{23}^2`$ which are necessary for the $`0\nu \beta \beta `$ decays, LMA MSW solar oscillation and atmospheric oscillation data. In all cases the masses of the right-handed Majorana neutrinos are above the weak scale. An interesting observation is that the prediction on the solar mixing angle is large but not maximal without any extra fine-tuning. Though the present analysis is a model independent one, it may be a useful guide for building models under the framework of grand unified theories with the extended flavour U(1) symmetry.

## Appendix

Diagonalisation of the mass matrix $`m_{LL}`$

Although the diagonalisation of $`m_{LL}`$ is trivial, we find it convenient to follow the general result given in Ref.. The neutrino mass matrix $`m_{LL}`$ of the general form

$$m_{LL}=\left(\begin{array}{ccc}a& b& c\\ b& d& e\\ c& e& f\end{array}\right)$$

where $`c=t_{23}b`$ and $`t_{23}=\mathrm{sin}\theta _{23}/\mathrm{cos}\theta _{23}`$, $`f=d+(t_{23}^1t_{23})e`$
can be diagonalised by $`V_{MNS}`$

$$V_{MNS}=\left(\begin{array}{ccc}\mathrm{cos}\theta _{12}& \mathrm{sin}\theta _{12}& 0\\ \mathrm{cos}\theta _{23}\mathrm{sin}\theta _{12}& \mathrm{cos}\theta _{23}\mathrm{sin}\theta _{12}& \mathrm{sin}\theta _{23}\\ \mathrm{sin}\theta _{23}\mathrm{sin}\theta _{12}& \mathrm{sin}\theta _{23}\mathrm{cos}\theta _{12}& \mathrm{cos}\theta _{23}\end{array}\right)$$

where we have taken $`\mathrm{sin}\theta _{13}=0`$. This mixing MNS matrix transforms $`|\nu _i>`$ with the masses ($`m_{\nu _1}`$, $`m_{\nu _2}`$,$`m_{\nu _3}`$) into $`|\nu _f>`$ via
$`|\nu _f>=V_{MNS}|\nu _i>`$, $`f=e,\mu ,\tau `$ and $`i=1,2,3.`$ The mass eigenvalues and $`\theta _{12}`$ are calculated as
$`m_{\nu 1}=a\frac{1}{2}\sqrt{\frac{b^2+c^2}{2}}(x+\eta \sqrt{x^2+8}),`$ $`m_{\nu 2}`$ = ($`\eta `$$``$$`\eta `$ in $`m_{\nu 1}`$),
$`m_{\nu 3}=d+t_{23}^2(da+x\sqrt{\frac{b^2+c^2}{2}})`$, $`\mathrm{sin}^22\theta _{12}=\frac{8}{8+x^2}`$
with $`x=(ad+t_{23}e)/(\sqrt{\frac{b^2+c^2}{2}})`$ where $`|m_{\nu 1}|<|m_{\nu 2}|`$ is always maintained by adjusting the sign of $`\eta (=\pm 1)`$.


text:
# Branching of negative streamers in free flight

## Abstract

We recently have shown that a negative streamer in a sufficiently high homogeneous field can branch spontaneously due to a Laplacian instability, rather than approach a stationary mode of propagation with fixed radius. In our previous simulations, the streamer started from a wide initial ionization seed on the cathode. We here demonstrate in improved simulations that a streamer emerging from a single electron branches in the same way. In fact, though the evolving streamer is much more narrow, it branches after an even shorter propagation distance.

Streamers are widely known phenomena in electric breakdown: long and rapidly growing channels of high ionization penetrating a region of low or vanishing ionization under the influence of a strong electric field . In experiments, streamers frequently are seen to branch. On a phenomenological level, the so-called dielectric break-down model (DBM) , a variation of diffusion limited aggregation , has been suggested as a stochastic process explaining the branching. Such models are inspired by earlier concepts of streamer propagation where the space charge was assumed to be smeared out over the full streamer head. Branching then would occur due to randomly distributed ionization avalanches around the streamer. Such stochastic phenomenological models on larger scale have not been related to particular discharge models.

However, simulations  have established that the space charge of the streamer is concentrated in a thin layer around the head rather than being smeared out over the full head. This leads to a different field distribution and a much faster propagation mode, and to a revival of the concept of an ideally conducting streamer formulated by Lozansky and Firsov . In particular, this charge distribution allows for a completely different “Laplacian” instability mechanism of streamers  that can operate even without any randomness quite like in viscous fingering, dendritic growth of solids into undercooled melts etc. It naturally emerges in the minimal streamer model (1) – (4). In a recent paper , we found that in contrast to previous expectations, this Laplacian instability mechanism can cause the spontaneous branching of streamers propagating in a strong homogeneous field. The system was identical to that of previous simulations  except that the field was twice as high. However, our simulation results were somewhat limited by numerical constraints. Because of the unexpected results  and some non-smooth structures in the figures, some researchers wondered whether the figures showed a numerical rather than a physical instability .

Therefore we here present new simulations with improved numerics, and we show more details of the evolution. In particular, we increased the number of grid points in the simulations from 1000$`\times `$1000 to 2000$`\times `$2000, which improves both the accuracy of the numerical output data and the quality of the contour plots, so the new plots are smooth. Furthermore, we changed the boundary condition on the cathode from homogeneous Neumann $`_z\sigma =0`$ to homogeneous Dirichlet $`\sigma =0`$ for the electron density $`\sigma `$. This means that while previously the electrons freely could flow from the metal of the electrode into the gas, this current is now suppressed. With these improvements, the choice of the initial conditions now ceases to be constrained by numerical considerations. We continue to use cylinder symmetry to calculate effectively on a 2-dimensional grid. Doing so, we keep relying on the analytical argument that the constraint of cylindrical geometry suppresses some instability modes, and that therefore a truly 3-dimensional system would become unstable earlier or at the same time as the system with symmetry constraint, but certainly not later.

The improvements of the numerical code allow now for the simulation of a streamer starting from a single electron on the cathode rather than from a wide initial seed as in . Our new figures show the streamer branch in free flight. The initial electron on the cathode creates first an avalanche, then a streamer, and finally the streamer splits. The situation resembles the historical experiments of Raether , except that the field is higher.

In detail, we investigate the minimal streamer model, i.e., a “fluid approximation” with local field-dependent impact ionization reaction in Townsend approximation  in a non-attaching and non-ionized gas. In dimensionless units , the model has the form:

$`_t\sigma \left(\sigma 𝐄+D\sigma \right)`$ $`=`$ $`\sigma f(|𝐄|),`$ (1)
$`_t\rho `$ $`=`$ $`\sigma f(|𝐄|),`$ (2)
$`\rho \sigma `$ $`=`$ $`𝐄,𝐄=\mathrm{\Phi },`$ (3)
$`f(|𝐄|)`$ $`=`$ $`|𝐄|e^{1/|𝐄|}.`$ (4)
Here $`\sigma `$ and $`\rho `$ are the densities of electrons and positive ions, $`𝐄`$ is the electric field, $`\mathrm{\Phi }`$ the electric potential and $`D`$ the dimensionless diffusion constant. The mobility of ions is neglected. The same model was investigated in . The translation to physical units depends on the type and the pressure $`p`$ of the gas. With an effective field-dependent impact ionization coefficient $`\alpha (E)=Ap\mathrm{exp}\left(Bp/|E|\right)`$ and electron mobility $`\mu _e=\overline{\mu }/p`$ as in , and with the parameter values for nitrogen as in , we get for the scales of length, time, field and particle density

$`l_0{\displaystyle \frac{2.3\mu \text{m}}{(p/1\mathrm{b}\mathrm{a}\mathrm{r})}},t_0{\displaystyle \frac{3\text{ps}}{(p/1\mathrm{b}\mathrm{a}\mathrm{r})}},E_0200{\displaystyle \frac{\text{kV}}{\text{cm}}}{\displaystyle \frac{p}{1\mathrm{b}\mathrm{a}\mathrm{r}}},`$ (5)
$`n_0{\displaystyle \frac{510^{14}}{\text{cm}^3}}\left({\displaystyle \frac{p}{1\mathrm{b}\mathrm{a}\mathrm{r}}}\right)^2210^5[N_2]{\displaystyle \frac{p}{1\mathrm{b}\mathrm{a}\mathrm{r}}},`$ (6)
where $`[N_2]`$ is the neutral gas particle density. The pressure dependent scaling relations are of interest for laboratory experiments, as well as for so-called sprite discharges  in the mesosphere above thunderclouds.

The simulation is performed with the same voltage of 1000 and electrode distance of 2000 as in , corresponding to a field of 100 kV/cm and an electrode distance of 4.6 mm for nitrogen at $`p=1`$ bar. As an initial condition, one electron was distributed over one cell of the numerical grid next to the cathode. Because of the diffusive broadening of the drifting avalanche, this “distribution” of the initial electron is permissible.

In the figures, we omit the avalanche phase and concentrate on the developed streamer and its branching at times $`t=500`$, 525, 550 and 575. The field is applied in the $`z`$ direction, and $`r`$ is the radial coordinate extending up to $`r=2000`$ to suppress artifacts from the lateral boundaries. The simulation shows a streamer not connected to the electrode. It has the conical shape of the initial avalanche created by a single electron in a homogeneous field as first found in Raether’s experiments .

Fig. 1 shows the electron density $`\sigma `$ within the whole streamer body at four snapshots of the evolution. Due to the no-current-boundary condition on the electrodes, the streamer is not connected to the cathode at $`z=0`$. The equidensity lines mark multiples of 0.15 in all snapshots. The electrodes at $`z=0`$ and $`z=2000`$ are outside the region shown. The qualitative evolution is the same as in our previous paper : while the streamer extends, it becomes wider, the electron density increases, and finally the instability develops.

Fig. 2 shows four different observables of the whole streamer at one time instant $`t=525`$. We plot the electron density $`\sigma `$, the ion density $`\rho `$, the space charge density $`\rho \sigma `$, and the electric potential $`\mathrm{\Phi }`$. In the propagating front at the tip of the streamer, an overshoot of drifting electrons creates a negatively charged layer, while the back of the streamer is depleted from electrons and positively charged by the essentially immobile ions left behind. Therefore the “ion streamer” is longer in the back, and the ion density has no overshoot in the propagating tip. The resulting space charge densities are shown in the third plot: a positive space charge with a maximum density of $`0.25`$ in the back, and a negative space charge with a minimum density of $`0.25`$ in the propagating tip. While the whole object is electrically neutral, the electric polarization leads to a suppression of the field in the interior. This can be read from the large distance between equipotential lines inside the ionized body visible in the last plot. The streamer is indeed approaching the Lozansky-Firsov limit of “ideal conductivity” .

Fig. 3 shows a zoom into the tip of the streamer at the same time steps as in Fig. 1. Again the electron densities $`\sigma `$ are plotted, and addionally the equipotential lines $`\mathrm{\Phi }`$. Already at time $`t=525`$, the onset of the instability can be concluded from the somewhat irregular shape of the level lines. The curvature of the tip decreases, the densities increase, the field in the interior decreases and the field ahead of the streamer increases. Consequently, the tip becomes unstable and splits.

We conclude with a few remarks:
1. It is quite remarkable that the finger of the new simulations branches after an even shorter travel distance than the one growing out of a wide, highly ionized seed in , though it is more narrow during its whole evolution.
2. When the initial electron is not placed on the electrode but somewhere in the gap, simulations show the same figures as in this paper, essentially unchanged except for a trivial shift in space.
3. The transition times from avalanche to streamer and then further from streamer into the instability depend quite nonlinearly on the externally applied field as we will discuss in a forthcoming paper; therefore the question whether branching also would occur in a lower field after a longer propagation distance, is open.
4. Finally, the discovery of sprite discharges  in the mesosphere gives a new impetus to the study of negative streamers in high fields, since a negative sprite discharge propagates into increasing heights, i.e., areas of decreasing pressure; the scaling laws (5) therefore imply that the effective field continues to grow along their path.

Acknowledgement: A.R. was supported by the Dutch Research School CPS, the Dutch Physics Funding Agency FOM, and CWI Amsterdam.