text:
# On discretizations of axisymmetric systems

## I Introduction

Axisymmetric systems are notorious for the difficulties they pose in numerical simulations when they are expressed in coordinates adapted to the symmetry. The problem is due to the singular nature of the coordinates on the axis, i.e., the set of fixed points for the symmetry transformations or, equivalently, the set where the infinitesimal generator of the symmetry – the Killing vector – vanishes. Adapted coordinates, i.e., an angle $`\varphi `$ along the Killing vector, a radius $`r`$ measuring ‘distance’ from the axis and a third coordinate $`z`$, yield coordinates $`(r,z)`$ on the space of orbits of the action of the symmetry group. The general theory of group actions Jänich (1968); Frauendiener (1987) yields the following facts. The orbit space contains two types of orbits111We ignore here the so called ‘exceptional orbits’ because we are only interested in the neighbourhood of the axis.: the ‘regular orbit’ is a circle around the axis. It has maximal dimension (one) and the union of all regular orbits is dense in the space of orbits. The other type of orbit is a fixed point which has dimension zero. The set of fixed points is a submanifold of codimension two, the ‘axis’. The drop in dimension of the orbits shows up in the topology of the space of orbits which acquires a boundary corresponding to the axis points and in a degeneracy of the adapted coordinates: since on the axis (for fixed $`z`$) all values of $`\varphi `$ address the same point the angular coordinate becomes irrelevant on the axis, so that the coordinate system degenerates there.

Invariant tensor fields transform in a rigid way (Lie transport) under the symmetry so that they are determined on an entire regular orbit once they are known on one of its points. Therefore, in order to determine the fields on the entire space it is enough to determine them on a hypersurface which is transversal to the orbits. Thus, the problem is reduced to a problem in a space with lower dimension. For invariant tensor fields at a fixed point the invariance implies a transformation law among the tensor components which must be interpreted as a regularity condition for the tensor field on the axis.

The advantage of using the adapted coordinates is, of course, the manifestation of this reduction. In these coordinates the components of the tensor fields do not depend on the angular coordinate along the Killing vector. Thus, there are only two essential coordinates left over and the problem reduces in complexity because an axisymmetric problem on some open 3-dimensional domain reduces to an equivalent 2D-problem on a domain with boundary. The boundary conditions are obtained from regularity conditions on the axis.

The disadvantage of using the adapted coordinates is that the reduced equations become formally singular on the axis. The regularity conditions on the axis guarantee that the individual terms in the equation have in fact a regular limit on the axis so there is not a real problem. Yet, numerically, one has to evaluate terms which are formally $`0/0`$ or $`\mathrm{}\mathrm{}`$ which does pose problems. Thus, it seems that one has the choice between a 3D system with regular, such as Cartesian, coordinates and regular equations or a 2D system with degenerate coordinates and singular equations.

In an attempt to combine the best of both worlds, Alcubierre et al. Alcubierre et al. (1999) proposed what they called the cartoon method. This scheme was designed to borrow from the singularity-free nature of full 3D Cartesian coordinates $`x`$, $`y`$, $`z`$ and to allow the treatment of axisymmetric problems without the memory constraints of a full 3D problem. The method uses the Killing transport equation along the orbits to find the field values on the grid points of a Cartesian grid from its values on the hypersurface $`y=0`$. These values are obtained by interpolation from the grid points lying in that hypersurface. For a second order method one needs to keep only three hypersurfaces $`y=const`$ in memory so that the problem scales quadratically with the number of points and not cubically.

The plan of the paper is as follows. In section II we discuss the cartoon method in more detail. In order to have a specific problem at hand we apply it to the axisymmetric linear wave equation in section III. We write this equation as a first order symmetric hyperbolic system. The reason for this apparent complication is that we have implemented this method for the axisymmetric conformal field equations which are a symmetric hyperbolic system and because the linear wave equation served as a model equation to test various implementations Frauendiener and Hein (2002). It turns out that the crucial ingredient to the method is the interpolation and we discuss various possibilities. Furthermore, we show that there is no real difference between the cartoon method using Cartesian coordinates and the formulation in adapted coordinates so that the issues which are relevant for the former apply to the latter as well. In section IV we analyse the stability properties of the cartoon method with different interpolation procedures as well as other possible discretizations when the method of lines is used for time evolution. We end with a brief discussion of the results.

## II The Cartoon method

In this section we present the cartoon method in more detail. We consider a 3-dimensional Euclidean space $`^3`$ with Cartesian coordinates $`𝐱=(x,y,z)`$ on which we define the action of the circle group $`SO(2)=G`$

$$G=\{g=\left(\begin{array}{cc}a& b\\ b& a\end{array}\right)a^2+b^2=1\}$$

in the usual way by

$$\begin{array}{cc}\hfill \mathrm{\Phi }:(G,^3)& ^3,\hfill \\ \hfill (g,𝐱)& \mathrm{\Phi }(g,𝐱)=(axby,bx+ay,z)=g𝐱.\hfill \end{array}$$
(1)
The infinitesimal generator of this action is the vector field $`\xi =x_yy_x`$ which vanishes on the axis $`x=y=0`$, the set of fixed points for the action. The action $`\mathrm{\Phi }`$ is the prototype for any axisymmetric system in the neighbourhood of the axis in the sense that for any of these systems one can find local coordinates in which the action takes the above form Jänich (1968); Frauendiener (1987). Therefore, we will focus here on flat space because we are mainly interested in the local properties in a neighbourhood of the axis. Curvature will not play an essential role.

Introducing adapted cylindrical coordinates $`(r,\varphi ,z)`$ the coordinate expression for the action is

$$(\left(\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right),(r,\varphi ,z))(r,\varphi +\theta ,z)$$

and the Killing vector is given by

$$\xi =_\varphi .$$

These coordinates are singular on the axis which can be seen e.g., from the fact that the Jacobi matrix of the transformation between cylindrical and Cartesian coordinates has vanishing determinant there.

A tensor field $`T`$ in $`^3`$ is axisymmetric if for any $`gG`$ it satisfies the condition

$$\mathrm{\Phi }_g^{}T=T$$
(2)
where $`\mathrm{\Phi }_g^{}`$ denotes the pull-back with $`\mathrm{\Phi }_g:^3^3,𝐱g𝐱`$. This condition expresses the invariance of $`T`$ with respect to the action of $`SO(2)`$. The infinitesimal version of this condition is the equation of Killing transport

$$_\xi T=0.$$
(3)
Suppose we want to solve a quasi-linear system of PDEs of the form

$$\dot{T}=F(_aT,T,t,𝐱)$$
(4)
where $`T`$ is some given time dependent axisymmetric tensor field on $`^3`$. Here, $`_a`$ denotes covariant differentiation. The usual way to treat this system is to express it in terms of the adapted coordinates and to use the invariance condition (3) which asserts, that the components of $`T`$ in the basis of adapted coordinates do not depend on $`\varphi `$. Thus, the essential coordinates are $`r`$ and $`z`$ and the problem is reduced to a 2-dimensional one. However, from the form of the covariant derivative operator, i.e., the Christoffel symbols of the metric expressed in cylindrical coordinates it is clear, that the equation will, in general, contain singular terms.

