text:
# Modeling of multicomponent radiatively driven stellar winds using a Newton-Raphson method

## 1. Introduction

Radiatively driven stellar winds are accelerated by the absorption of radiation mainly in the resonance lines of such elements like C, N, O or Fe (see Castor, Abbott, & Klein 1975, hereafter CAK; Abbott 1982). These wind components have much lower density than the stellar wind itself, which is composed mainly by hydrogen and helium. The momentum obtained by the low-density absorbing wind component is transfered to the high-density nonabsorbing wind component via collisions between electrically charged particles (Springmann & Pauldrach 1992; Babel 1995). Clearly, multicomponent radiatively driven winds of hot stars have a multicomponent nature. This multicomponent nature is important mainly for the low-density stellar winds, which shall be described by the multicomponent hydrodynamic equations.

Krtička & Kubát (2000; 2001a; 2001b) calculated models of such multicomponent stellar winds using Newton-Raphson method. In this paper we describe in detail method of solution of these multicomponent hydrodynamic equations.

## 2. Model equations

Each wind component is described by a density $`\rho _a`$, radial velocity $`v_{r}^{}{}_{a}{}^{}`$, temperature $`T_a`$ and electrical charge $`q_a`$. Typically, we take into account only three wind components (ie. absorbing ions, nonabsorbing ions and electrons).

Stationary spherically symmetric multicomponent radiatively driven stellar winds are described by the continuity equation for each component $`a`$ of the wind

$$\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2\rho _av_{r}^{}{}_{a}{}^{}\right)=0,$$
(1)
(note, that in the case of continuity equation of electrons a right-hand term occurs due to the change of the wind ionization in the flow, see Krtička & Kubát, 2001b), momentum and energy equations for each component $`a`$

$`v_{r}^{}{}_{a}{}^{}{\displaystyle \frac{\mathrm{d}v_{r}^{}{}_{a}{}^{}}{\mathrm{d}r}}=g_a^{\mathrm{rad}}g{\displaystyle \frac{1}{\rho _a}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}\left(a_a^2\rho _a\right)+{\displaystyle \frac{q_a}{m_a}}E+{\displaystyle \underset{ba}{}}g_{ab}^{\mathrm{fric}},`$ (2)
$`{\displaystyle \frac{3}{2}}v_{r}^{}{}_{a}{}^{}\rho _a{\displaystyle \frac{\mathrm{d}a_a^2}{\mathrm{d}r}}+{\displaystyle \frac{a_a^2\rho _a}{r^2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}\left(r^2v_{r}^{}{}_{a}{}^{}\right)=Q_a^{\mathrm{rad}}+{\displaystyle \underset{ba}{}}Q_{ab}^{\mathrm{ex}}+{\displaystyle \underset{ba}{}}Q_{ab}^{\mathrm{fric}},`$ (3)
where $`a_a`$ is isothermal sound speed, $`E`$ is charge separation electric field, $`g`$ is gravitational acceleration, radiative acceleration $`g_a^{\mathrm{rad}}`$ acts on free electrons and absorbing ions. Radiative acceleration acting on absorbing ions is taken in the CAK approximation with the finite disk correction factor (Pauldrach, Puls, & Kudritzki 1986) and with force multipliers after Abbott (1982), slightly modified for the multicomponent case (Krtička & Kubát 2000). Frictional acceleration $`g_{ab}^{\mathrm{fric}}`$, frictional heating $`Q_{ab}^{\mathrm{fric}}`$ and exchange of the heat $`Q_{ab}^{\mathrm{ex}}`$ between wind components $`a`$ and $`b`$ are given by Burgers (1969; see also Krtička & Kubát 2001b).

The only transitions which contribute to the radiative heating/cooling term $`Q_a^{\mathrm{rad}}`$ in the static atmosphere are bound-free and free-free transitions (Kubát, Puls, & Pauldrach 1999; see also Kubát, this volume). Because these transitions deposit energy directly to electrons, this radiative heating/cooling term is considered in electron energy equation. Mean intensity $`J_\nu `$ at the base of the wind is taken as an emergent radiation from a spherically symmetric static hydrogen model atmosphere (Kubát 2001). In the case of moving media additional heating/cooling term occurs, so-called Gayley-Owocki heating, which is caused by the dependence of a radiative force on a velocity via Doppler effect (Gayley & Owocki 1994; Krtička & Kubát 2001b).

The electrical charge $`q_a`$ is obtained using the approximate ionization equilibrium described by Mihalas (1978, eq. \[5-46\]).

The system of hydrodynamic equations (1, 2, 3) is closed using an equation for polarisation electric field in the form of

$$\frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2E\right)=4\pi \rho =4\pi \underset{a}{}\frac{q_a\rho _a}{m_a}.$$
(4)
Finally, the system of equations solved is supplemented by appropriate boundary conditions. We start to calculate our models at the critical point of nonabsorbing component. We also use zero current condition, the condition of quasi-neutrality and assume that the flow at the inner boundary is in radiative equilibrium and that the boundary temperature of all components is the same. Boundary density of absorbing ions is determined using abundance ratios.

## 3. Method of solution

Velocities, densities and other unknown variables are calculated at grid points. We selected logarithmically spaced grid in radii with accumulation of grid points near the stellar photosphere,

$$r_i=qr_{i1},q=\mathrm{exp}\left(\frac{\mathrm{ln}\frac{r_{\mathrm{N}\mathrm{R}}}{r_1}}{\mathrm{N}\mathrm{R}1}\right),$$
(5)
where $`i`$ is a number of a given grid point and $`\mathrm{N}\mathrm{R}`$ is number of grid points. We use typically 50–300 grid points.

For the approximation of the derivative of expression $`X`$ we selected