A different procedure is to employ any coordinate system which is regular in a neighbourhood of the axis and then to make explicit use of the invariance condition (2) or its infinitesimal version (3Thornburgh (2002). So we may choose Cartesian coordinates with the axis at $`x=y=0`$. Let $`T`$ be an arbitrary axisymmetric tensor field on $`^3`$. From the invariance condition (2) we can find its values at any point in $`^3`$ from its values on the half-plane $`=\{y=0,x0\}`$. To illustrate this let us first assume that $`T`$ is a scalar field. From the invariance condition we have

$$T(𝐱)=T(g𝐱)$$

for any rotation matrix $`g`$ and each point $`𝐱=(x,y,z)`$. Thus, in particular, choosing

$$g=\left(\begin{array}{cc}\frac{x}{\sqrt{x^2+y^2}}& \frac{y}{\sqrt{x^2+y^2}}\\ \frac{y}{\sqrt{x^2+y^2}}& \frac{x}{\sqrt{x^2+y^2}}\end{array}\right)$$

for $`x0`$, $`y0`$ we have $`\mathrm{\Phi }_g(𝐱)=(r,0,z)`$ with $`r=\sqrt{x^2+y^2}`$ so that

$$T(𝐱)=T(r,0,z).$$
(5)
Similarly, if $`T=T_1dx+T_2dy+T_3dz`$ is a 1-form, we have

$$\mathrm{\Phi }_g^{}T(𝐱)=T_1(g𝐱)d(\mathrm{\Phi }_g^{}x)+T_2(g𝐱)d(\mathrm{\Phi }_g^{}y)+T_3(g𝐱)d(\mathrm{\Phi }_g^{}z).$$

Inserting the explicit form of the matrix $`g`$ yields

$$\begin{array}{c}\mathrm{\Phi }_g^{}T(𝐱)=T_1(g𝐱)\left(adxbdy\right)+T_2(g𝐱)\left(bdx+ady\right)+T_3(g𝐱)dz\hfill \\ \hfill =\left(aT_1(g𝐱)+bT_2(g𝐱)\right)dx+\left(aT_2(g𝐱)bT_1(g𝐱)\right)dy+T_3(g𝐱)dz\end{array}$$

so that, in this case, the invariance condition implies the three equations

$$T_1(𝐱)=aT_1(g𝐱)+bT_2(g𝐱),$$

$$T_2(𝐱)=bT_1(g𝐱)+aT_2(g𝐱),$$

$$T_3(𝐱)=T_3(g𝐱).$$

Choosing, as before, $`a=x/\sqrt{x^2+y^2}`$ and $`b=y/\sqrt{x^2+y^2}`$, we get the relations

$$T_1(𝐱)=\frac{1}{r}\left(xT_1(r,0,z)yT_2(r,0,z)\right),$$

$$T_2(𝐱)=\frac{1}{r}\left(yT_1(r,0,z)+xT_2(r,0,z)\right),$$
(6)
$$T_3(𝐱)=T_3(r,0,z).$$

Similar relations hold for tensor fields with different type. This shows that we can get information about the values at all points outside the axis from values at the points on $``$ except for the axis. *On* the axis these relations imply regularity conditions for the tensor fields. E.g., in the case of a 1-form above we get $`T_1(0,0,z)=T_2(0,0,z)=0`$ which follows by taking limits from various directions towards the axis points.

How are we to use this analytical background in numerical applications? Following Alcubierre et al. (1999) we cover the half-plane $``$ by a regular Cartesian grid

$$x_i=i\mathrm{\Delta }x,z_k=k\mathrm{\Delta }z,\text{ with }i,k,i0$$

which we think of as being a part of a full 3D cartesian grid with points given by $`x_i,y_j,z_k`$ with $`y_0=0`$. Each component of a tensor field $`T`$ yields a grid function $`T_{ijk}`$. The axisymmetry implies that we should be able to determine the functions $`T_{ijk}`$ from their values on $``$ alone.

In order to treat equations like (4) above one needs to compute approximations to the spatial derivatives of the components. This is straightforward in the case of $`x`$\- and $`z`$-derivatives. However, for the $`y`$-derivatives we need to find the values of the grid functions at points outside $``$.

The situation is illustrated in Fig. 1. It shows the view from the positive $`z`$-axis down onto the $`(x,y)`$-plane. The black dots indicate grid points in $``$, the point on the left with the cross being the axis. The open circles indicate points outside $``$ with $`y=\pm \mathrm{\Delta }y`$. The values of the fields at these points have to be computed from the invariance equation (2). Following the orbits through these outside points one can connect them with points in $``$, indicated by the open squares, and one can find the values at the outside points from the values at these points on $``$. Thus, e.g., applying (5) in this situation we obtain for a scalar

$$T(x_i,\pm \mathrm{\Delta }y,z_k)=T(\sqrt{x_i^2+\mathrm{\Delta }y^2},0,z_k),$$

while the corresponding transformation for a 1-form follows from (6). Note, that the squares are on $``$ but they are not, in general, grid points. Thus, having reduced the problem to $``$ the remaining question is how to determine the values at these points. The natural choice is to use some form of interpolation as it is also suggested in Alcubierre et al. (1999). However, it is not clear what kind of interpolation should be used and it is the purpose of this work to point out that not all interpolation schemes perform well.

## III The axisymmetric wave equation

To have a concrete example at hand we will focus now on the axisymmetric scalar wave equation in flat space. It is clear from the considerations above that the $`z`$-dependence is irrelevant for this analysis. Therefore, we simplify even further by assuming that the field is in fact axisymmetric and $`z`$-independent. Then the equation is

$$\ddot{\varphi }=\varphi _{xx}+\varphi _{yy}$$

where $`\varphi =\varphi (t,x,y)`$ is axisymmetric. Introducing the 1-form $`d\varphi =udt+vdx+wdy`$ we can derive from this equation a symmetric hyperbolic first order system

$$\dot{u}=v_x+w_y,$$
(7)
$$\dot{v}=u_x,$$
(8)
$$\dot{w}=u_y$$
(9)
with the constraint $`v_yw_x=0`$. The function $`\varphi `$ can be recovered from a solution of this system by either evolving $`\dot{\varphi }=u`$ or by solving $`\mathrm{\Delta }\varphi =v_x+w_y`$ at each time step. Since $`\varphi (t,x,y)=\varphi (t,\sqrt{x^2+y^2},0)=\varphi (t,\pm x,\pm y)`$ we have, in particular, $`w(t,x,0)=0`$ and $`v(t,0,0)=0`$. Furthermore, taking derivatives we obtain $`u_y(t,x,0)=0`$ and $`v_y(t,x,0)=0`$. Hence, we can forget about (9) and the constraint, because these equations are identically satisfied on $``$. So we are left with the system

$$\dot{u}=v_x+w_y,$$
(10)
$$\dot{v}=u_x.$$
(11)
Before we go into the more numerical details it is worthwhile to write down the system in adapted (i.e., cylindrical) coordinates $`(r,\varphi ,z)`$

$$\dot{u}=v_r+\frac{v}{r},\dot{v}=u_r,$$
(12)
which follows by expressing the Laplace operator in these coordinates, using the independence of $`z`$ and $`\varphi `$ and then introducing the first derivatives as new dependent variables.

We solve these equations numerically by applying the method of lines, i.e., we discretize the spatial derivatives to obtain a system of ODEs. This can then be solved with standard ODE solvers. We put the equations on a grid as indicated above. Now $``$ is in effect a 1-dimensional grid with points $`x_i=i\mathrm{\Delta }x`$ and we have grid functions $`u_i(t)=u(t,i\mathrm{\Delta }x,0)`$, $`v_i(t)=v(t,i\mathrm{\Delta }x,0)`$. According to the idea of the cartoon method we compute the term $`w_y`$, which of course does not vanish on $``$, by following the orbits. From the transformation law (6) for a 1-form we get

$$w(t,x,y)=\frac{y}{\sqrt{x^2+y^2}}v(t,\sqrt{x^2+y^2},0).$$
(13)
The discretization of the term $`w_y`$ using centered differences to get a second order accurate approximation yields

$$w_y(t,x_i,0)=\frac{w(t,x_i,\mathrm{\Delta }y)w(t,x_i,\mathrm{\Delta }y)}{2\mathrm{\Delta }y}=\frac{v_{}}{x_{}},$$

where $`x_{}=\sqrt{x_i^2+\mathrm{\Delta }y^2}`$ and $`v_{}=v(t,x_{},0)`$. This, then, yields the final system obtained by application of the cartoon method

$$\dot{u}=v_x+\frac{v_{}}{x_{}},\dot{v}=u_x.$$
(14)
It is instructive to compare this system with the system (12) obtained in adapted coordinates. Expanding in powers of $`\mathrm{\Delta }y`$ we get

$$\frac{v_{}}{x_{}}=\frac{v(x_{})}{\sqrt{x_i^2+\mathrm{\Delta }y^2}}=\frac{v_i}{x_i}+\frac{1}{2}\left(\frac{v}{x}(x_i)\frac{v_i}{x_i}\right)\frac{\mathrm{\Delta }y^2}{x_i^2}+𝒪(\mathrm{\Delta }y^4).$$

Hence, *up to the accuracy* $`\mathrm{\Delta }y^2`$ for which the derivation of (14) is valid, *the two systems agree*.

One can, in fact, say even more. The infinitesimal invariance condition (3) or, equivalently, taking the $`y`$-derivative of (13) and evaluating at $`y=0`$ yields

$$w_y(t,x,0)=\frac{v(t,x,0)}{x}$$

for all $`x>0`$. Hence, whenever we use the cartoon method to approximate $`w_y`$ to a certain order in the grid spacing $`\mathrm{\Delta }y`$ we approximate $`v/x`$ to the same order. Thus, there is no difference in that order of accuracy between the systems (14) and (12). And, therefore, the discussions in the following section apply to both.

## IV Numerical issues

In this section we discuss the stability properties of various implementations of the systems (12) and (14). In the spirit of what was said in the previous section we regard the term $`\frac{v_{}}{x_{}}`$ as an approximation to $`\frac{v}{x}`$ which is not necessarily located on a grid point. Then, the main question is how to approximate this term on $``$. We implement several possibilities and look at their stability properties.

In order to have a well defined problem we have to worry about boundary conditions at the origin $`x=0`$ and at an outer boundary which we put at $`x=1`$. The conditions at the origin follow from the transformation properties of the fields under the symmetry. They are

$$v(t,x,0)=v(t,x,0)v(t,0,0)=0,$$
(15)
$$u(t,x,0)=u(t,x,0)\frac{u}{x}(t,0,0)=0.$$
(16)
The fact, that $`x=0`$ is the fixed point of the symmetry is reflected in the simultaneous vanishing of $`v`$ and $`x`$. This also implies that at the origin we have $`\frac{v}{x}=v_x`$ so that at $`x=0`$ the time evolution of $`u_0=u(t,0,0)`$ is given by the equation $`\dot{u}_0=2v_x(t,0,0)`$.

Since we are interested in the properties of the scheme in the interior we implement some sort of periodic boundary conditions. Since we want to avoid a second origin at $`x=1`$ we imagine that our wave equation lives on a 2-sphere and is symmetric under rotations of the sphere around its poles. Now we impose the condition that it be also symmetric under reflection across the equator. This discrete symmetry allows us to consider the equation only along a meridian from the north pole down to the equator. The symmetry implies conditions for the fields at the equator which we put at $`x=1`$:

$$v(t,1x,0)=v(t,1+x,0)v(t,1,0)=0,$$
(17)
$$u(t,1x,0)=u(t,1+x,0)\frac{u}{x}(t,1,0)=0.$$
(18)
In the context of the method of lines we write down a system of ODEs for the grid functions $`v_i`$ and $`u_i`$. This will be a linear system of the form $`\dot{f}=𝐀f`$ and we can find its stability properties by looking at the spectrum of the matrix $`𝐀`$ which we determine numerically.

With these preparations out of the way we can now write down the first of our discretizations for the system (12). This is the first one which comes to mind, namely

$$\begin{array}{cc}\hfill \dot{u}_0& =\frac{2v_1}{h},\hfill \\ \hfill \dot{u}_i& =\frac{v_{i+1}v_{i1}}{2h}+\frac{v_i}{ih},\hfill \\ \hfill \dot{u}_N& =\frac{v_{N1}}{h}\hfill \end{array}\begin{array}{cc}\hfill \dot{v}_0& =0,\hfill \\ \hfill \dot{v}_i& =\frac{u_{i+1}u_{i1}}{2h}\hfill \\ \hfill \dot{v}_N& =0.\hfill \end{array}$$
(19)
Here we have chosen $`h=1/N`$ where $`N+1`$ is the number of grid points so that $`x_i=ih`$. The matrix $`𝐀`$ which corresponds to this discretization has dimensions $`2(N+1)\times 2(N+1)`$ and the form

$$𝐀=\frac{1}{h}\left(\begin{array}{cc}0& 𝐂\\ & \\ 𝐃& 0\end{array}\right)$$

where $`𝐂`$ is the matrix

$$𝐂=\left(\begin{array}{cccccc}0& 2& 0& \mathrm{}& & \\ \frac{1}{2}& 1& \frac{1}{2}& & & \\ & \frac{1}{2}& \frac{1}{2}& \frac{1}{2}& & \\ & & \mathrm{}& \mathrm{}& \mathrm{}& \\ & & & \frac{1}{2}& \frac{1}{N1}& \frac{1}{2}\\ & & \mathrm{}& 0& 1& 0\end{array}\right)$$

and $`𝐃`$ is

$$𝐃=\left(\begin{array}{ccccc}0& 0& \mathrm{}& 0& 0\\ \frac{1}{2}& 0& \frac{1}{2}& & \\ & \mathrm{}& \mathrm{}& \mathrm{}& \\ & & \frac{1}{2}& 0& \frac{1}{2}\\ 0& 0& \mathrm{}& 0& 0\end{array}\right).$$

This system of ODEs is solved by a standard ODE solver, such as a Runge-Kutta scheme or a multistep scheme like the Adams-Bashforth schemes. We also consider the ‘iterated Crank-Nicholson’ (ICN) scheme which has become popular among numerical relativists (see Teukolsky (2000) for an analysis of this scheme). Relevant for deciding about the stability of a difference approximation for ODEs is the region of absolute stability (see e.g. Hairer et al. (2000); Hairer and Wanner (1996)). This is a closed set $`\mathrm{\Omega }`$ in the complex plane which consists of those $`z=\lambda \mathrm{\Delta }t`$ for which the method applied to the equation $`\dot{u}=\lambda u`$ with step size $`\mathrm{\Delta }t`$ produces bounded approximations. The stability regions for the schemes mentioned above are shown in appendix A and the stability region(s) for ICN is determined in appendix B.

Thus, in the context of the method of lines, the stability conditions can be obtained by determining the eigenvalues of the discretization matrix $`𝐀`$ and checking whether it is possible to scale it so that it fits entirely into the stability region of the ODE solver. The scale factor determines the maximal admissible time step. Strictly speaking, this consideration applies only to the linear case. For non-linear equations this check should be understood as a ‘rule of thumb’ Trefethen (2000).

The spectrum $`\sigma _N(h𝐀)`$ of $`h𝐀`$ is shown for $`N=100`$ in Fig. 2.

We point out some of its properties which are valid for all other discretizations below as well. The spectrum is symmetric with respect to reflections about the real axis as well as about the imaginary axis. While the former symmetry is due to the reality of the equation, the latter is a consequence of the fact that the underlying equation is a second order wave equation which has in- and outgoing modes. Furthermore, the spectrum is concentrated mostly on the imaginary axis with a few ‘outliers’. It can be seen that the eigenvalues with the largest, in absolute value, real part remain unchanged when we change $`N`$ which implies that the spectral radius of $`𝐀`$ is proportional to $`N`$.

An interesting feature of the spectrum are the two ‘handles’ which are located near $`\pm \mathrm{i}`$. These, and the outliers are due to the $`v/x`$ term in the equation and, depending on the particular discretization scheme, they may or may not be present as we will see below. Since the eigenvalues of the corresponding continuous system are given by $`\pm \mathrm{i}\lambda _k`$, where $`\lambda _k`$ is the k-th zero of the Bessel function $`j_1`$, all the eigenvalues with nonvanishing real part are spurious, i.e., they do not have an analogue in the continuous problem. Thus, the corresponding eigenvectors are parasitic modes because they do not approximate a regular mode.

When combined with a time integrator this spectrum and the stability region of the integrator decide about the stability of the overall scheme. The stability regions of some ODE solvers are sketched in appendices A and B. Some of these include an interval on the imaginary axis while others include only the origin. All the stability regions extend further into the left half plane than into the right. This implies that if a spectrum contains points in the right plane (i.e., modes with positive real part) then these modes tend to either unnecessarily diminish the time step because they have to be scaled into the stability region (if the stability region of ODE solver extends into the right half plane int he first place) or else they create instabilities.

Another possibility to overcome instabilities is the addition of numerical dissipation. This has the effect that the spectrum of $`𝐀`$ is shifted towards the left. However, the presence of the handles and of the outliers implies that the amount of dissipation necessary for stabilising the scheme is rather high compared to cases where the spectrum is located entirely on the imaginary axis. This means that the evolution cannot be followed accurately because ultimately the waves will be damped out. Thus, the appearance of handles and outliers seems to leave us with the choice between the Scylla of instability and Charybdis of inaccuracy. Therefore, we take the appearance of these modes as a sign for an inadequate scheme.

The discretizations which we have looked at can be divided into three groups. Discretizations of the system (14) obtained from the cartoon method of computing the transversal derivatives, 2nd order discretizations of the system (12) and three somewhat non-standard discretizations of that system. In the list given below we only show the formula for the inner points of the grid. The values at the boundary points are obtained by using the appropriate symmetry conditions. Furthermore, we omit the equation for $`\dot{v}_i`$ in those cases where its discretization is the same as for the derivative term in the $`\dot{u}_i`$ equation.

1. Using the cartoon approximation we discretize the system (14) to second order in $`h`$. The derivatives are taken to be centered in all the three different schemes, while we use different interpolations for the term $`\frac{v_{}}{x_{}}`$.

1. Two point interpolation (see Fig. 3):

    $$\dot{u}_i=\frac{v_{i+1}v_{i1}}{2h}+\frac{(x_{i+1}x_{})v_i+(x_{}x_i)v_{i+1}}{hx_{}}.$$
    2. Three point interpolation (see Fig. 3):

    $$\dot{u}_i=\frac{v_{i+1}v_{i1}}{2h}+\underset{j=1}{\overset{1}{}}p_{i+j}(x_{})v_{i+j},$$

    where $`p_{i1}(x)=(xx_i)(xx_{i+1})/2h^2`$, $`p_i(x)=(xx_{i+1})(xx_{i1})/h^2`$ and $`p_{i+1}(x)=(xx_i)(xx_{i1})/2h^2`$.
    3. Centered two point interpolation (see Fig. 4):

    $$\dot{u}_i=\frac{v_{i+1}v_{i1}}{2h}+\frac{(x_{i+1}x_{})v_{i1}+(x_{}x_{i1})v_{i+1}}{2hx_{}}.$$

It is obvious that the interpolation in schemes $`1a`$ and $`1b`$ is not useful because of the appearance of the handles and the outliers. This suggests that for higher order Lagrange interpolation one has to expect the same kind of behaviour. In fact, in our studies related to Frauendiener and Hein (2002) we find that e.g. 5th order Lagrange interpolation is as unstable as interpolation using Chebyshev polynomials. However, the fact that scheme $`1c`$ yields a purely imaginary spectrum (up to round-off error) suggests that the use of the centered interpolation is advantageous. In view of the fact that the cartoon method provides an approximation of the $`1/x`$ term in (12) we now look at this system directly.
2. Three second order discretizations of the system (12).

1. The scheme discussed already above (see Fig. 2)

    $$\dot{u}_i=\frac{v_{i+1}v_{i1}}{2h}+\frac{v_i}{ih}.$$
    2. A centered approximation for the numerator of the second term (see Fig. 5)

    $$\dot{u}_i=\frac{v_{i+1}v_{i1}}{2h}+\frac{v_{i+1}+v_{i1}}{2ih}.$$
    3. A staggered grid discretization based on the variables $`u_i=u(t,ih)`$, $`v_{i+1/2}=v(t,ih+h/2)`$ (see Fig. 6)

    $$\dot{u}_i=\frac{v_{i+1/2}v_{i1/2}}{h}+\frac{v_{i+1/2}+v_{i1/2}}{2ih},\dot{v}_{i+1/2}=\frac{u_{i+1}u_i}{h}.$$

Again, we find that discretizing in a centered and compatible way yields a purely imaginary spectrum (again up to round-off error). The staggering in scheme $`2c`$ is suggested by the structure of the equations which is ultimately related to the meaning of the variables. While $`u`$ is a time derivative located at the grid points $`v`$ is a spatial derivative which should be located between the grid points.
3. The final three discretization schemes of (12) are somewhat different from the ones above:

1. This scheme takes the idea of centering to fourth order and again yields a purely imaginary spectrum (see Fig. 7)

    $$\dot{u}_i=\frac{v_{i2}8v_{i1}+8v_{i+1}v_{i+2}}{12h}+\frac{v_{i2}+4v_{i1}+4v_{i+1}v_{i+2}}{6ih}.$$
    2. The next scheme makes use of further information coming from the theory of group actions. While the radial coordinate $`r`$ is *not* a smooth coordinate on the orbit space its square is smooth. Hence, it seems natural to express the variables as functions of the coordinate $`x=r^2`$. Then the minor redefinition $`\stackrel{~}{u}(x)=2u(r)`$ and $`\stackrel{~}{v}(x)=rv(r)`$ of the variables yields the (regular) system

    $$_t\stackrel{~}{u}=_x\stackrel{~}{v},_t\stackrel{~}{v}=x_x\stackrel{~}{u},$$

    which is discretized in the standard way (see Fig. 8)

    $$\dot{\stackrel{~}{u}}_i=\frac{\stackrel{~}{v}_{i+1}\stackrel{~}{v}_{i1}}{2h},\dot{\stackrel{~}{v}}_i=\frac{i}{2}(\stackrel{~}{u}_{i+1}\stackrel{~}{u}_{i1}).$$
    3. The same system discretized using a staggered grid (see Fig. 9)

    $$\dot{\stackrel{~}{u}}_i=\frac{\stackrel{~}{v}_{i+1/2}\stackrel{~}{v}_{i1/2}}{h},\dot{\stackrel{~}{v}}_{i+1/2}=(i+1/2)(\stackrel{~}{u}_i\stackrel{~}{u}_{i1}).$$

## V Conclusion

We discussed in this work several schemes for treating axisymmetric systems which are difficult to handle numerically due to the degenarcy of adapted coordinates and hence the formal singularity of the ensuing equations. We applied various schemes to the axisymmetric wave equation written as a symmetric hyperbolic system for the derivatives and determined the spectrum of the ensuing matrices. As a consequence of the shape of the stability regions of common ODE solvers we take the appearance of handles and outliers described in section IV as an indication for a ‘bad’ discretization scheme. These eigenvalues correspond to grid modes which, when evolved with a standard ODE solver, will lead to instabilities. In order to avoid them one needs to be careful in discretizing the equations because the derivative term and the $`1/x`$ term must be treated in a compatible way. This is due to their common geometric origin. In the general case, where tensors of higher rank are involved different terms may have to be combined. Exactly which terms these are is indicated by the behaviour of the tensor components under the symmetry transformation.

One particular method we discussed is the cartoon method. It turns out that this method is not really different from a direct implementation of the axisymmetry in terms of the adapted coordinates. Here the problem is to make the interpolation compatible with the discretization in the $`y=0`$ half plane so that parasitic modes are excluded. However, the advantage in using the cartoon method is that it admits a quick and easy extension of a 2D axisymmetric code to a full 3D code because the underlying grid structure and the equations need not be changed.

Clearly, this investigation cannot deliver the ultimate solution to the problem. It should be regarded as providing a guide to eliminate ‘bad’ discretizations and suggesting possible ‘good’ ones. For instance, we have implemented the schemes $`1c`$ and $`3a`$ in our axisymmetric code Frauendiener and Hein (2002) for solving the conformal field equations. While these schemes perform well for the axisymmetric wave equation they still tend to produce slowly growing instabilities in the case of the quasilinear conformal field equations. This behaviour might be related to the fact that in most cases the real parts of the eigenvalues vanish only up to round-off error and that these small real parts start to grow due to non-linear interactions. It would be useful to have a scheme in which the spectrum is *exactly* on the imaginary axis and whose only non-vanishing eigenvalues are the physical ones. However, such a scheme is probably impossible to find.

###### Acknowledgements.

This work was supported in part by NATO Collaborative Linkage grant PST.CLG.978726.

## Appendix A Stability regions

We illustrate the procedure for determining the stability region of a particular time evolution scheme for an explicit Runge-Kutta scheme. Let $`\dot{u}=f(t,u)`$ be the (system of) ODE(s) which we intend to solve. An explicit Runge-Kutta scheme with $`s`$ ‘stages’ to compute an approximation $`u_1`$ of $`u(t_0+h)`$ is formulated as follows:

$`k_1`$ $`=f(t_0,u_0),`$
$`k_2`$ $`=f(t_0+c_2h,u_0+ha_{21}k_1),`$
$`k_3`$ $`=f(t_0+c_3h,u_0+ha_{31}k_1+ha_{32}k_2),`$
$`\mathrm{}`$
$`k_s`$ $`=f(t_0+c_sh,u_0+ha_{s1}k_1+ha_{s2}k_2+\mathrm{}+ha_{ss1}k_{s1}),`$
$`u_1`$ $`=u_0+h(b_1k_1+\mathrm{}+b_sk_s).`$
The particular scheme is characterized entirely by the coefficients $`a_{ik}`$, $`b_i`$ and $`c_k`$ which are conveniently presented as a tableau

$$\begin{array}{cccccc}0& & & & & \\ c_2& a_{21}& & & & \\ c_3& a_{31}& a_{32}& & & \\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & \\ c_s& a_{s1}& a_{s2}& \mathrm{}& a_{ss1}& \\ & & & & & \\ & b_1& b_2& \mathrm{}& b_{s1}& b_s\end{array}$$

The method is said to be of $`p`$-th order if for sufficiently smooth $`f`$ the condition

$$u_1u(t_0+h)Ch^{p+1}$$

holds.

To get to the stability region one applies the scheme to the specific model equation $`\dot{u}=\lambda u`$ where $`\lambda `$ is arbitrary. Going through the procedure above it is clear that the final result is of the form

$$u_1=P(\lambda h)u_0$$

where $`P(z)`$ is a polynomial of degree $`s`$, the so called stability function of the method. Since $`u_1=|P(\lambda h)|u_0`$ we find that the method will produce bounded approximations if and only if $`|P(\lambda h)|1`$. This motivates the definition of the stability region as the set

$$\mathrm{\Omega }=\{z||P(z)|1\}.$$

In case of an implicit scheme the stability function, obtained in the same way by applying the scheme to the model equation, is not a polynomial but a rational function. In a similar way, the stability regions of linear multistep schemes are obtained. Here, the simplest one is the leapfrog scheme which updates $`u_{n+1}`$ from the values on the previous time levels $`u_n`$ and $`u_{n1}`$ as $`u_{n+1}=u_{n1}+2hf(t_0+nh,u_n)`$. It is easy to see that the stability region for the leapfrog scheme is the part of the imaginary axis between $`\mathrm{i}`$ and $`\mathrm{i}`$.

The stability regions for explicit Runge-Kutta schemes of order $`p=s`$ for $`s=1,2,3,4`$ are shown in Fig. 10.

The stability regions for the second and third order Adams-Bashforth methods are shown in Fig. 11

as an example for stability regions of multistep methods. For the stability regions of other multistep schemes we refer to Fornberg (1996); Gustafsson et al. (1995). Qualitatively, these regions extend into the left half plane of the complex plane including possibly some part of the imaginary axis. For higher order methods the stability region may even extend into the right half plane.

## Appendix B The ICN time evolution scheme

The ICN scheme was described in Teukolsky (2000) in application to a one-dimensional advection equation. Viewing the ICN scheme as a time stepper to solve systems of ODEs $`\dot{u}=f(t,u)`$ one can deduce from this description the following algorithm to compute $`u_1u(t_0+h)`$

$`k_1`$ $`=f(t_0,u_0),`$
$`k_2`$ $`=f(t_0+\frac{1}{2}h,u_0+\frac{1}{2}hk_1),`$
$`k_3`$ $`=f(t_0+\frac{1}{2}h,u_0+\frac{1}{2}hk_2),`$
$`\mathrm{}`$
$`k_s`$ $`=f(t_0+\frac{1}{2}h,u_0+\frac{1}{2}hk_{s1}),`$
$`u_1`$ $`=u_0+hk_s.`$
Cleary, this is an explicit Runge-Kutta scheme with the tableau

$$\begin{array}{cccccc}0& & & & & \\ \frac{1}{2}& \frac{1}{2}& & & & \\ \frac{1}{2}& 0& \frac{1}{2}& & & \\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & \\ \frac{1}{2}& 0& 0& \mathrm{}& \frac{1}{2}& \\ & & & & & \\ & 0& 0& \mathrm{}& 0& 1.\end{array}$$

Note, that $`k`$ iterations correspond to $`s=k+1`$ stages of the corresponding Runge-Kutta scheme. Application of this scheme to the model equation $`\dot{u}=\lambda u`$ yields the stability function

$$P_s(z)=1+z\underset{k=0}{\overset{s1}{}}\left(\frac{z}{2}\right)^k=1+z+\frac{z^2}{2}+\frac{z^3}{4}\underset{k=0}{\overset{s3}{}}\left(\frac{z}{2}\right)^k.$$

The fact that the first three terms agree with the Taylor expansion of $`\mathrm{exp}(z)`$, the stability function of the exact time evolution, shows that this scheme is second order accurate Hairer et al. (2000). The stability regions are obtained as the level set $`\{P_s(z)P_s(\overline{z})=1\}`$. The regions for some values of $`s`$ are shown in Fig. 12.

We see that the imaginary axis intersects the stability regions for $`s=3,4,7,8,\mathrm{}`$ in a finite interval while for the other values they have only the origin in common with the imaginary axis. To explain this somewhat unexpected result we look at the behaviour of the modulus of the stability function on the imaginary axis in a neighbourhood of the origin. Writing

$$P_s(z)=1+z\underset{k=0}{\overset{s1}{}}\left(\frac{z}{2}\right)^k=1+z\frac{1\left(\frac{z}{2}\right)^s}{1\frac{z}{2}}$$

and inserting $`z=\mathrm{i}y`$ into $`Q(z)=P_s(z)P_s(\overline{z})`$ we find

$$Q(\mathrm{i}y)=1+2\frac{(\mathrm{i}y/2)^{s+1}(1+(1)^{s+1})+(\mathrm{i}y/2)^{s+2}(1+(1)^s)+2(y/2)^{2s+2}}{1+(y/2)^2}.$$

This shows that for all $`s>1`$

$$Q(\mathrm{i}y)=\{\begin{array}{cc}1\left(\frac{\mathrm{i}}{2}\right)^sy^{s+2}+\mathrm{}\hfill & s\text{ even}\hfill \\ 1+\left(\frac{\mathrm{i}}{2}\right)^{s+1}y^{s+1}+\mathrm{}\hfill & s\text{ odd}.\hfill \end{array}$$

The behaviour of the stability regions is determined by the coefficient of the lowest order term. These lead to the sequence $`1`$, $`1/4`$, $`1/4`$, $`1/16`$, $`1/16`$, $`\mathrm{}`$ for $`s=1,2,\mathrm{}`$. The negative coefficients imply a decrease of $`Q(iy)`$ along the imaginary axis close to $`y=0`$ and hence they correspond to the cases where the stability region contains a part of the imaginary axis. The other cases, when the coefficient of the dominating term is positive, are those when the stability region just touches the imaginary axis at the origin. In these cases, systems with eigenvalues on or to the right of the imaginary axis such as in particular the advection equation and other hyperbolic equations cannot be stably evolved with ICN. Thus, we come to the same conclusion as in Teukolsky (2000), namely that ICN($`k`$) for $`k=2,3,6,7,\mathrm{}`$ is unstable for hyperbolic equations. Furthermore, since the method is always second order accurate regardless of how many iterations are performed, there is no use in performing more than two. In fact, in most practical cases it is probably more efficient to use a higher order Runge-Kutta method instead.

Another interesting aspect of the ICN scheme is its relationship with the Crank-Nicholson method. This is a second order accurate, implicit time stepping scheme whose stability function is given by

$$P(z)=\frac{1+z/2}{1z/2}$$

and whose stability region is the entire left plane. Taking the limit of $`P_s(z)`$ for $`s\mathrm{}`$ we find that in fact $`P_s(z)P(z)`$. However, absolute convergence holds only in the circle $`|z|<2`$. This has the following implication: iterating the ICN scheme until convergence to produce approximations to the exact Crank-Nicholson method works only within that circle. Thus, this scheme (which might be termed ‘the full ICN scheme’) also has a stability barrier and hence a limitation of the time step in contrast to the exact CN scheme which is unconditionally stable. Thus, we conclude that the ICN scheme is not an efficient alternative to higher order Runge-Kutta or implicit schemes like CN.


format: text

text:
# 1 Number, from 800 samples, of cases where still a consensus was reached, for 𝐿=17 and 31. For 𝐿=101, 301, and 1001 only 80 samples were run and the resulting numbers thus multiplied by 10.

Frustration from Simultaneous Updating in Sznajd Consensus Model

Dietrich Stauffer

Institute for Theoretical Physics, Cologne University

D-50923 Köln, Euroland

In the Sznajd model of 2000, a pair of neighbouring agents on a square lattice convinces its six neighbours of the pair opinion iff the two agents of the pair share the same opinion. Now we replace the usual random sequential updating rule by simultaneous updating and find that this change makes a complete consenus much more difficult. The larger the lattice is, the higher must be the initial majority for one of the two competing opinions to become the consensus.

Key Words: opinion dynamics, computer simulation, sociophysics

The application of cellular automata, Ising models and other tools of (computational or statistical) physics has a long tradition (Majorana 1942, Schelling 1971, Sakoda 1971, Callen and Shapero 1974, Galam et al 1981, Schweitzer 1997, Weidlich 2000). In particular one would like to know the conditions to reach a consensus out of an initially diverging set of opinions (Deffuant et al 2000, Kobayashi 2001, Hegselmann and Krause 2002). Most models assume that every agent is influenced by its neighbours and takes, for example, the opinion of the majority of them, or of a weighted average. The Sznajd model (Sznajd-Weron and Sznajd 2000; for a review see Stauffer 2002), on the other hand, assumes that every agent tries to influence its neighbours, without caring much about what they think first. Thus in the Sznajd model the information flows outward to the neighbourhood, as in infection or rumour spreading (Noymer 2001), while in most other models the information flows inward from the neighbourhood. Also, the Sznajd model takes into account the well-known psychological and political fact that “united we stand, divided we fall”; only groups of people having the same opinion, not divided groups, can influence their neighbours.

On the square lattice, where every site is occupied by an agent having one of two possible opinions +1 and $``$1, the most-studied Sznajd rule is: A pair of nearest neighbours convinces its six nearest neighbours of the pair opinion if and only if both members of the pair have the same opinion; otherwise the pair and its neighbours do not change opinion. Initially the opinions are distributed randomly, +1 with probability $`p`$ and $`1`$ with probability $`1p`$. This standard model then gave always a consensus, which for large lattices was that opinion which initially had a majority; if $`p=1/2`$ initially, then half of the cases ended with everybody having opinion +1, and the other half of the cases with the opposite opinion.

In these simulations random sequential updating was used, i.e. one of the $`L\times L`$ agents in the square lattice was selected randomly, and then one of its four neighbours to check if they share the same opinion. One time step was completed if on average each of the $`L\times L`$ agents was selected once as the first member of the pair. A. Iosselevitch (private communication) suggested to compare this rule with the simultaneous updating traditional for cellular automata, used also in the consensus model of Kobayashi (2001): Each pair is judged by its opinion at time $`t`$ to gives its six neighbours their possibly new opinions at time $`t+1`$. Now we can go through the lattice like a typewriter to find the first member of the pair; only for the second member of the pair a random selection is still needed. Going through the whole lattice once constitutes one time step.

In both random sequential as well as simultaneous updating, an agent can belong to the neighbourhoods of several convincing pairs. For random sequential updating, the agent then follows each pair in the order in which it receives orders, just like civil servants followed their various governments in Germany during the 20th century. For simultaneous updating, on the other hand, it does not know what to do if one pair has opinion +1 and another also neigbouring pair has the opposite opinion. It then feels frustrated and does nothing. i.e. it stays with its old opinion. (Similar frustration effects are known from some models of magnetism.) This frustration then hinders the development of a consensus.

With up to 800 samples, and $`L13`$ we never found a consensus at $`p=1/2`$. One has to use small lattices, or $`p`$ different from 1/2, to find all agents at the end having the same opinion. Fig.1 shows how the number (among 800) of samples without a consensus even after 10,000 time steps varies with lattice size $`L`$ and initial concentration $`p`$; if there was a consensus it was in favour of the initial majority. (For $`7L<13`$ also at $`p=1/2`$ rare cases of a consensus were found.) Fig.2 shows how the $`p`$ needed to get a consensus in half the cases varied for varying lattice size; the problem is by definition symmetric about $`p=1/2`$, and only $`p<1/2`$ is thus plotted in our figures.

Of course, reality differs from a square lattice. If everybody interacts with everybody else equally, without regard of a geometrical distance (Kobayashi 2001), this “mean field” problem may be suited for analytical solution or description by differential equation. More realistic may be a network where a few agents have lots of neighbours, and most have only a few neighbours, without a sharp boundary between celebrities and common folks (Albert and Barabasi 2002). Here an Ising model gave already an unusual phase transition (Aleksiejuk et al 2002), and a different Sznajd model with many possible opinions agreed with Brazilian election results.

Our simultaneous updating corresponds to formal committee meetings at times fixed for all participants, while random sequential updating corresponds to informal meetings of subgroups at various times. Our simulations then indicate that informal meetings have a higher chance to lead to a consensus.

I thank A. Iosselevitch for suggesting this work and the supercomputer center in Jülich for time on their Cray-T3E.

R. Albert, A.L. Barabasi, Statistical mechanics of complex networks. Reviews of Modern Physics 74: 47-97 (2002).

A. Aleksiejuk, J.A. Hołyst and D. Stauffer. Ferromagnetic Phase Transition in Barabási-Albert Networks. Physica A 310: 260-266 (2002)

E. Callen and D. Shapero, A Theory of Social Imitation, Physics Today, July 1974, 23-28.

G. Deffuant, D. Neau, F. Amblard and G. Weisbuch, Mixing beliefs among interacting agents, Adv. Complex Syst. 3: 87-98 (2000).

S. Galam, Y. Gefen and Y. Shapir, Sociophysics: A mean behavior model for the process of strike, J. Mathematical Sociology 9: 1-13 (1982).

R. Hegselmann and M. Krause, Opinion dynamics and bounded confidence, for Journal of Artificial Societies and Social Simulation 5 (2002).

J. Kobayashi, Unanimous Opinions in Social Influence Networks. J. Mathematical Sociology 24: 285-297 (2001).

E. Majorana, Il valore delle leggi statistiche nella fisica e nelle scienze sociali, Scientia 36: 58-66 (1942).

A. Noymer, The transmission and persistence of ‘urban legends’: Sociological application of age-structured epidemic models. J. Mathematical Sociology 25: 299-323 (2001).

F. Schweitzer, (ed.) Self-Organization of Complex Structures: From Individual to Collective Dynamics, Gordon and Breach, Amsterdam 1997.

J.M. Sakoda, The checkerboard model of social interaction, J. Mathematical Sociology 1: 119-132 (1971).

T.C. Schelling, Dynamic Models of Segregation, J. Mathematical Sociology 1: 143-186 (1971).

D. Stauffer, Monte Carlo Simulations of the Sznajd model, Journal of Artificial Societies and Social Simulation 5, No.1, paper 4 (2002) (jasss.soc.surrey.ac.uk).

K. Sznajd-Weron and J. Sznajd, Opinion Evolution in Closed Community. Int. J. Mod. Phys. C 11: 1157-1166 (2000).

W. Weidlich, Sociodynamics; A Systematic Approach to Mathematical Modelling in the Social Sciences. Harwood Academic Publishers, 2000.


text:
# High-performance planar light-emitting diode

## Abstract

Planar light-emitting diodes (LEDs) fabricated within a single high-mobility quantum well are demonstrated. Our approach leads to a dramatic reduction of radiative lifetime and junction area with respect to conventional vertical LEDs, promising very high-frequency device operation. Devices were fabricated by UV lithography and wet chemical etching starting from p-type modulation-doped Al0.5Ga0.5As/GaAs heterostructures grown by molecular beam epitaxy. Electrical and optical measurements from room temperature down to 1.8 K show high spectral purity and high external efficiency. Time-resolved measurements yielded extremely short recombination times of the order of 50 ps, demonstrating the relevance of the present scheme for high-frequency device applications in the GHz range.

Widespread diffusion of optical communication systems demands low-cost and high-performance light sources. High-radiance, wide-modulation-bandwidth light-emitting diodes (LEDs) would be promising candidates to replace laser sources, especially for short-haul in-fiber optic links. These sources offer many advantages thanks to the simplicity of the device structure, ease of fabrication, high reliability and simplified biasing arrangement. Conventional vertical LEDs, however, suffer from intrinsic bandwidth and efficiency limitations due to minority-carrier lifetime, internal capacitance, and non-radiative recombination processes.

Several fabrication schemes were proposed to improve performance, mainly based on heavily-doped materials lyon ; lyon2 ; chen . Introducing impurities strongly reduces the total recombination time leading to increased modulation bandwidth, but it also affects luminescence efficiency by introducing non-radiative recombination channels. The device performance is thus a compromise between modulation bandwidth and radiative efficiency. High-speed operation (in the GHz range) was reported with typical efficiencies of few $`\mu `$W/mA lyon ; lyon2 ; chen .

Here we demonstrate a different approach based on the reduction of junction-dimensionality and on a planar geometry.

Our approach takes advantage of the extremely short recombination times that can be obtained in a two-dimensional modulation-doped quantum well (QW)sanvitto . The doping scheme ensures an increased radiative efficiency by drastically reducing the non-radiative recombination channels. Additionally, our planar configuration leads to a drastic reduction of the parasitic capacitance. These aspects make planar-junction devices very promising from the point of view of modulation bandwidth and spectral purity. Furthermore, when an appropriate cavity is included in the design, our scheme may allow laser-emission by cold-electron injection into the active region. The planar configuration imposes no restrictions on cavity geometry, allowing the fabrication of both horizontal- and vertical-cavity devices.

The few existing realizations of pn junctions with lateral geometries saiki ; vaccaro ; north  are based on the amphoteric nature of Si in GaAs. Under optimized growth conditions, Si acts as an acceptor on the (311)A-oriented plane and as a donor on the higher-index (n11)A-oriented planes. Therefore, it is possible to obtain p- and n-type adjacent doped regions on a properly processed Si-doped GaAs substratemiller . This approach leads to a significant reduction of the junction area and to good performance in terms of electro-luminescence, but the presence of doping impurities in the conduction layers negatively impacts transport properties and the optical performance of the devices.

Our fabrication scheme is shown in the upper inset of Fig. 1(a). The starting system was a p-type modulation-doped Al0.5Ga0.5As/GaAs heterostructure with a two-dimensional hole gas (2DHG) confined within the GaAs quantum-well. This 2DHG constitutes the p-type portion of the final device. The heterostructure was processed into mesas and p-type Au/Zn/Au (5/50/150 nm) Ohmic contacts were evaporated and annealed (60 s at 460 C) in nitrogen overpressure.

The fabrication of the n-type region involved two processing steps: an ortophosphoric-acid-based etching solution ($`\mathrm{H}_3\mathrm{PO}_4:\mathrm{H}_2\mathrm{O}_2:\mathrm{H}_2\mathrm{O}=3:1:50`$ for 95 seconds) was used to completely remove the Be-doped Al0.5Ga0.5As layer from a portion of the mesa to obtain an intrinsic region within the QW. Following the etching procedure, a self-aligned, n-type contact consisting of AuGe(eutectic)/Ni/Au (120/20/140 nm) was evaporated. A 95 s annealing at 450 C in nitrogen overpressure completed the procedure. The n-type contact introduces donors into the host semiconductor, creating an electron gas within the GaAs layer below the metal pad, adjacent to the 2DHG.

Several pn junctions were fabricated following the above procedure, from a molecular-beam-epitaxy-grown heterostructure whose layer sequence is shown in the lower inset of Fig. 1(a). Nominal Be concentration was 1018 cm-3. Free-hole concentration in the 20 nm-wide QW was numerically calculated by a Poisson-Schrödinger solver which yielded a density of $`8.8\times 10^{10}`$ holes/cm2.

Devices were tested by current-voltage (I-V) measurements from room temperature down to 1.8 K. All fabricated devices presented very reproducible rectifying characteristics, with threshold around 1.5 V and negligible reverse-bias current ($`<20`$ nA at room temperature) in all the explored voltage range (see Fig. 1(a)). The diode-like behavior and threshold values consistent with the GaAs energy-gap represent a first evidence of the actual formation of the planar pn junction.

Spatially, spectrally and temporally-resolved electro-luminescence was measured in a wide range of temperatures in order to characterize the optical properties of the devices. Light emission was obtained from all the devices measured for forward biases above the turn-on value. Data were taken from room temperature down to 1.8 K. The spatial distribution of the emitted light was measured using an experimental set-up based on a low-vibration cold-finger cryostat. The light-collection system consisted of a 100$`\times `$ objective coupled to a multi-mode optical fiber. The system was mounted on a computer-controlled motorized x-y translation stage to allow spatially-resolved measurements. The signal was detected by a cooled photomultiplier after spectral filtering by a single-grating monochromator. Figure 1(b) shows an electro-luminescence profile at room temperature under 6 V forward bias superimposed onto the image of one device. The darker region between the contacts is the mesa, from which emission is observed. The spatial distribution of the emission was found to be markedly temperature and drive current sensitive. At the lowest temperatures, the emission was found to originate from around the p-type contact at moderately high biases, demonstrating that the recombining electrons are able to travel a considerable distance before recombination.

Figure 2 shows electro-luminescence spectra taken at 15 K under a 2 V forward bias. The emission peak is centered at 1522 meV with a full width at half maximum (FWHM) of 7.4 meV. Data show remarkably high spectral purity, much higher than available with conventional vertical or even with planar junctions obtained on patterned GaAs (311)A-oriented substrates.

The observation of intense electro-luminescence and its peculiar spectral features demonstrates the functionality of the pn planar-junction scheme.

Emitted electroluminescence power was measured as a function of injected current and temperature (from 2.5 K to room temperature). Assuming isotropic emission we were able to estimate device efficiency as the ratio of the integrated collected power over the injected current. Maximum efficiency was obtained at low temperatures and low injected currents (I $`0.187`$ mA) and was found to reach $`13\mu \mathrm{W}/\mathrm{mA}`$. By increasing the current, efficiency drops to about 1 $`\mu \mathrm{W}/\mathrm{mA}`$ (I $`2`$ mA) and then remains nearly constant up to 8 mA. At room temperature the current dependence of the emission efficiency is less pronounced and varies from 4.1 $`\mu \mathrm{W}/\mathrm{mA}`$ to 2.5 $`\mu \mathrm{W}/\mathrm{mA}`$ in the range of currents investigated.

The temporal evolution of the electroluminescence after short electrical excitation pulses was studied using a similar set-up consisting of a microscope objective, single grating spectrometer and a Si avalanche photodiode detector. Time-resolved electroluminescence was measured at 5 K after excitation with 150 ps voltage pulses. An exponential fit to the curves yielded typical decay times of around 300 ps, close to the resolution limit of the avalanche photodiode detector. Although our set-up did not allow the temporal characteristic to be fully resolved, the measurements demonstrate that the emission can be modulated on sub-ns time scales. The measurements also showed that emission from spatial positions away from the n-contact is delayed relative to emission next to the contact. We suppose this delay, which is of the order of a few nanoseconds, is due to the propagation delay of electrons travelling along the mesa.

High-temporal-resolution photoluminescence measurements were also made on the 2DHG using a Streak camera set-up. For these measurements the luminescence was excited by 1 ps optical pulses at a repetition rate of 80 MHz and detected using a spectrometer and Streak camera in synchroscan mode. The density of the hole gas was varied by applying a bias between an Ohmic contact and a semi-transparent Schottky contact on the mesa surface. Figure 3 shows how the photoluminescence lifetime varies with the applied Schottky gate bias. The lifetime of the photo-excited carriers was found to be extremely short, ranging from 52 ps to 242 ps. Lifetime increases sharply as the QW is depleted. In the presence of a dilute density of excess holes, photo-induced formation of positively charged excitons (X+) from free excitons (X) becomes importantshields . X+ excitons decay much more rapidly because of the presence of a massive particle in the final state which easily allows momentum conservation. By lowering the hole density, the formation of X+ excitons is suppressed resulting in an increase in the population of X. In this regime, which we observed for positive gate voltages, decay of neutral excitons dominates radiative recombination and leads to the increased measured lifetime. These measurements demonstrate the high quality of our QW structures and point out the suitability of this approach for the implementation of high-frequency devices in the GHz range.

In conclusion, we fabricated LEDs using a novel planar-design scheme. Lateral junctions between n- and p-type high-quality two-dimensional systems were fabricated starting from modulation-doped AlGaAs/GaAs heterostructures containing a 2DHG within a GaAs QW. Devices were characterized in a wide range of temperatures by current-voltage measurements, spatially and spectrally resolved electro-luminescence and emission-power measurements. Very good performance was obtained in terms of spectral purity and high efficiency. The small junction area and the extremely short minority-carrier lifetimes make these devices ideally suitable for high-frequency operation.

This work was supported in part by the European Commission through the FET Project SAWPHOTON. M.C. acknowledges support by the C.N.R.


text:
# 1 Introduction

## 1 Introduction

CP violation may have played a crucial role in the formation of matter in the universe. It must have, if the inflationary cosmology is right, thus ruling out the (very unpalatable) possibility that the baryon asymmetry was an initial condition of the Big Bang. If inflation took place, as most of us believe based on the growing observational evidence, it made the universe devoid of matter and set the stage for reheating. The baryon asymmetry then must have been produced at some later point through CP-violating processes.

## 2 CP violation and electroweak baryogenesis

That CP, as well as C and B must be broken for baryogenesis to work, was first noted by Sakharov . In 1967 the only reason for considering the baryon number violation was theorists’ ambitions. Neither experiment, nor favored theoretical models supported this hypothesis. In addition to breaking the symmetries, one needs an out-of-equilibrium state of the universe to generate the asymmetry.

The advent of the Standard Model and ’t Hooft’s discovery of a baryon number violating instanton provided a requisite source of $`B`$ violation. The Standard Model also incorporates C, P, and CP breaking. Although the baryon number violation is highly suppressed in the Standard Model at zero temperature, Kuzmin, Rubakov, and Shaposhnikov  realized that at high temperatures the baryon number violating processes could go unsuppressed. In addition, the universe might be out of thermal equilibrium at the time of the electroweak phase transition. Hence, Kuzmin, Rubakov, and Shaposhnikov  put forth a very plausible and appealing possibility that the baryon asymmetry might arise from the Standard Model physics, at the time of the electroweak phase transition .

Unfortunately, in the minimal Standard Model, this scenario does not work. First, the phase transition is too weak, unless the Higgs mass is $`M__H<45`$ GeV, which is ruled out by experiment. Second, the CP violation from the CKM matrix makes a vanishing contribution to baryon asymmetry because it is suppressed by a product of Yukawa couplings .

Baryogenesis in the supersymmetric extensions of the Standard Model is less problematic because the additional scalar states can provide both the new sources of CP violation  and a way to make the phase transition more strongly first-order . The phase transition is stronger if one of the stops is very light . The new sources of CP violation may come, for example, from the chargino mass matrix:

$$\overline{\psi }_RM_\chi \psi _L=(\overline{\stackrel{~}{w}^^+},\overline{\stackrel{~}{h}_2^^+})_R\left(\begin{array}{cc}m_2& gH_2(x)\\ gH_1(x)& \mu \end{array}\right)\left(\begin{array}{c}\stackrel{~}{w}^^+\\ \stackrel{~}{h}_1^^+\end{array}\right)_L+\mathrm{h}.\mathrm{c}.$$
(1)
As long as $`m_2`$ and $`\mu `$ are complex, spatially varying phases in the bubble wall provide a source of (spontaneous) CP violation . The remaining window for electroweak baryogenesis in the MSSM is very narrow ; several parameters must be adjusted to maximize the resulting baryon asymmetry (in particular, one must assume that the wall is very thin, take $`\mathrm{tan}\beta <3`$, and choose the “optimal” bubble wall velocity $`v_w0.02`$), as shown in Fig. 1.

CP violation from the chargino sector (1) may enhance the $`B_s`$ mixing as compared to the Standard Model value , especially if the stop is light (Fig. 2). In practice, however, this effect is observable only if the CKM matrix elements are known to a very high precision. In particular, one would need to reduce theoretical uncertainties in $`V_{ub}`$ to 5-10% and in $`\mathrm{sin}2\beta `$ to a few percent.

## 3 CP violation and leptogenesis

If a lepton asymmetry of the universe formed after inflation but before the electroweak phase transition, the sphalerons, which violate $`B`$ and $`L`$ but preserve $`(BL)`$, would convert (roughly, a half of) the lepton asymmetry into a baryon asymmetry. This observation gave rise to an extremely appealing scenario of leptogenesis . The relevant CP violation may reside in the neutrino mass matrix, which, in general, has a number of complex phases .

## 4 Transient CP violation, baryogenesis, and dark matter

In addition to CP violation hard-wired into the lagrangian, some manifestations of CP non-conservation may occur for a short period of time in the early universe, during which time the baryon asymmetry might have formed. A well-know example is the Affleck-Dine scenario for baryogenesis , in which CP violating seeds may be effectively amplified by the motion of the scalar condensate. Similarly, electroweak baryogenesis at preheating  can take advantage of CP-violating motions of time-dependent condensates during preheating . CP violation of this kind is poorly constrained by experiment because it becomes small after thermalization.

### 4.1 Affleck-Dine baryogenesis and dark-matter-genesis

In models with low-energy supersymmetry inflation can lead to formation of the Affleck-Dine condensate  with a large VEV. A high-scale physics undoubtedly violates B and CP through higher-dimension operators. Hence, the motion of the scalar condensate after inflation is not B and CP symmetric. Thus, the universe acquires a baryon asymmetry. The Affleck-Dine scenario  is simple, appealing, and flexible in that the final baryon asymmetry can easily be made consistent with the data.

In addition, the Affleck-Dine baryogenesis can produce dark matter as well. In general, the Affleck-Dine condensate does not remain homogeneous and can break up  into SUSY Q-balls . This affords a number of interesting possibilities for generating baryons and dark matter simultaneously .

Since baryons and dark matter arise form the same physical process,

the ratio of $`\mathrm{\Omega }_{\mathrm{dark}}`$ to $`\mathrm{\Omega }_B`$ may have a natural explanation in some models .

### 4.2 Electroweak baryogenesis at preheating

Several viable scenarios for electroweak baryogenesis at preheating were presented in Ref. . For example, a modified spontaneous baryogenesis a la Cohen, Kaplan, and Nelson  becomes very efficient at preheating. Their original scenario used the variation of the Higgs field inside a wall of a bubble formed in a first-order phase transition. A similar effect can occur at preheating uniformly in space, on the horizon scales . One can obtain the desired baryon asymmetry in a Standard Model supplemented by an additional Higgs doublet and an inflaton sector . The difference with the scenario proposed by Cohen, Kaplan and Nelson  is that in our case CP violation occurs homogeneously in space, not only inside a bubble wall. In addition, the final prediction for the baryon asymmetry in the CKN scenario was very far from the equilibrium value because the sphaleron rate was slow on the time scales associated with the growth of bubbles. In the case of preheating, the Higgs parameters change slowly in time while the baryon number non-conservation is rapid. This allows a slow adiabatic adjustment of the baryon number to that which minimizes the free energy.

Several additional sources of CP violation might affect the physics of preheating and facilitate baryogenesis .

## 5 Strong CP violation and the axion cosmology

The QCD vacuum is a superposition $`|\theta =_n\mathrm{exp}\{in\theta \}|n`$ of topologically distinct vacuum states $`|n`$. As a result, the QCD Lagrangian can be written as

$$_{QCD}=_{\mathrm{pert}}+\overline{\theta }\frac{g^2}{32\pi ^2}F\stackrel{~}{F},$$
(2)
where

$$\overline{\theta }=\theta +\mathrm{arg}detM$$
(3)
Experimentally, the value of $`\overline{\theta }`$ must vanish to a high precision, $`\overline{\theta }10^{10}`$. However, if all the quarks have non-zero masses, there is no (simple) reason why the phase in the Yukawa matrix should cancel the QCD vacuum phase, unless $`\overline{\theta }`$ relaxes to zero dynamically, by a VEV of a scalar field. This elegant solution to the strong CP problem was proposed by Peccei and Quinn . The breaking of a global U(1) symmetry gives rise to a light scalar field, the axion .

Several models can accommodate an axion consistent with the existing experimental bounds . A light, weakly interacting axion makes a good candidate for dark matter. The present experimental limits are shown in (Fig. 3).

## 6 Conclusions

There is every reason to believe that matter-antimatter asymmetry is a consequence of CP non-conservation in particle physics. On the other hand, CP violation from the quark mixing is not sufficient for baryogenesis. This implies the existence of new, yet undiscovered, sources of CP violation in nature.

## 7 Acknowledgments

The author thanks J. Cline and H. Murayama for helpful discussions. This work was supported in part by the DOE grant DE-FG03-91ER40662.