$$\frac{\mathrm{d}X}{\mathrm{d}r}|_{r=r_i}\{\begin{array}{cc}y_i\frac{X_{i+1}X_i}{\mathrm{\Delta }r_{i+1}}+(1y_i)\frac{X_iX_{i1}}{\mathrm{\Delta }r_i},& i<\mathrm{N}\mathrm{R},\\ \frac{X_iX_{i1}}{\mathrm{\Delta }r_i},& i=\mathrm{N}\mathrm{R},\end{array}$$
(6)
where

$`\mathrm{\Delta }r_i`$ $`=r_ir_{i1},`$ $`\overline{r}_i`$ $`={\displaystyle \frac{1}{2}}\left(r_i+r_{i1}\right),`$ $`y_i`$ $`={\displaystyle \frac{r_i\overline{r}_i}{\overline{r}_{i+1}\overline{r}_i}}.`$ (7)
An appripriate approximation of the derivative is important, numerical tests showed that approximation (6) gives the best convergence of the models.

An approximation of hydrodynamical equations (1, 2, 3) together with the equation for a polarisation electric field (4), equation of ionization equilibrium, and boundary conditions can be formally written as

$$𝖯𝝍=0,$$
(8)
where the vector describing the solution has the form of

$$𝝍=(𝝍_1,𝝍_2,\mathrm{},𝝍_{\mathrm{N}\mathrm{R}})^\mathrm{T},$$
(9a)
$$𝝍_i=(\underset{a}{}\rho _{a,i},\underset{a}{}v_{ra,i},\underset{a}{}T_{a,i},\underset{a}{}z_{a,i},E_i,\mathrm{\Delta }v_{r,i})$$
(9b)
and $`𝖯`$ consists from all equations mentioned above.

For the solution of the nonlinear system of model equations (8) we selected Newton-Raphson method. The solution of Eqs.(8) can be obtained using iterative procedure in the form of

$$𝖩^n\delta 𝝍^{n+1}=𝖯^n𝝍^n,$$
(10)
where $`𝝍^n`$ denotes solution in the $`n`$-th iteration, $`\delta 𝝍^{n+1}`$ is a correction of the solution and the Jacobian is

$$J_{kl}^n=\frac{P_k}{\psi _l}.$$
(11)
The Newton-Raphson method is similar to the well-known method of complete linearization (Auer & Mihalas 1969). The analytical form of Jacobi matrix $`J_{kl}^n`$ can be easily obtained from Eqs.(8) and for the case of three-component radiatively driven stellar wind was calculated by Krtička (2001). However, the value of corrections $`\delta 𝝍^{n+1}`$ in each iterative step is limited by $`0.7𝝍^n`$ during the calculation of a model below the critical point and by $`0.1𝝍^n`$ during the calculation of a model above the critical point. The velocity differences between wind components are limited similarly.

Another problem for the calculation of wind models is the inclusion of critical point condition. Critical point is a point where the wind velocity is equal to the sound speed. Generally, there are two possibilities of inclusion of the critical point condition. Nobili & Turola (1988) showed that critical point condition can be directly included into the set of linearised equations. However, due to some numerical problems, we selected the so called ”shooting method”. We change the density at the base of the wind in order to fulfil the critical point condition. When the base density and model upstream the critical point are obtained, we calculate model downstream the critical point. The initial model for the Newton-Raphson iterations is selected in such a way that it converges to the correct branch of CAK solution appropriate for the flow downstream the critical point (see CAK).

For the solution of the system of linearised equations (10) we use numerical package LAPACK. We also applied the Gaussian elimination, however we found that the package LAPACK is faster and even more accurate.

Newton-Raphson method represents quite powerful tool for the solution of stationary hydrodynamic equations. The iterations converge very fast (typical relative change is of order $`10^{10}`$ after 20 iterations), the number of equations solved can be easily extended and it is possible to calculate multidimensional models. On the other hand the model convergence depends on initial estimate and we should also check the dynamical stability of obtained solution.

### Acknowledgments.

The author is greatly indebted to Dr. J. Kubát for discussion and cooperation on this topic. This work was supported by grants GA ČR 205/01/0656, 205/02/0445, and by projects K2043105 and AV0Z1003909.

## References

Abbott, D. C., 1982, ApJ, 259, 282

Auer, L. H., Mihalas, D., 1969, ApJ, 158, 641

Babel, J., 1995, A&A, 301, 823

Burgers, J. M.,1969, Flow equations for composite gases, Academic Press, New York

Castor, J. I., Abbott, D. C., & Klein, R. I., 1975, ApJ, 195, 157 (CAK)

Gayley, K. G., & Owocki, S. P., 1994, ApJ, 434, 684

Krtička, J., 2001, Ph.D. thesis, Masaryk University Brno

Krtička, J., Kubát, J., 2000, A&A, 359, 983

Krtička, J., Kubát, J., 2001a, A&A, 369, 222

Krtička, J., Kubát, J., 2001b, A&A, 377, 175

Kubát, J., 2001, A&A, 366, 210

Kubát, J., Puls, J., & Pauldrach, A. W. A., 1999, A&A, 341, 587

Mihalas, D., 1978, Stellar Atmospheres. 2nd ed., W.H.Freeman & Comp., San Francisco

Nobili, L., & Turola, R., 1988, ApJ, 333, 248

Pauldrach, A., Puls, J., & Kudritzki, R. P., 1986, A&A, 164, 86

Springmann, U. W. E., & Pauldrach, A. W. A., 1992, A&A, 262, 515


format: text

text:
# 1 The ratio R≡|S₁ᶢ(M_H₁²)/Sᶢ_SM(M_H_SM²)| assuming M_H_SM=M_H₁ in the CPX scenario [] with M_SUSY=0.5 TeV, Φ_CP=90^∘, and (a)tanβ=4 and (b)tanβ=5.

## Acknowledgements

We should like to thank Beate Heinemann, Valery Khoze and Misha Ryskin for many helpful discussions. BEC and JSL thank PPARC for their support in funding part of this research.


text:
# Genus Topology of the Cosmic Microwave Background from WMAP

## 1 Introduction

The greatly anticipated results from the WMAP project (Bennett et al. 2003a) have redefined the state of the art in Cosmic Microwave Background (CMB) science. The WMAP project has provided full-sky coverage of the CMB at unprecedented angular resolution. The improvement over the only previous full-sky dataset (COBE, Smoot et al. 1992) is of order 500 in terms of resolution elements, and the sky coverage improvement over balloon and ground-based experiments (e.g., de Bernardis et al. 2000 \[Boomerang\], Stompor et al. 2001 \[Maxima\], Mason et al. 2002 \[CBI\], Kovac et al. 2002 \[DASI\]) is several hundred. The experiment has resulted in dramatically improved constraints on cosmological parameters, such as the matter density, $`\mathrm{\Omega }_m`$, Hubble Constant, $`H_0`$, and cosmological constant, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ (Spergel et al. 2003).

The tool of choice for assessing the cosmological parameters, is the power spectrum, where one calculates the products of the spherical harmonic coefficients $`a_\mathrm{}m`$ and their complex conjugates (for each $`\mathrm{}`$ and $`m`$ in the $`Y_\mathrm{}m`$ expansion of the CMB sky, summed over $`m`$ values to form $`C_{\mathrm{}}`$). This $`a_\mathrm{}ma_\mathrm{}m^{}`$ product, however, explicitly removes phase information. And while the power-spectrum is a very powerful tool for assessing cosmological parameters such as $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$, it does (explicitly) remove any phase information contained in the $`a_\mathrm{}m`$ coefficients themselves. These phases, however, contain critical information for characterizing the primordial density fluctuations. Namely, standard inflation (e.g., Guth 1981, Albrecht & Steinhardt 1982, Linde 1982, Linde 1983) predicts that the temperature fluctuations in the CMB, at the resolution measured by WMAP, will be characterized by spherical harmonic coefficients with Gaussian distributed amplitudes and random phases. The WMAP data provide our best opportunity to date to test that hypothesis.

The genus topology method developed by Gott, Melott & Dickinson (1986) directly tests for the Gaussian random-phase nature of a density (or temperature) distribution in 3 dimensions (Adler 1981; Gott, Melott & Dickinson 1986; Hamilton, Gott & Weinberg 1986; Gott, Weinberg & Melott 1987), or in 2 dimensions (Adler 1981; Melott et al. 1989). Coles (1988) independently developed an equivalent statistic in 2 dimensions. The 2 dimensional case has been studied for a variety of cosmological datasets: on redshift slices (Park et al. 1992; Colley 1997), on sky maps (Gott et al. 1992; Park, Gott, & Choi 2001; Hoyle, Vogeley & Gott 2002); and on the CMB, in particular (Smoot et al. 1992; Kogut 1993; Kogut et al. 1996; Colley, Gott & Park 1996; Park, C-G. et al. 2001).

The WMAP team has recently measured the genus of the WMAP sky (Komatsu et al. 2003), and demonstrated that the WMAP results are consistent with the Gaussian random-phase hypothesis. To do this, they carried out a large number of simulations of the CMB, in which the spherical harmonic coefficients were drawn from a Gaussian random-phase distribution. They then used their known beam profiles, to simulate the results in each frequency, and applied the Kp0 (Bennett et al. 2003b) mask, just as one would with the real dataset. For each of these simulations, they computed the the genus (as defined by Melott et al. ) and compared it to the to that of the real dataset. They found that the real dataset does not depart significantly from the genus of the simulated Gaussian random-phase datasets, very much in agreement with the results predicted for WMAP by Park et al. (1998). That group further explored Gaussianity by testing other Minkowski functionals (Minkowski 1903) and the bispectrum (related to the three-point correlation function), and found constraints on Gaussianity similar to those provided by the genus study.

We seek, first, to confirm this result using our own methods, but also to compare the genus measured from the data directly to the theoretical genus curve for structures on a sphere. Rather than simulating Gaussian random-phase realizations of the CMB, we compare directly to the theoretical prediction for the genus in two dimensions. Furthermore, we provide some details about the nuances of carrying out the genus calculation on the HEALPix map projection and the stereographic map projection.

## 2 WMAP Observations

The Wilkinson Microwave Anisotropy Probe (WMAP)111The WMAP homepage is at http://map.gsfc.nasa.gov/ project team has released its one-year dataset in several different skymaps measured at varying frequency and angular resolution. Namely, the project has measured the anisotropy at frequencies of 23, 33, 41, 61, and 94 GHz at angular resolutions of 0.82, 0.62, 0.49, 0.33, and 0.21 degrees (FWHM beamwidths), respectively.

The WMAP team has released the data in a unique map projection called the HEALPix projection (Górski et al. 2000). We have chosen to re-present the WMAP results in Fig. 1 for a few reasons. First, we will be studying the genus on stereographic projections (sections 6 and 7), so it is useful to see the full sky map in that projection. Second, we have chosen a much different color scheme that is more relevant to genus studies than the original color scheme used by the WMAP team. The WMAP color scheme is pretty, informative and impressive, but is not symmetric with respect to hot and cold spots. In our color scheme, the mean temperature contour is white. Higher temperatures than the mean are represented as linearly redder with temperature. Lower temperatures than the mean are represented as linearly bluer with temperature. Therefore, the amount of blue ink per pixel on the page is linearly proportional to the coldness of a cold spot, and the amount of red ink per pixel is linearly proportional to the the hotness of a hot spot. The scale runs from $`200\mu \text{K}`$ (solid blue) to $`+200\mu \text{K}`$ (solid red), the same range as in the temperature maps produced by the WMAP team (Bennett et al. 2003a). Since the two-dimensional genus of a Gaussian random-phase field is all about the symmetry between hot spots and cold spots, we find this color map highly instructive for our purposes.

## 3 Genus Topology on a Sphere

The properties of the genus are well-known in three dimensions (3D) (Gott et al. 1986, Hamilton et al. 1986, Gott et al. 1987), but require some explanation in the two-dimensional (2D) case (Melott et al. 1989), particularly in the case of the sphere (Gott et al. 1990).

For 2D topology on a plane the 2D genus of a microwave background map is defined as (Melott et al. 1989):

$$g_{2D}=\text{Number of hot spots}\text{Number of cold spots}$$
(1)
For a Gaussian random field,

$$g_{2D}\nu \mathrm{exp}(\nu ^2/2),$$
(2)
where $`\nu `$ is a parameter that measures the area fraction in the hot spots:

$$f=(2\pi )^{1/2}_\nu ^{\mathrm{}}\mathrm{exp}(x^2/2)𝑑x.$$
(3)
So for $`\nu >0`$ ($`f<0.5`$) there are more hot spots than cold spots, while for $`\nu <0`$ ($`f>0.5`$) there are more cold spots than hot spots. The genus is also equal to the integral of the curvature around the temperature contour divided by $`2\pi `$. If we were to drive a truck around an isolated hot spot we would have to turn a total angle of $`2\pi `$ as we complete an entire circuit around the hot spot. Driving a truck around an isolated cold spot, we would turn a total angle of $`2\pi `$ with the opposite sign, with a negative turn angle defined as one that is a turn to the left when the hot region is on your right. Thus, in a temperature field that is divided into pixels we may define a pixel as hot if it is above the contour threshold and cold if it is below the contour threshold. Then the contour line is actually composed of a series of line segments with turns occurring at vertices in the pixel map. There is a certain turn occurring at each vertex (which, divided by $`2\pi `$, is equal to the contribution to the genus). Around each vertex there are four pixels, the contribution to the total turn (and genus) at that vertex depends on how many of the four surrounding pixels are hot and cold:

| 0 hot, 4 cold | turn = 0, genus = 0 | the contour line does not intersect the vertex |
| --- | --- | --- |
| 1 hot, 3 cold | turn = $`90^{}`$, genus = $`+1/4`$ | the contour line makes a right angle turn |
| 2 hot, 2 cold | turn = 0, genus = 0 | no turn, or turns that add to zero on average |
| 3 hot, 1 cold | turn = $`90^{}`$, genus = $`1/4`$ | right angle turn around a cold spot |
| 4 hot, 0 cold | turn = 0, genus = 0 | contour line does not intersect the vertex. |

This is carried out by an identical algorithm to that in the CONTOUR2D program as described in Melott et al. (1989).

Now we wish to define rigorously the 2D genus on a spherical surface. Here we will follow our previous derivation (Gott et al. 1990). The 2D genus is defined to be equal to minus the 3D genus of the hot regions confined within a thin spherical shell, measured in 3D. Use lead paint to paint the hot regions on the surface of a balloon—burst the balloon and you will have solid curved lead shapes that will have a certain 3D genus—take the minus of this number and that will be the 2D genus as we will define it.

Recall that we have defined the 3D genus as the number of (doughnut) holes minus the number of isolated regions. Our 3D genus is equal to the integral of the Gaussian curvature over the contour surface divided by $`4\pi `$. Thus, a sphere in 3D has a genus of $`1`$, because it is one isolated region. The Gaussian curvature is $`+(1/r^2)`$ while the area of the sphere is $`4\pi r^2`$, so the integral of the Gaussian curvature over the sphere is $`4\pi `$, independent of the size of the sphere, and the genus is $`1`$. A cube has a genus of $`1`$ also, because it has the same topology as a sphere. In this case, the curvature of the surface entirely consists of 8 delta functions at the eight vertices where three squares meet at a point introducing a contribution to the Gaussian curvature integral of $`\pi /2`$ at each vertex, equal to the conical deficit angle at the vertex. (If we sandpapered off the corners and edges of the cube, we would produce at each vertex an octant of a sphere containing a contribution to the Gaussian curvature integral of 1/8th that of a sphere or $`\pi /2`$. The edges would be sanded off to quarter cylinders which have no Gaussian curvature and the faces would have no Gaussian curvature, so all the curvature is located at the vertices). A doughnut has a genus of 0, because it has one hole and consists of one isolated region. A sphere with two handles has a genus of $`+1`$ on our definition because it consists of one isolated region and has two holes.

So using these definitions we can define the genus of a microwave background map on the celestial sphere. Suppose for example we have one hot spot in the north polar region, and the rest of the sphere is cold. Then the genus would be $`+1`$, because the hot spot cap is one isolated region. Now imagine the hot region covers the northern hemisphere while the southern hemisphere is cold. The genus would still be $`+1`$, because a hemispherical bowl is one isolated region. Suppose the hot region covers all of the sphere except for a cold spot in the south polar region. The genus would still be $`+1`$, because this would look like a sugar bowl with out any handles, which is also one isolated region. The topology in each case is identical since one can be deformed into the other, so the genus should be the same in all three cases. If the hot region ran around the equatorial regions and there were cold spots at the north and south polar regions, the 2D genus would be 0 because in 3D the equatorial band (confined in a thin spherical shell) is a doughnut in 3D which has a 3D genus of 0. If the sphere had $`N`$ isolated hot spots (polka dots) on a cold background, the 2D genus would be $`N`$. If the sphere had $`N`$ isolated cold spots on a hot background the genus would be $`2N`$. (Remember that 1 isolated hot spot on a cold background \[a hot spot centered on the north pole; the rest cold\] is the same topologically as an isolated cold spot centered on the south pole \[the rest hot\], because the latitude of the contour line can be simply moved downward to transform one into the other.) Suppose we have a hot spot covering the northern hemisphere with the southern hemisphere cold. This is a spherical cap of area $`2\pi `$ steradians. Its 2D genus is $`+1`$ because in 3D it is one isolated region. Considered as a spherical cap (part of a thin spherical shell) this isolated region must have an integral of the Gaussian curvature on its surface equal to $`4\pi `$ by the Gauss-Bonet theorem. There are several contributions to this: there is a contribution of $`2\pi `$ from the outside (spherical) cap surface, and a contribution of $`2\pi `$ from the inside (spherical) cap surface (a spherical shell has an inside and an outside boundary both spherical), and a net contribution of 0 from the delta function curvatures on the inside and outside edges at the equator marking the boundaries of the inner and outer spherical caps—the contribution of the outer edge is positive, while the contribution of the inner edge is saddle-shaped and therefore negative. The equator is a geodesic, so we could drive a truck around the equator without turning; hence, the integral of the turning around the contour boundary is zero. Since we are circling an isolated hot spot, why is the total turning not $`2\pi `$ as it is on the plane? Because we are circling a curved region, and if we parallel transport (without turning) a vector around the boundary it will suffer a deflection equal to the integral of the Gaussian curvature in the interior of the boundary, which in this case is $`2\pi `$, the value we expect for a closed curve in the plane. The vector must return to its initial position after circling, so the total of the turning integral and the parallel transport deflection due to circling a curved region must add to give $`2\pi `$.

Now consider a hot spot centered at the north pole and whose boundary extends to latitude $`30^{}\text{N}`$. This is a spherical cap with area $`\pi `$. It is an isolated region so the 3D genus is $`1`$, and the sum of the integral of the Gaussian curvature over the spherical cap in 3D must be $`4\pi `$. The integral on the outside of the spherical cap is $`\pi `$, and on the inside of the spherical cap is also $`\pi `$, and the boundary surface connecting the inner and outer spherical cap shells is locally flat (part of a cone), so the contribution from the inside and outside edges of the boundary must be $`2\pi `$, which must be equal to twice the turning integral for a truck driving on the spherical surface (since the total Gaussian integral, $`4\pi `$, is twice the $`2\pi `$ value for circling an isolated boundary \[turning plus deviation from parallel transport\]). If we were to parallel transport a vector around this boundary, we would find a deflection of $`\pi `$ because that is the area in the spherical cap; to return to its original orientation, a rotation of $`2\pi `$, the truck must turn an additional angle of $`\pi `$. Indeed, if you would like to drive counterclockwise around the sphere at latitude $`30^{}\text{N}`$, you will have to keep your steering wheel turned to the left, since this not a geodesic, and the total of that turning integral is $`\pi `$. Now, locally it is the turning integral that is being calculated from the formula for a Gaussian random field. So let us define the effective genus as

$$g_{2D,eff}=g_{2D}2f,$$
(4)
where $`f`$ is the fraction of the area of the sphere in the hot spots. In our example, where the one hot spot was the entire northern hemisphere, $`g_{2D}=1`$, $`f=1/2`$ and $`g_{2D,eff}=0`$. The value of $`g_{2D,eff}`$ is equal to the turning integral of a truck driving around all the individual components of the boundary contour (around all the hot spots) divided by $`2\pi `$; in this case it is 0 because the truck driving around the equator drives straight ahead without turning. For the case where the hot spot is centered at the north pole and extends to latitude $`30^{}\text{N},f=1/4`$, so the value of $`g_{2D,eff}=1/2`$, because the turning integral in this case is $`\pi `$. Now for a Gaussian random field on the sphere

$$g_{2D,eff}\nu \mathrm{exp}(\nu ^2/2),$$
(5)
because the Gaussian random field locally behaves on the sphere as on the plane to produce this contribution to the turning integral. Thus, in comparing the WMAP data to the random-phase formula we will use $`g_{2D,eff}`$ defined rigorously as described above.

## 4 Genus in the HEALPix Projection

The WMAP data are plotted using an unusual map projection of the sphere, called the HEALPix projection (Górski et al. 2000). The sphere is effectively projected onto a rhombic dodecahedron. (This semi-regular polyhedron is the cell for the face centered cubic crystal lattice, the polyhedron representing the boundary of the set of points closer to a particular atomic nucleus than to any other.) It has 12 faces, four diamonds that meet at the top like the four sides of a pyramid, four diamonds that circle the equator, and four that meet at the bottom like the sides of an inverted pyramid. If we were to set this rhombic dodecahedron up aligned with the cardinal directions we could label the 12 faces as follows: TN (top north), TE (top east), TS (top south), TW (top west), NE (northeast, on the equator), SE (southeast, on the equator), SW (southwest, on the equator), NW (northwest, on the equator), BN (bottom north), BE (bottom east), BS (bottom south), BW (bottom west). There are 6 vertices where 4 diamonds meet at a point (i.e., TN, TE, TS, TW all meet at the top, and similarly TN, NE, NW, BN meet at a vertex). There are 8 vertices where 3 diamonds meet at a point (e.g., TN, TE, NE meet at a vertex). Each diamond may be subdivided into diamond shaped pixels. This is done by factors of 4, subdividing each diamond into 4 diamonds, and repeating by similarly subdividing each diamond. (In a similar way we may produce a checkerboard of 64 square pixels by dividing a square first into 4 quadrant squares, dividing each of these into 4 sub-squares, and each of these sub-squares into 4 checkerboard squares). In the WMAP case these are all diamonds (rhombuses) instead of squares, but the topology is the same. The subdivisions are made so that all the pixels are equal area, and diagonal rows of pixels touching point to point have centers that lie on circles of latitude. With $`N`$ successive subdivisions of the diamonds into 4 sub-diamonds, we have a total of $`12\times 4^N`$ pixels of equal area covering the sphere. For example, WMAP uses $`N=9`$, where there are 3,145,728 diamond shaped pixels covering the sphere, each with dimensions of approximately $`0.11^{}\times 0.11^{}`$, which is adequate resolution to show the WMAP data which has an angular resolution of approximately $`0.3^{}`$. To study the topology properly, one needs pixels that are at least 2.5 times smaller than the smoothing length.

Each of these diamond pixels on the sphere has angles which would be complicated to calculate, so a version of CONTOUR2D for the sphere would seem to be complex and require many calculations to calculate the turning angles. However, we may simplify this greatly by using the topological invariance of $`g_{2D}`$ to projection (first from the sphere to the rhombic dodecahedron, then from the rhombic dodecahedron to a cube). This cube has six faces, which we will label T (top), N (north), E (east), S (south), W (west), B (bottom). Each diamond of the rhombic dodecahedron is oriented with its short diagonal along an edge of the cube, and each of the 12 rhombuses is mapped into two triangles that appear on adjacent faces of the cube. For example the rhombus TN is divided onto two triangles, one of which is mapped onto the T map on the cube and the other is mapped onto the N face of the cube. In fact, the letter designation of each rhombus tells us which of the two cubic faces it is mapped onto. For example, the SW rhombus is divided into two triangles which are mapped onto the S and W faces of the cube. The T face of the cube thus is a square which is divided into 4 right triangles whose hypotenuses form the sides of the square, and which meet in the center; these 4 triangles represent the upper halves of the four rhombuses TN, TE, TS, TW which meet at the top.

In such a projection we paint the hot regions with lead paint on the surface of the cube, and the 3D genus of these lead shapes is the same as if they were painted onto the sphere, because they have only been distorted in shape, not changed in topology. The total integrated Gaussian curvature over the cube, which totals $`4\pi `$ is concentrated entirely in 8 delta functions at the 8 vertices, where there is an angle deficit of $`\pi /2`$ at each, as 3 squares instead of four meet at each of these 8 vertices. Each of the 12 rhombuses of the original dodecahedron is mapped onto two right triangles meeting at their hypotenuses along an edge of the cube; these can be flattened out to make a square, with square pixels. So we have equivalently, 12 square maps of $`4^N`$ square pixels each. Thus, we can consider each rhombus as a square map with square pixels and we can calculate $`g_{2D}`$ for each using CONTOUR2D. Vertices that are within the rhombus are included. Vertices that are along an edge with the next rhombus are shared according to a prescription where the vertices along the top and right edges of the rhombus are assigned to it when laid out properly, and the bottom and left edges are assigned to adjacent rhombuses, since two rhombuses meet along each edge. That accounts for all pixel vertices except for those at the corners of the rhombuses. Each rhombus has 2 corners where 4 rhombuses meet (such as TN, TE, TS, TW meet at the top) and 2 corners where 3 rhombuses meet at a vertex (such as TN, TE, NE) which occurs at a vertex of the cube where there is a $`\pi /2`$ angle deficit in our cubic map and where 3 pixels only meet at a point. The corners where 4 rhombuses meet at a point are like ordinary vertices in the plane where 4 square pixels meet at a point, so this is handled with CONTOUR2D. The value of the contribution to the genus from that vertex must then be divided and shared equally among each of the rhombuses that meet there. The prescription for vertices where 3 pixels meet at a point (at the corner of the cube in our projection) is as follows:

| 0 hot 3 cold | genus = 0 the contour does not intersect the vertex |
| --- | --- |
| 1 hot 2 cold | genus = 1/4 added to hot rhombus, 0 added to others |
| 2 hot 1 cold | genus = 0 |
| 3 hot 0 cold | genus = 1/12 added to each of the three rhombuses. |

If we were to smooth the edges of the cube by sanding, the edges would turn into quarter cylinders (with zero Gaussian curvature) and the 8 vertices (where three squares meet at a point) turn into tiny octants of a sphere each with an integrated Gaussian curvature of $`\pi /2`$. The faces of the cube are flat so all the curvature is located in the 8 vertices. The edges of the three rhombuses meet at $`120^{}`$ angles in the spherical octant. So in the case 1 hot and 2 cold there is a turn of $`60^{}`$ there in the contour plus a curvature of $`30^{}`$ (equal to 1/3 of the integrated curvature over the octant of the sphere) so the genus contribution (turn plus parallel transport deviation = $`90^{}`$) should correspond to $`+1/4`$ (or $`\pi /2`$ divided by $`2\pi `$). This is then added to the hot rhombus as this is the hot spot we are circling. In the case 2 hot, 1 cold, there is a $`60^{}`$ turn, plus two regions of $`30^{}`$ curvature contributing to the parallel transport deviation, making a total contribution of 0 to the genus. In the case 3 hot, 0 cold, the vertex lies entirely within the hot region but it adds $`\pi /2`$ to the curvature integral within that region and this must be shared equally with all three rhombuses to add $`+1/12`$ to the genus in each. The total genus for the sphere $`g_{2D}`$ is calculated by adding the genus from all 12 rhombuses. From $`g_{2D}`$ we can calculate $`g_{2D,eff}`$ by subtracting $`2f`$. We can calculate $`f`$ as simply the fraction of the diamond shaped pixels that are hot, since these are all equal area. Thus, from the WMAP data we can calculate the 2D topology and compare it with that expected from a Gaussian random-phase distribution.

## 5 The Genus of the CMB, Measured by WMAP

The most direct data product from the WMAP team is the “internal linear combination” (ILC) map of the CMB, which is given in the HEALPix format. This ILC map uses the optimum linear combination of the skymaps at the different frequencies to remove the Galaxy and some other foregrounds (Bennett et al. 2003b). This is the primary map distributed as the best rendering of the CMB and has a resolution (beam width) of $`0.3^{}\times 0.3^{}`$. Though this map presents a few problems for directly computing the genus (which we will address shortly), it is worth checking the genus on this map directly.

The HEALPix projection provides a natural division of the dataset (sky) into 12 independent regions. We therefore measured the genus in each of these regions and computed the total for the whole sky, which we plot in Fig. 2. In this figure $`\nu `$ is computed in terms of the pixel-wise temperature mean and standard deviation measured over the whole sky (method 1). This method would be equivalent to the area fraction method described by equation (3) if the pixel histogram of temperature were strictly Gaussian. Fig. 3 is identical, except that $`\nu `$ is computed by area fraction in each individual HEALPix rhombus (method 2); see equation (3). Using the area fraction method measures more directly the random-phase nature of the distribution, since it separates that from any departures from a Gaussian histogram of the temperature (which could be measured directly). In method 1, we use $`\nu `$ defined in terms of the standard deviation in temperature over the whole sky; in method 2, we define $`\nu `$ by area fraction in each rhombus, thus treating each of the 12 rhombuses completely independently. The two methods provide a contrast of data treatment, and yet, as we shall see, give essentially identical results.

The inner errorbars in Figs. 2 and 3 are the standard deviation of the mean in the genus for the whole sky, estimated from those the 12 independent regions of the sky (adding the values of $`g_{2D,eff}`$ obtained from each). The standard one-sigma errorbars are shown as the inner errorbars. The outer errorbars show the $`95.4\%`$ (Gaussian two-sigma) confidence interval for a Student’s-$`t`$ variable with 12 degrees of freedom, which we will explain shortly.

The solid curves in Figs. 2 and 3 give the Gaussian random-phase genus curve, according to equation (2), with the best fit amplitude applied. The amplitudes, $`A`$, \[$`g_{2D,eff}=A\nu \mathrm{exp}(\nu ^2/2)`$\], in Figs. 2 and 3 are $`A=3432`$ and $`A=3657`$, respectively. These are by far the highest in any genus measurements of cosmological structure to date. As expected, the first value of $`A`$ is lower than the second, because the power at large scales (quadrupole and octopole) lifts some regions (rhombuses) to higher temperature, congealing some hot spots therein, but lowers other regions (rhombuses), congealing some cold spots therein. Method 2 sets a median in each rhombus separately and thus discovers more structures.

The Student’s-$`t`$ formulation is necessary since we are estimating the standard error with the distribution itself. Namely, if the true genus for the whole sky at each point on the curve were $`\mu `$ and the true standard deviation in the mean expected for that number of structures were $`\sigma `$ (we expect the distribution to be Gaussian for Gaussian random-phase fields \[Gott et al. 1990\]), then our estimator for $`\mu `$ is the arithmetic average estimated from the 12 rhombuses (the values from each rhombus are multiplied by twelve, then they are averaged), $`\overline{x}`$, and our estimator for $`\sigma `$ is $`s/\sqrt{n1}`$, where $`s`$ is the root-mean-square difference from the mean (as usual). We expect that $`\mu `$ is given exactly by equation (2) and that our best fit for the amplitude has negligible error compared to any one point on the curve. The expression $`(\overline{x}\mu )/(s/\sqrt{n1})`$ is defined to be a Student’s-$`t`$ variate (Lupton 1993).

On average, we would expect only one of the 21 points in each of the genus curves in Figs. 2 and 3 to miss the curve outside of the 95.4% confidence errorbar. In fact, not even one misses. From the properties of the Student’s-$`t`$ distribution with 12 d.o.f., one would expect to find that 33.7% ($`7`$) of the 21 points miss the curve by more than the inner errorbars. In both cases, the agreement is again better than expected, with only three misses in each figure. This better-than-expected agreement can easily arise by chance from correlations between the genus measurements at neighboring values of $`\nu `$, as shown by Colley (1997) (these correlations mean that the 21 data points are not completely independent—neighboring genus values measure many of the same structures). Nonetheless, the genus of temperature fluctuations in the CMB, as observed by WMAP, is certainly consistent with the Gaussian random-phase hypothesis.

## 6 Genus on a Stereographic Map Projection

In the ILC HEALPix data release, there are a couple issues that must be addressed. The first issue is foreground emission. Bennett et al. (2003b) discuss in quite a bit of detail the foregrounds in the WMAP results. They provide a very useful sky mask that should be applied to avoid foreground contaminants, the principal one of which is the Galaxy. In a second study, we have therefore used their Kp0 (most conservative) sky mask, which also removes some 200 point sources. Furthermore, we have restricted ourselves to Galactic latitudes $`|b|>18^{}`$. Since the mask must remove much of the area near the Galactic Equator, the HEALPix projection, with its 4 equatorial diamonds, ceases to be ideal. Instead, we re-project the map, using the stereographic projection for both the Northern and Southern Galactic Hemispheres (as in Fig. 1, in which the $`|b|>18^{}`$ cut occurs at 72.7% of the total radius in each hemisphere).

We can calculate $`g_{2D}`$ from the north and south galactic caps by using stereographic map projections of the caps. Stereographic map projections of a hemisphere are conformal and preserve angles. A stereographic projection of the north galactic cap is a circle, in which azimuthal angle in the map is the galactic longitude in the cap. Radius in the map, $`r=2\mathrm{tan}[(90^{}b)/2]`$ (See Fig. 1). The outside boundary is a circle of radius $`r_b=2\mathrm{tan}[(90^{}18^{})/2]`$. The contours on the sphere can be approximated as spherical polygons—geodesic segments joined at vertices where turning on the sphere occurs. The stereographic projection preserves angles, so the turning that occurs at vertices is mapped properly. But, geodesics (great circles) are mapped into arcs of circles rather than straight lines, which turn on the map but not on the sphere so the total genus $`g_{2D}`$ on the sphere is not automatically calculable from the total turning that comes from circling the density contour on the map.

So instead of calculating $`g_{2D}`$ for the sphere from the stereographic map directly, which would be complicated, let us instead calculate the relation between $`g_{2D,eff}`$ on a spherical cap, and $`g_{2D,\text{map}}`$ measured by the CONTOUR2D program on the flat stereographic map of the cap. \[Of course once we have $`g_{2D,eff}`$ we can calculate the contribution to $`g_{2D}`$ due to the structures in the north galactic cap by adding $`2f\text{Area}_{\text{cap}}/(4\pi )`$\]. But it is really $`g_{2D,eff}`$ we are interested in anyway because it is what will be compared to the theoretical formula. Thus, we can skip a step by going directly after $`g_{2D,eff}`$.

To begin, let us excise the spherical cap from the sphere by cutting it off and separating it from the rest of the sphere. This will add boundary terms. Each hot spot (think curved lead sheet) that intersects the boundary of the cap will be closed off by a curve along the boundary that will make it an isolated hot spot. Think about these hot spots in 3D when the cap is excised from the rest of the sphere: any hot spot that hits the boundary will become isolated, floating in space, because beyond the boundary there is nothing. This spherical isolated cap will have a value of $`g_{2D,\text{excised cap}}=g_{2D,\text{excised map}}`$; the total genus of the isolated flat stereographic map sitting isolated in the plane—think of flat lead sheets lying in a plane; because one can be continuously deformed into the other. We can calculate the terms contributing to $`g_{2D,\text{excised cap}}`$:

$$\begin{array}{cc}\hfill g_{2D,\text{excised cap}}=& (2\pi )^1(\text{turns on spherical cap})\hfill \\ & +f\text{Area}_{\text{cap}}/(2\pi )\hfill \\ & +(2\pi )^1(\text{turns at cap boundary})\hfill \\ & +f_b[1\text{Area}_{\text{cap}}/(2\pi )].\hfill \end{array}$$
(6)
The first term, $`(2\pi )^1(\text{turns on spherical cap})`$, is equal to $`g_{2D,eff,\text{cap}}`$ for the interior of the spherical cap. It is the total turns made driving on the sphere around all the contours within the spherical cap, and excludes any boundary terms produced by excising the spherical cap. Now $`g_{2D,eff,\text{cap}}\nu \mathrm{exp}(\nu ^2/2)`$, so it is what we want to calculate to compare with the theoretical genus formula. The term $`f\text{Area}_{\text{cap}}/(2\pi )`$, where $`\text{Area}_{\text{cap}}`$ is measured in steradians, is the integrated Gaussian curvature within the hot spots in the cap divided by $`2\pi `$, which contributes to the genus by the 3D argument presented above. The term $`(2\pi )^1(\text{turns at cap boundary})`$ is the sum of the turns that are made when you are driving along a contour and you encounter the cap boundary and then must turn a sharp corner there to continue driving along the cap boundary to complete circling that hot spot as an isolated region. In general, if along the boundary there are $`N`$ hot segments, there will be $`2N`$ such turns, as we enter and then exit each of these N hot segments along the boundary. Each of these $`2N`$ turns will be some angle between 0 and $`2\pi `$. The final term, $`f_b[1\text{Area}_{\text{cap}}/(2\pi )]`$, represents the total turning on the sphere that occurs driving along the $`N`$ hot segments of the boundary of the cap. If the entire boundary were hot then driving around it would have a total turn of $`2\pi \text{Area}_{\text{cap}}`$. (For example, if this were a hemispherical cap the area of the cap would be $`2\pi `$ and the total turn circling the boundary would be zero, because the boundary is the equator and it is a geodesic so no turning is required to circle the boundary). As the area of the cap goes to zero, the curvature within the cap becomes unimportant to parallel transport and the total turn required to circle the boundary approaches $`2\pi `$. In general, the angle deflection suffered by parallel transport around the closed boundary curve is equal to the integral of the Gaussian curvature inside, which in this case is equal to the area of the cap. When we have completed the circuit, we have returned to where we have started so we have suffered a total rotation of $`2\pi `$, so the turning we do on the sphere circling the boundary plus the area of the cap is $`2\pi `$. Thus, turning = $`2\pi \text{Area}_{\text{cap}}`$. The number $`f_b`$ is the fraction of the boundary that is hot, and so it is this fraction of the boundary that we drive around. Dividing this by $`2\pi `$, gives the contribution to the genus derived by driving along the hot segments of the boundary: $`f_b[1\text{Area}_{\text{cap}}/(2\pi )]`$. Using the Copernican Principle, one would not expect the boundary to be special, so on average the fraction of the cap boundary that was hot should be equal to the fraction of the sphere that is in hot spots. On average then, we expect $`f_b=f`$. Substituting we find:

$$g_{2D,\text{excised cap}}=g_{2D,eff,\text{cap}}+\frac{f\text{Area}_{\text{cap}}}{2\pi }+\frac{1}{2\pi }(\text{turns at cap boundary})+f\left(1\frac{\text{Area}_{\text{cap}}}{2\pi }\right).$$
(7)
$$g_{2D,\text{excised cap}}=g_{2D,eff,\text{cap}}+f+\frac{1}{2\pi }(\text{turns at cap boundary}).$$
(8)
Now let us measure the genus on the excised stereographic map ($`g_{2D,\text{excised map}}`$), recalling that the excised stereographic map is topologically equivalent to the excised cap.

$$g_{2D,\text{excised map}}=\frac{1}{2\pi }(\text{turns on stereographic map})+\frac{1}{2\pi }(\text{turns at map boundary})+f_b.$$
(9)
The first term, $`(2\pi )^1(\text{turns on stereographic map})=g_{2D,\text{map}}`$, is the quantity measured by CONTOUR2D on the flat map. CONTOUR2D examines only the Cartesian vertices within the stereographic map region and does not include any boundary effects. There is no $`f\text{Area}`$ term because the map is flat and there is zero curvature inside the hot spots. The next term, $`(2\pi )^1(\text{turns at map boundary})`$, adds the turns that are taken when the contour hits the outer circular boundary of the map. If we use a stereographic projection, which is conformal and preserves angles, each and every turn at the boundary in the map will be equal that encountered on the sphere, so:

$$\frac{1}{2\pi }(\text{turns at map boundary})=\frac{1}{2\pi }(\text{turns at cap boundary}).$$
(10)
The final term, $`f_b`$, is the total turn taken on the $`N`$ hot segments of the boundary of the map divided by $`2\pi `$. If we circled the entire boundary we would have to turn a total of $`2\pi `$, because the map is flat and Euclidean geometry applies. The $`N`$ hot segments cover a fraction $`f_b`$ of the boundary (galactic longitude is mapped onto azimuth in the stereographic map, so the fraction of the map boundary occupied by the $`N`$ hot segments is exactly the same as on the sphere). Thus, the total turn is $`2\pi f_b`$, and dividing this by $`2\pi `$ gives $`f_b`$. As we have remarked above, if the boundary is not special we expect $`f_b=f`$. Since the map is excised and sits alone on the plane, these outer boundary segments must be included to create isolated hot spots in the plane out of the hot regions that hit the boundary. Equating $`g_{2D,\text{excised cap}}`$ and $`g_{2D,\text{excised map}}`$ (because the cap can be deformed into the map), and substituting from the above equations we find:

$$g_{2D,\text{excised cap}}=g_{2D,\text{excised map}}.$$
(11)
$$g_{2D,eff,\text{cap}}+f+\frac{1}{2\pi }(\text{turns at cap boundary})=g_{2D,\text{map}}+\frac{1}{2\pi }(\text{turns at map boundary})+f.$$
(12)
$$g_{2D,eff,\text{cap}}=g_{2D,\text{map}}.$$
(13)
So “what you see is what you get.” If we measure the 2D topology in the interior of the stereographic map using CONTOUR2D (which ignores the boundary) our result will equal, on average, $`g_{2D,eff,\text{cap}}`$ which is what we were looking for. This works for any radius spherical cap, since terms proportional to the area of the cap cancel out. And since the map is conformal, the turns at the boundary are equal in both cases and cancel out as well. Finally, since we expect the boundary to be typical, its coverage with hot spots should be equal to the fraction $`f`$ of the sphere covered by hot spots, and this causes other terms to cancel as well, leaving the above simple result.

## 7 Genus Topology on Stereographic WMAP Results

We may now compute the 2-D genus on a stereographic re-projection of the WMAP data, in which we address the issues of foreground contamination and noise variability arising from the scanning strategy of the WMAP spacecraft.

First, to alleviate contamination from the Galaxy, we completely exclude the equator region up to $`|b|<18^{}`$. The Kp0 mask does include removal of local regions for point source avoidance, and some additional Galaxy masking above the $`|b|>18^{}`$ cut, particularly near the Galactic Center, the Gum Nebula and Orion Complex. Any vertex in the stereographic map projection that falls within the mask is simply ignored by CONTOUR2D.

We proceed by dividing the stereographic projection into octants which correspond to octants on the sphere. These octants form 8 independent regions in which to measure the genus. A slight adjustment for the actual solid angle in each octant is made, because the Kp0 mask changes the solid angle of each octant very slightly (a few percent variation). We normalize the genus to the mean solid angle of the masked octants, and proceed treating each masked octant as one of 8 independent measurements of the genus per average masked octant.

In carrying out the re-projection to the stereographic map, we have applied a slight smoothing to a total smoothing FWHM of $`0.35^{}`$. Figs. 4 and 5 show the genus for the masked stereographic projection, computed by methods 1 and 2 (same methods as in the HEALPix genus computation). The errorbars are computed now for 8 octants instead of 12 diamonds, which changes somewhat the relative size of the 95.4% errorbars and the direct one-sigma errorbars. The best fit theoretical genus curve (equation 2) is plotted in each figure, with respective amplitudes of $`A=2057`$ and $`A=2186`$ (as with the HEALPix genus curves, the first value is slightly less than the second, as expected). The amplitudes have come down somewhat from the HEALPix maps due to the exclusion of the equatorial band (so the total genus covers only a fraction of the whole sphere), the additional masked pixels in the Kp0 prescription, and the additional smoothing applied in the stereographic projection.

In each of the figures, the agreement with the theoretical genus curve is striking. Quantitatively, we would expect on average that one of the 21 points sampling the curve would miss the theoretical fit outside is 95.4% errorbar. With method 1 (Fig. 4), no points miss the theoretical curve outside their 95.4% errorbar, while with method 2 (Fig. 5), two points miss by a small margin. The one-sigma (inner errorbars) for the 8 degrees of freedom of the octants reside at the 65.3% confidence interval, so on average, one would expect approximately 7 of the 21 points on each of Fig. 4 and Fig. 5 to miss the theoretical curve outside their one-sigma errorbars. Using method 1, we find that 3 points miss at this margin, and 9 miss using method 2. Though the method 2 results show slightly more than the average number of misses expected at both the one-sigma and 95.4% level, the excess is certainly within the uncertainty arising from small number statistics in each case.

The second contaminant in the ILC HEALPix map is that the noise behavior is very complicated, and differs substantially pixel to pixel. First, the beams at each frequency in the combination have different widths. Second, the scan method of the WMAP satellite yields many more samples at the ecliptic poles than at the ecliptic equator (see Bennett et al. 2003a, Fig. 3). To alleviate some of these problems, we have downloaded from the WMAP website maps in the same 5 frequencies as used in the ILC map, but ones that have been pre-smoothed to the same effective beam-width ($`0.82^{}`$). These are then added with the same weights as used in producing the ILC map. At cost of number of structures, these maps significantly enhance the signal-to-noise of remaining structures, thus reducing the significance of the noise differences arising from the scanning strategy. This is, therefore our most conservative map.

In projecting these maps into stereographic maps, we have smoothed to $`0.9^{}`$, slightly more than the $`0.82^{}`$ smoothing of the pre-smoothed maps. Figs. 6 and 7 show the genus measured in these stereographic projections. The best-fit amplitudes have dropped to $`A=772`$ for method 1 and $`A=831`$ for method 2, due to the smoothing.

From these final stereographic projections, the genus measured at the 21 points should also be expected to miss the theoretical curve outside of the 95.4% errorbars once on each figure. What we find is that with method 1 (Fig. 6), no points miss the curve outside the 95.4% errorbar, while with method 2, (Fig. 7), one point misses by a small amount. The points miss the curve outside their one-sigma errorbars 2 times for method 1, and 4 times for method 2, again a bit better than expected. The final figure, Fig. 7, represents our most conservative estimate of the true genus available in the WMAP distribution, and shows excellent agreement with the Gaussian random-phase hypothesis.

## 8 Discussion

As we have shown, topologically, within the errors, the CMB measured by WMAP appears to be consistent with Gaussian random-phase. Indeed, the overall visual impression is strikingly random-phase. The stereographic projection is conformal, so shapes are preserved locally, though the scale grows as one approaches the the outer boundary (by a factor of two from the pole to the equator, as seen in Fig. 1). We may ask if there are any features visible in the map at all that look non-random-phase in any way.

Our choice of map projection has a miscellaneous advantage over the Molleweide projection used by the WMAP team, which is that it shows the Galactic polar regions very well. In the south polar region, there is one possibly non-random feature, a fairly narrow red (hot) feature stretching from approximately the 8:30 position toward the 2:30 position curving like a bow below the South Galactic Pole. This feature is more visually apparent when viewed in landscape orientation. It has the appearance of a scanning error approximately along a great circle, but is not coincident with any particular locus in the WMAP scanning geometry.

We have considered several possibilities to explain such a feature. First, there is the simple phenomenon of the “canals” on Mars, the linear features reported in naked-eye observations of the red planet in the late $`19^{\text{th}}`$ century. We now know that there are no such canals, and that our eyes’ highly evolved ability to detect linear features or edges was simply creating an illusion. Most likely the same thing is going on here, particularly since the feature is more apparent when the figure is viewed in landscape rather than in portrait orientation. In fact, if one had a true Gaussian random-phase distribution and there were no features that looked unusual to the eye, that would be unusual. So this feature may be merely coincidental.

What are the other possibilities? It is not coincident with the Magellanic Stream measured seen in 21cm (Brooks 2000). Furthermore, we would expect the stream to have similar microwave “color” to the gas in the Galaxy, and hence be subtracted about as well as the Galaxy is in the ILC map. To test this, we constructed a three color map, using the K band as red, Q as green and W as blue, and found that while the foregrounds were distinctly pink to orange to yellow, the CMB was quite grey, and so was the feature.

Another possibility that comes to mind is cosmic strings. Cosmic strings, however, do not amplify temperature along a line. Instead, a moving cosmic string produces a cliff, hotter on one side (the trailing side) than the other side (the leading side). Thus a hot stripe would require two parallel cosmic strings, having just passed each other (an interesting theoretical possibility \[Gott 1991\], but not a likely one, since we are unlikely to find the two most prominent cosmic strings within the horizon in such an alignment by chance). Also, there is no obvious evidence for cosmic string cliffs anywhere else in the WMAP data, which we have tested for by making gradient maps as suggested by Gott et al. (1990).

The feature could also be due to the integrated Sachs-Wolfe effect from a great wall seen edge-on, but in looking at several galaxy catalogs and X-ray maps, we could detect no obvious feature with appropriate geometry.

This list of possibilities is by no means exhaustive, and, in the end, it may be most likely that we are simply falling prey to “Martian canals” in an overall Gaussian random-phase map. Since there are countless possibilities for what one may see and deem non-random (circles, squares, triangles, even letters or words), the a posteriori statistical significance of any such weak feature is difficult to estimate.

The overall consistency with of the genus topology in the WMAP data with that predicted a priori from a Gaussian random-phase distribution is quite spectacular. WMAP’s unprecedented combination of angular resolution and all-sky coverage of the Cosmic Microwave Background has presented by far the best confirmation to date of the standard inflationary prediction of Gaussian random-phase initial conditions.

In standard inflation, the structures seen in the CMB are fossil remains of random quantum fluctuations occurring just $`10^{33}`$ seconds after the Big Bang. In galaxy clustering surveys, these fossils have become somewhat distorted due to non-linear effects and biased galaxy formation. In the CMB, however, we are seeing these fluctuations while they are still in the linear regime, making the WMAP measurements of the CMB a very attractive dataset for studying these primordial perturbations.

With the hundreds of structures observed even at the $`0.9^{}`$ scale, the WMAP results have provided the largest and best dataset for studying the primordial structures. The genus topology test, in particular, tests the dataset for the standard inflationary prediction that the fluctuations should derive from a Gaussian random-phase distribution. All of our tests have confirmed this prediction dramatically.

## 9 Acknowledgments

The WMAP datasets were developed by Princeton University and NASA’s Goddard Space Flight Center, and graciously provided to the public at http://map.gsfc.nasa.gov/. JRG thanks the National Science Foundation for its support under grant AST 99-00772.


text:
# Rigid vortices in MgB2

## Abstract

Magnetic relaxation of high-pressure synthesized MgB2 bulks with different thickness is investigated. It is found that the superconducting dia-magnetic moment depends on time in a logarithmic way; the flux-creep activation energy decreases linearly with the current density (as expected by Kim-Anderson model); and the activation energy increases linearly with the thickness of sample when it is thinner than about 1 mm. These features suggest that the vortices in the MgB2 are rather rigid, and the pinning and creep can be well described by Kim-Anderson model.

The recently discovered type-II superconductor magnesium diboride (MgB2) Nagamatsu ; Finnemore  exhibits noticeably different magnetic properties and flux dynamics as compared with the high-temperature cuprate superconductors and conventional superconductors. As a $`s`$-wave superconductor Kohen ; Plecenik , it shows odd behaviors in the temperature dependence of the lower critical field $`H_{c1}(T)`$ and the penetration depth $`\lambda (T)`$ Li ; the irreversibility field $`H_{\mathrm{irr}}(T)`$ is remarkably lower than the upper critical field $`H_{c2}(T)`$ near zero temperature; the evidently low magnetic relaxation rate depends weakly on the temperature in the low and intermediate temperature region Wen1 ; Wen2 , which makes the magnetic moment rather stable Thompson , and it is interesting that the relaxation rates of the high-pressure synthesized bulk and the film Wen1 ; Zhao1  differ by about one order of magnitude, implying that the flux-creep activation energy might depend on the thickness of the sample.

In view of these, the systematic magnetic relaxation measurements of MgB2 bulks in the shape of cylinder are carried out with a superconducting quantum interference device (SQUID) in this work, where the height (thickness) of the cylinder is decreased at each time. It is found that the superconducting dia-magnetic moment depends on time in a logarithmic way, and the exponent $`\mu `$ in the universal equation of activation energy Blatter ; Yeshurun  is roughly $`1`$ as expected by the Kim-Anderson model for a vortex segment with a fixed length.

The high-pressure synthesized bulk sample of MgB2, which is in the shape of a cylinder with the diameter of 1.5 mm and the original height of 3 mm, was fabricated at 950C for 0.5 hour under the pressure of 6 GPa. The synthesis procedures were published elsewhereRen . The X-ray diffraction (XRD) analysis shows that it is nearly in a single phase with the second phase less than 1 wt.%. The height of the cylinder is decreased by sawing and polishing after each experiment, from the original value 3 mm to 0.16 mm at last.

The measurements of magnetic relaxation are carried out by a SQUID (Quantum Design MPMS, 5.5 T), with the temperature varying from 2 K to 32 K and in the applied magnetic field of 0.5 T , 1 T and 1.5 T, respectively. The axis of the cylinder is aligned with the field direction during the magnetic measurements.

The data of the magnetic relaxation for the sample with thickness of 1.15 mm and in the magnetic field of 0.5 Tesla, are plotted in the inset of Fig. 1. It shows that the absolute value of the dia-magnetic moment $`M`$ decreases in a logarithmic way with time, i. e., $`M\mathrm{ln}t`$.

We start from the model of thermally activated flux motion (TAFM) in which the dissipation is described by

$$E=v_0B\mathrm{exp}\left(\frac{U(j_s,T)}{k_\mathrm{B}T}\right),$$
(1)
where $`j_s`$ is the superconducting current density, $`E`$ is the electric field due to TAFM over the activation energy $`U(j_s,T)`$, $`v_0`$ is the average velocity of the flux motion, $`B`$ is the magnetic induction. For the activation energy $`U(j_s,T)`$, many models have been proposed Anderson ; CP ; Fisher ; Zeldov , among which the collective pinning CP  and vortex glass Fisher  models have received most attention in the study on the flux dynamics of cuprate superconductors. To guarantee that $`U(j_s,T)|_{j_sj_c}=0`$, the following interpolation was proposed Malozemoff  and has been widely used:

$$U(j_s,T)=\frac{U_c(T)}{\mu }\left[\left(\frac{j_c(T)}{j_s}\right)^\mu 1\right]$$
(2)
where $`U_c`$ is the activation energy of vortices. Although Eq. (2) was proposed based on the flux dynamics in cuprate superconductors, it can be used as a general equation which describes the flux dynamics under specific situation by choosing different $`\mu `$ values. For example the linear $`U(j)`$ dependence requested by the Kim-Anderson model Anderson  to describe the creep of a vortex segment with fixed length can be obtained by choosing $`\mu =1`$. From the above two equations, one can derive the superconducting current density,

$$j(t)j_c\left[1+\frac{\mu k_\mathrm{B}T}{U_c}\mathrm{ln}\left(1+\frac{t}{t_0}\right)\right]^{\frac{1}{\mu }}.$$
(3)
When choosing $`\mu =1`$, one can easily see that $`j`$ depends on time in a logarithmic way which is consistent with the data shown in the inset of Fig. 1. This feature may suggest that the flux dynamics can be described by the Kim-Anderson model rather than the collective pinning model as suggested recently Qin . The temperature dependence of the relaxation rate $`S=d\mathrm{ln}M/d\mathrm{ln}t`$ is then derived, as shown in Fig. 1. The extrapolation of relaxation rate does not approach zero as $`T0`$, which has been interpreted by either quantum tunneling of vortices or the two-gap feature of MgB2 Wen1 ; Wen2 .

In order to get the activation energy, we use the relation Wen3

$$T/S=U_c/k_\mathrm{B}+\mu T\mathrm{ln}(1+t/t_0),$$
(4)
where $`\mathrm{ln}(1+t/t_0)`$ can be treated as a constant within a fixed time window. The inset of Fig. 2 shows a big hump in each curve of $`T/S`$ versus $`T`$. It is shown in Fig. 1 that the relaxation rate $`S`$ keeps flat up to a temperature of about 15 K. Therefore the rising part of the big hump in the low temperature region is induced by the weak temperature dependence of $`S`$. When $`T`$ is above 15 K, $`T/S`$ decreases with temperature which may be induced by the decreasing of $`U_c(T)`$ or by a negative $`\mu `$ value. One can have an approximate value for $`U_c`$(15 K), which is about 7000 K for 0.5 T. This value is about 10 times greater than that in the high-temperature cuprate superconductors. We will show below that this larger activation energy is induced by the strong pinning to a very long rigid vortex.

One important factor characterizing the vortex dynamics is the $`U(j)`$ dependence. For collective pinning, $`U(j)`$ is non-linear and the $`\mu `$ value in Eq. (2) is positive. However, as forementioned, for a vortex segment with fixed length, Kim-Anderson model proposes a linear $`U(j)`$ relation. In our present work, Maley’s method Yeshurun ; Maley  is used to determine the dependence of the activation energy on the current density. The results are shown in Fig. 3. The linear $`U(j)`$ relation suggests an Kim-Anderson type effective pinning barrier $`U(j)=U_c(1j/j_c)`$ with $`\mu =1`$. At low temperatures, the curves deviate from linearity, which indicates that other mechanism may set in, leading to the rather stable $`S`$ value in low temperature region. It is interesting to note that in high-temperature cuprate superconductors, such as YBa2Cu3O7, $`U(j)`$ diverges in the small current limit which is resulted by the elastic property of the vortices.

Based on the assumption $`U_c(T)T/S`$, we calculate the values of $`U_c`$. Fig. 2 shows that the values of $`U_c`$ is proportional to the thickness of sample and saturates at about 1 mm. It also depends on the disorder strength which is different between bulks and films Zhao2 . In the scenario of weak collective pinning where the pinning barrier is caused by the weak point disorder, vortices will break up into segments Blatter  due to the competition between the elastic and pinning energy, the longitudinal size of which is the collective pinning length $`L_c(\epsilon _0^2\xi ^2/\gamma )^{1/3}`$, with $`\epsilon _0`$ the basic energy scale, $`\xi `$ the coherence length, and $`\gamma `$ a parameter of disorder strength. Thus the linear dependence of $`U_c`$ on the sample thickness suggests that the flux line in MgB2 bulks will behave as an individual segment if its length is below $`L_c1`$ mm. Combining with the significant large value of $`U_c`$ ($`10^3`$ K), we conclude that the vortices are remarkably rigid over the whole sample.

For a long and rigid vortex, the intrinsic pinning energy $`U_c`$ is large, therefore the thermal activation barrier is high, vortices may move to the next pinning site by quantum tunneling. This naturally interprets the rather stable $`S`$ in the temperature region up to 15 K ($`0.4T_c`$) and also the large gap between $`H_{c2}(0)`$ and $`H_{\mathrm{irr}}(0)`$. However, the big gap between $`H_{c2}(0)`$ and $`H_{\mathrm{irr}}(0)`$ may be explained by an alternative model based on two-gap feature $`\mathrm{\Delta }_1(0)=1.8`$ meV and $`\mathrm{\Delta }_2(0)=6.8`$ meV Choi , that the upper critical field $`H_{c2}(T)`$ is determined by $`\mathrm{\Delta }_2(T)`$, while the irreversibility field $`H_{\mathrm{irr}}(T)`$ may include the contribution of both gaps. The thickness dependence of $`U_c`$ even holds below 15 K, implying that some further research are needed to clarify whether the quantum tunneling or the two-gap feature dominates.

In conclusion, the magnetic relaxation of high-pressure synthesized MgB2 cylinder is measured by SQUID in a series of experiments, with the height of sample decreased each time. The linear $`M`$-$`\mathrm{ln}t`$, $`U`$-$`j`$ relations and the thickness dependence of the activation energy are found. These features imply a Kim-Anderson type flux creep and the rather rigid vortices with the collective pinning length $`L_c`$ about 1 mm in bulk MgB2, in contrast to the cuprate superconductors and the previous results of collective-pinning. The thickness dependence of $`U_c`$ even at low temperatures also suggest the quantum tunneling to be still as an open question.

This work is financially supported by the National Natural Science Foundation of China (Grant No. NSFC 19825111) and the Ministry of Science and Technology of China (Grant No. NKBRSF-G1999064602).


text:
# Time-scale dependence of correlations among foreign currencies

## 1 Introduction

It is widely known that foreign exchange rates sometimes show very violent fluctuations characterized by power law distributions \[T. Mizuno\]\[K. Matia\]. A very important report appeared recently that arbitrage chances exist by considering three currency markets simultaneously \[Y. Aiba\]. This result is due to the fact that the distribution of return of rotation transaction of three currencies show a similar fat-tail. As statistical laws among foreign exchange rates in the short time scale are not known sufficiently, we investigate correlations among foreign exchange rates in this paper.

Table.1 The set of Daily data for 25 kinds of foreign exchange we examined.

| Base Currency : | U.S. Dollar |  |  |
| --- | --- | --- | --- |
| Target Currencies | : |  |  |
| Australian Dollar, | Brazilian Real, | British Pound, | Canadian Dollar, |
| Chilean Peso, | Colombian Peso, | Czech Koruna, | Euro, |
| Hungarian Forint, | Indonesian Rupiah, | Japanese Yen, | Mexican Peso, |
| New Zesland Dollar, | Norwegian Kroner, | Peruvian New Sole, | Philippines Peso, |
| Polish Zloty, | Russian Ruble, | Singapre Dollar, | Slovakian Koruna, |
| South African Rand, | South Korean Won, | Swedish Krona, | Swiss Franc, |
| Thai Baht |  |  |  |

Correlations among financial stocks in markets have been actively discussed. The strength of synchronization among stocks was analyzed by using the ultrametric spaces and they discussed about the portfolio of stock \[R. N. Mantegna, H. E. Stanley\]. The interaction of companies was investigated by analyzing return of many stocks and a directed network of influence among companies was defined \[L. Kullmann\]. It will be noticed in the present paper that the correlation among foreign exchange rates resembles the stock case.

In the following section we discuss the cross-correlation between exchange rates with no time difference to show the relation between synchronization of exchange rate. Then, we show the maximum correlation between exchange rates is observed with nonzero time shift, namely, the direction of influence is discussed.

## 2 The correlation among foreign exchange rates

We first analyze a set of daily data provided by Exchange Rate Service \[ERS\] for 25 exchange rates for about 3 years from January ‘99 to August ‘01 as listed in Table 1. We first estimate cross-correlation functions for these exchange rates measured by USD (United States Dollar). The largest correlation value (=0.95) is observed for a pair of CHF(Swiss Franc)/USD and EUR(Euro)/USD. As demonstrated in Fig.1 it is evident that these exchange rates are remarkably synchronized. There are cases with negative correlation values as found for the pair of MXP(Mexican-Peso)/USD and CHF(Swiss Franc)/USD whose correlation value is —0.23. It is known that correlations between geographically closer currencies tend to have larger correlation and there exist key currencies in each area, for example, Euro for West Europe, Yen for Asia, Hungarian Forint for East Europe and Australian Dollar for Oceania \[H. Takayasu and M. Takayasu\].

Although we can find such large correlations among currencies in daily data, we can expect low correlations in short time scale as it is common that dealers in major banks tend to work with a single foreign exchange market. In order to clarify this tendency we examine tick-by-tick data provided by Reuters for about 4 months from March ‘02 to July ‘02. In Fig.2 we plot the correlation value as a function of coarse-grained time-scale for a pair of CHF/USD and EUR/USD. From this figure

we notice that the correlation vanishes if we observe the high-resolution data with the precision of seconds and the correlation value is about 0.5 in the time scale of 5 minutes (300 sec.). From these results it is understood that these two currency markets, CHF/USD and EUR/USD, are working independently in very short time-scale.

## 3 The cross-correlation with a time shift

In order to clarify the nature of short time interaction among currencies we calculate the cross-correlation with a time shift, namely, we observe the correlation of these two markets with a time difference by the following equation,

$$C\left(dt\right)=\frac{dp_A\left(t\right)dp_B\left(t+dt\right)dp_A\left(t\right)dp_B\left(t+dt\right)}{\sigma _A\sigma _B},$$
(1)
where, $`dp_A\left(t\right)`$ is the rate change in the market $`A`$ at time $`t`$, $`dp_B\left(t+dt\right)`$ is the rate change in the market $`B`$ at time $`t+dt`$, $`\sigma `$ is the standard deviation of rate changes in each market. In Fig.3, we show the correlation value between CHF/USD at time $`t+dt`$ and EUR/USD at time $`t`$ as a function of time difference $`dt`$. Here, we show two plots for different coarse-graining time-scales, 60 sec. and 120 sec.. In both cases it is found that the largest value of correlation is observed around $`dt=10`$ seconds, which implies that in an average sense the EUR/USD market is going about 10 second ahead and the CHF/USD market is following it.

## 4 Currency correlation in short time scale

Here, we discuss about value of currency correlation in each foreign exchange market in a short time scale. From Fig.2, it is noticed that the correlation between foreign exchanges is very low in a short time scale. Namely, each exchange rate is changing rather independently. In order to clarify this property, we analyze the exchange rate of Yen-Dollar by analyzing a set of tick-by-tick data provided by CQG for about 2 years from February ‘99 to March ‘02.

We introduce two definitions of Yen-Dollar rate: One definition of JPY/USD is the usual transaction rate and the other JPY/USD is defined through Euro as {EUR/USD}$`\times `${JPY/EUR}. Here, all exchange rates are given by the middle rate (=(Bid rate + Ask rate)/2). In Fig.4 we plot the cross-correlation value of these exchange rates in different time scales. The correlation value is not unity in the time scale less than about 1 hour. This result clearly shows that the value of a currency differs in different markets in the time scale less than an hour. This is a kind of self-contradiction of markets causing the occurrence of triangular arbitrage opportunity as shown in Fig.5 \[Y. Aiba\].

## 5 Discussion

We have clarified the detail properties of correlation among foreign exchanges; the short time correlation is generally very small even between the pair of currencies showing large correlation in daily data. Within the time scale of a minute we can observe the direction of influence from one currency market to others. In very short time scale we can find contradiction of exchange rates between JPY/USD and {EUR/USD}$`\times `${JPY/EUR}. For example, if you observe carefully the two rates, JPY/USD and {EUR/USD}$`\times `${JPY/EUR}, then you can buy Yen cheaper in one market and can sell it with higher rate in the other market even taking into account the effect of the spread (=Ask rate – Bid rate) that is about 0.05%. Although, no time lag is considered regarding the actual transactions, it is now clear how triangular arbitrage opportunity appears.

Acknowledgement

We would like to thank Mr. Hiroyuki Moriya of Oxford Financial Education for providing us with the CQG data, Prof. Tohru Nakano and Prof. Mitsuo Kono for stimulus discussions.

References

Mizuno, T., Kurihara, S., Takayasu, M., Takayasu, H., Analysis of high-resolution foreign exchange data of USD-JPY for 13 years, Physica A, in press.

Matia, K., et al (2002) Non-Levy Distribution of Commodity Price Fluctuations, Phys. Rev. E 66, 045103.

Aiba, Y., et al (2002) Triangular arbitrage as an interaction among foreign exchange rates, Physica A 310, 467-479.

Mantegna, R. N., Stanley, H. E., (2000) An Introduction to Econophysics: Correlation and Complexity in Finance, Cambridge University Press, Cambridge.

Kullmann, L., Kertesz, J., Kaski, K., (2002) Time-dependent cross-correlations between different stock returns: A directed network of influence, Phys. Rev. E 66, 026125.

Exchange Rate Service: http://pacific.commerce.ubc.ca/xr/ .

Takayasu, H., Takayasu, M., (2001) Econophysics-Toward Scientific Reconstruction of Economy, Nikkei, Tokyo. (in Japanese)

## Index


text:
# Local and Reversible Change of the Reconstruction on Ge(001) Surface between c(4×2) and p(2×2) by Scanning Tunneling Microscopy

## Abstract

The reconstruction on Ge(001) surface is locally and reversibly changed between c(4$`\times `$2) and p(2$`\times `$2) by controlling the bias voltage of a scanning tunneling microscope (STM) at 80 K. It is c(4$`\times `$2) with the sample bias voltage $`V_b`$ \- 0.7 V. This structure can be kept with $`V_b`$ 0.6 V. When $`V_b`$ is higher than 0.8 V during the scanning, the structure changes to p(2$`\times `$2). This structure is then maintained with $`V_b`$ \- 0.6 V. The observed local change of the reconstruction with hysteresis is ascribed to inelastic scattering during the electron tunneling in the electric field under the STM-tip.

Local modification of surface structure in atomic scale has been demonstrated using scanning tunneling microscopes (STM’s). Surface adatoms were reversibly manipulated first on a Pt surface at low temperature Eigler . On the other hand on semiconductor surfaces, the modifications by temporal voltage pulses, such as the removals of surface atoms Lyo  and adsorbed hydrogenBecker ; Lyding , have been performed. Up to now, however, there has been no report on the reversible modification of local surface reconstructions by using STM.

Such modification can be expected on the ground state c(4$`\times `$2) reconstructed surface of Ge(001) as well as Si(001) because the energy difference between the ground state and the p(2$`\times `$2) structure is estimated to be a few meV/dimer Needels ; Inoue ; Yoshimoto . On the both surfaces, the neighboring two atoms form a buckled dimer, and their reconstructions are characterized by the ordering of the buckled dimers Chadi . Experimentally, fluctuation between c(4$`\times `$2) and p(2$`\times `$2) structures was induced from defects on Si(001)Shigekawa  or from Ag adsorbates on Ge(001)Naitoh  surfaces. Recently, the reconstructed structure of a clean Si(001) surface was studied in detail at low temperatures by STM Hata . One of the three structures, 2$`\times `$1, p(2$`\times `$2) and c(4$`\times `$2) is observed, depending both on the sample bias voltage and on the dopant species. This suggests that the band bending at the subsurface plays an important role to fix the structure as well as the electric field between the STM tip and the surface. It has been argued that the ground state energy of Si(001) surface is partly determined by the electrostatic energy among the electric dipole moments at the dimers Z . In contrast to the extensive studies on Si(001), the ground state of Ge(001) has not yet thoroughly investigated while c(4$`\times `$2) reconstruction was generally observed below 150 K Kevan ; Kubby . Up to now, the transition induced by the bias voltage change of STM was reported only from c(4$`\times `$2) to p(2$`\times `$2) on a vicinal Ge(001) surfaceRottger .

In the present letter, we show reversible change between c(4$`\times `$2) and p(2$`\times `$2) structures on a clean Ge(001) surface, which can be locally controlled by the sample bias voltage $`V_b`$ of STM with hysteresis below 80 K. The c(4$`\times `$2) structure is observed at negative $`V_b`$ $``$ \- 0.7 V while the p(2$`\times `$2) at positive $`V_b`$ $``$ 0.8 V. The both structure can be kept even under scanning the surface with $`|V_b|`$ $``$ 0.6 V at 80 K. These findings are important not only for the fundamental understanding of the stability of the reconstructed surface, but also for the application as a nanoscale memory.

All experiments were performed in an ultrahigh vacuum (UHV) system with a base pressure of below $`1\times 10^8`$ Pa consisting of a commercial variable temperature scanning tunneling microscope (OMICRON, LT-STM) and a surface preparation chamber Naitoh . Germanium specimens with the size of $`8\times 3\times 0.4`$ mm3 were cut from a Ge(001) wafer (Sb-dope, 0.35 $`\mathrm{\Omega }`$$``$cm at room temperature, 0.4 mm thick). The Ge(001) clean surfaces were obtained by several repetitions of Ar ion bombardment (1 keV, 2.5 $`\mu `$A$``$cm-2, 10 min) and annealing at 980 - 1000 K for 10 min by passing the DC current directly through the specimen in the UHV preparation chamber. Then, the specimen was transferred to the STM, and the surface was observed at 80 K or 10 K with a tungsten tip in a constant current mode.

Figure 1 shows a series of STM images obtained successively with three values of $`V_b`$ on the same area of Ge(001) surface at 80 K. In the case of $`V_b`$ = - 2.0 V the surface superstructure is c(4$`\times `$2) as shown in Fig. 1(a). We can keep this structure when $`V_b`$ is less than 0.6 V. Figure 1(b) shows the image with $`V_b`$ = - 0.2 V. When we further increase the bias to 1.2 V, the structure is p(2$`\times `$2) as shown in Fig. 1(c). Then, this structure can be maintained with $`V_b`$ down to - 0.6 V across the zero-bias as demonstrated in Fig. 1(d) for $`V_b`$ = - 0.2 V. The bias dependence of the surface reconstruction at 80 K is summarized in Fig. 1(e). We can fix either of the two structures at the same bias by changing the bias from $`V_b`$ 0.8 V or $`V_b`$ \- 0.7 V to $`|V_b|`$ 0.6 V. It is noted that the both structures can be kept without electric field due to the STM-tip; even after keeping the tip apex at 500 nm away from the area of the interest for several minutes, we observed again the same local structure at the area. The surface structure at 10 K depends on $`V_b`$ in the same way as that at 80 K. The observed bistability at $`|V_b|`$ 0.6 V indicates an energy barrier between the c(4$`\times `$2) and p(2$`\times `$2) structures. At the structural change, a flip process of the buckled dimers occurs in c(4$`\times `$2) or p(2$`\times `$2) domains. The energy required for such a process in c(4$`\times `$2) domain was theoretically estimated to be 416.7 meV based on a first-principles calculationKawai .

When we increase $`V_b`$ in several msec from - 2.5 V to + 1.2 V and continue to scan the surface, the surface gradually transforms from c(4$`\times `$2) to p(2$`\times `$2) through a random series of sudden flip of the dimers at a single Ge dimer-row as in Figs. 2(a-c). Here, we show selected successive STM images after the $`V_b`$ change. Arrows in the figures indicate the sudden changes of the dimer phase during the scanning of the surface from the bottom of the image to the top. When the buckling direction of a dimer changes, all dimers in the same dimer-row change simultaneously. That is, the change of the dimer phase continues to the outside of the scanning area along the dimer row until the dimer row is terminated by steps or any defects. The transformation of the structure was also observed from p(2$`\times `$2) to c(4$`\times `$2) when we decrease $`V_b`$, for example, from + 2.0 V to - 1.2 V. In contrast to the change from c(4$`\times `$2) to p(2$`\times `$2), however, this change is limited to a few nm away from the area of the scanning along the dimer-row direction. At the boundaries between the c(4$`\times `$2) and p(2$`\times `$2) areas, there form pairs of the adjacent dimers along the dimer row direction with the same buckling direction.

Figure 2(d) shows a scanning tunneling spectrum spatially averaged on a clean Ge(001) surface at 10 K. The observed peaks at -1.25 V and + 1.05 V are consistent with the previous result at RT Kubby . All the peaks appearing between -1.0 and 1.0 V are attributed to the dangling bond $`\pi `$ and $`\pi `$ states as in Si(001) surface Okada . The positions of the peaks in the tunneling spectra move toward $`V_b`$ = 0 at 80 K. This is mainly ascribed to the temperature dependence of the band bending nearby the surface.

The structures c(4$`\times `$2) and p(2$`\times `$2) are characterized by the phase difference between two adjacent dimer rows. If the dimers are all in phase, the domain structure is p(2$`\times `$2), and out of phase c(4$`\times `$2). Here we define “c(4$`\times `$2) length” on an area as the total length of the phase-mismatched adjacent dimers along the dimer-row direction. We estimated the change rate from c(4$`\times `$2) to p(2$`\times `$2) by measuring the c(4$`\times `$2) length in successive images of the same area as a function of time after a sudden $`V_b`$ increase at several values of the tunneling current, $`I_t`$. The results are shown in Fig. 3 for the $`V_b`$ change from - 2.0 V to 1.2 V. It took 6 min to obtain each image of 50 $`\times `$ 50 nm2. The c(4$`\times `$2) length after the first scanning decreases with increasing $`I_t`$. The buckling phases of a few dimer rows on this surface could not be changed by the scanning with $`V_b`$ = 1.2 V because the phase is firmly fixed by defects on the surface. Apart from these dimer rows, the characteristic times for the structural change are less than 6 min for $`I_t`$ = 1.6 nA and 12 min for 0.8 nA. We confirmed that the tip moves toward the surface just by 0.06 nm with increasing the tunneling current from 0.4 nA to 2 nA at $`V_b`$ = 1.2 V. This value is much shorter than the tip-surface distance in the tunneling condition. Consequently, the observed structure transformation is little influenced by the changes of the electric field and the atomic force accompanied by the tunneling current change. All these results indicate that the driving force for the change of the buckling phase is mainly attributed to the inelastic electron tunneling.

The phase of the dimer row suddenly changes during the scanning as shown in Fig. 2. The change rate of the structure decreases with decreasing the absolute value of $`V_b`$ or temperature. We can keep the both superstructures for more than 2 hours during the scanning on the surface, for example, with $`I_t`$ = 1 nA and $`V_b`$ = $`\pm `$ 0.6 V at 80 K. The change rate with $`I_t`$ = 1 nA and $`V_b`$ = 0.8 V decrease to about 1/3 by decreasing temperature from 80 K to 10 K. According to the tunneling spectrum shown in Fig. 2(d), electrons are effectively injected into or removed from the dangling bond $`\pi `$ or $`\pi `$ states by tunneling when $`|V_b|`$ 0.7 V.

Figures 4 demonstrate the successive formations and annihilations of local p(2$`\times `$2) areas by scanning the small area shown as white rectangles at 80 K. The images were observed with $`V_b`$ = - 0.2 V. After taking the image shown in Fig. 4(a), the area A in the figure was scanned with $`V_b`$ = 0.8 V to change from c(4$`\times `$2) to p(2$`\times `$2). We reduced $`V_b`$ as low as possible to avoid any influence to the dimer rows outside of the scanning area. The result is shown in Fig. 4(b). The change along the dimer row continued to the outside of the scanning area as described above. Then we scanned the area B in Fig. 4(b) to form another p(2$`\times `$2) area. The result is shown in Fig. 4(c). The reverse change from p(2$`\times `$2) to c(4$`\times `$2) was done by scanning the area C in Fig. 4(c) with $`V_b`$= - 0.7 V. In this case, we had to scan all of the desired p(2$`\times `$2) area because the change to c(4$`\times `$2) was limited nearby the area of the scanning. The annihilation of p(2$`\times `$2) is shown in Fig. 4(d).

The ground state of Ge(001) surface has been believed to be c(4$`\times `$2) structure. Both LEED observationKevan  and STM studyKubby  reported c(4$`\times `$2) structure at low temperature consistently with the theoretical studies Needels ; Inoue ; Yoshimoto . The surface superstructure of Ge(001) is determined by minimizing two competing energies. One is the strain energy at subsurface latticeYoshimoto , and the other is mutual electric-dipole energy among the Ge dimersZ . On this surface, the second-layer atoms under the upper-atom of the dimer are pulled closer while those under the lower-atom are pushed apart to minimize the energy of lattice distortion due to the asymmetric dimer. In the c(4$`\times `$2) structure, the upper-atom of a dimer is adjacent to the upper-atom in its adjoining dimer row while it is adjacent to the lower-atom in the p(2$`\times `$2) structure. Thus, the strain energy of the c(4$`\times `$2) structure at subsurface is lower than that of the p(2$`\times `$2). On the other hand, the dipole interaction energy in the p(2$`\times `$2) structure is smaller than that in the c(4$`\times `$2) because the dipoles are oriented to the same direction in p(2$`\times `$2). We note that the upper-atom of the dimer is negatively charged and the lower positively. The c(4$`\times `$2) structure can be the ground state because the energy gain of the lattice-strain dominates the loss of the dipolar energy in this structure. The difference of the energy between p(2$`\times `$2) and c(4$`\times `$2) is as small as 3 meV/dimerYoshimoto .

When the surface is positively biased under the STM tip, the negatively charged upper-atoms of the dimers are electro-statically pushed toward the subsurface by the tip, and the lower-atoms of the dimers are pulled toward the surface. Consequently, the dimer buckling becomes smaller. This makes the difference of the subsurface strain energy between p(2$`\times `$2) and c(4$`\times `$2) structures smaller than that of the dipole energy among buckled dimers. Thus for the positively biased surface, the p(2$`\times `$2) structure is favored, and can be the ground state. On the other hand, under the negative $`V_b`$, the buckling of the dimers are increased by the electric field, and the c(4$`\times `$2) structure is more favored. This mechanism is qualitatively consistent with the stable structures observed at positive and negative biases by STM.

At present, however, we can quantitatively infer neither the variation of the stability difference between c(4$`\times `$2) and p(2$`\times `$2) structures nor the flipping probability of a dimer by inelastic electron tunneling under the electric field due to the tip-sample bias. The first-principles calculation of these parameters seems very difficult because of the rather long-scale band bending of Ge subsurface possibly requires extremely large scale calculations, for example. To clarify the transition mechanism in detail, we should tackle the above two, which are largely dependent on the local electric field and the local current at the dimer.

In case of Si(001) surface, the ground state has been considered to be c(4$`\times `$2) as in Ge(001). Recently, however, both 2$`\times `$1 and p(2$`\times `$2) structures were observed by STM below 50 K. Hata et al.Hata  claimed that the observed 2$`\times `$1 structure at 10 K is attributed to the STM-tip induced flip-flop of the Si dimers. In this case, p(2$`\times `$2) structure appears only in n-type Si, and c(4$`\times `$2) only in p-type Si. Moreover, the latter was observed only in a very small range of $`V_b`$. These are quite in contrast to the present results at Ge(001) surface. We have never observed 2$`\times `$1 structure or symmetric dimers on Ge(001) surface below 80 K except at the dimer rows adjacent to defects on the surface. The both c(4$`\times `$2) and p(2$`\times `$2) structures are observed in a wide $`V_b`$ range, and are stable during the STM observations with $`V_b`$ between - 0.6 V and 0.6 V. The difference from Si(001) surface is ascribed to the large energy barrier between c(4$`\times `$2) and p(2$`\times `$2) structures in Ge(001) Inoue ; Yoshimoto .

In conclusion, we have demonstrated the local modification of the Ge(001) reconstruction between c(4$`\times `$2) and p(2$`\times `$2) reversibly by $`V_b`$ change using STM below 80 K. The p(2$`\times `$2) structure can be the ground state under an electric field, and inelastic tunneling is attributed to the driving force to the local transformation of the structure. The two different structures coexist with $`|V_b|`$ 0.6 V, which is applicable to a nanoscale memory.