text:
# Geometrical Dependence of High-Bias Current in Multiwalled Carbon Nanotubes

## Abstract

We have studied the high-bias transport properties of the different shells that constitute a multiwalled carbon nanotube. The current is shown to be reduced as the shell diameter is decreased or the length is increased. We assign this geometrical dependence to the competition between electron-phonon scattering process and Zener tunneling.

PACS numbers: 73.63.Fg, 73.50.Fq, 72.10.Di

Building electronic devices at the molecular scale has motivated intense research for the last years. Carbon nanotubes have been used as building blocks for devices such as nanotube-nanotube junctions  or field-effect transistors  that have been integrated in logic circuits . Recently nanotubes have also attracted widespread attention as future interconnects, since their current-carrying capacities are several orders of magnitude larger than in present-day interconnects. The current in metal single-wall nanotubes (SWNT) saturates at $`2025`$ $`\mu `$A, which corresponds to a current density exceeding $`10^9`$ A/cm2 .

An important question that remains to be addressed is how the current saturation grows with the section of an interconnect made of nanotubes. The interconnect may be made larger by selecting SWNTs assembled in a bundle. However, the intertube transmission is weak and leads to a bad current distribution over the whole bundle when the current is injected in one tube. Another approach is to choose multiwalled nanotubes (MWNTs) with many concentric shells. The intertube transmission within a single MWNT is larger. The interface of two neighbor shells is maximum and equal to the entire shell surface. Recently, it has been shown that most of the shells contribute to the current at large bias  and that the current of each shell saturates like in SWNTs . Remarkably, the saturation value was reported to be close to what is obtained for SWNTs and to be independent on the shell diameter. This is very surprising since more subbands are expected to carry current in large shells. MWNT shells can be much larger than SWNTs, reducing substantially the energy separation between subbands.

In this Letter, we present studies of the high-bias transport in MWNTs with different electrode separations $`L`$. We find that the current saturation per shell varies between $`10`$ and $`60\mu `$A. For the shortest separation ($`200`$ nm), the current is shown to increase as a function of diameter. For longer separations ($`1\mu `$m), the diameter dependence is much weaker and more difficult to observe. We propose that this geometrical dependence is due to a weak variation of the number of current carrying subbands within each shell that results from the competition in the electron transmission between the electron-phonon scattering process and Zener tunneling. The Landauer formula is applied that incorporates these two processes and the numerical calculations reproduce well the experimental results. Interestingly, the model describes the experimental observation that the high-bias current does not depend on the metal or the semiconducting character of the shell at large diameters.

The MWNTs were synthesized by arc-discharge evaporation and carefully purified . The MWNTs are sonicated in dichloroethane and dispersed onto an oxidized Si substrate. Atomic force microscopy is used to image, select and locate nanotubes that are then contacted with Cr/Au electrodes by electronic lithography. An example of a device is shown in Fig 1(a). Typical two-point resistances at low-bias vary in the range $`5`$ to $`30`$ k$`\mathrm{\Omega }`$.

The electronic properties of different shells in a MWNT are probed using the electrical-breakdown method developed by Collins $`et`$ $`al.`$ . The shells are selectively removed stepwise by injecting a large current in the MWNT. The bias voltage is increased until a sharp step in the current is detected that corresponds to the failure of one shell. Then, the bias voltage is quickly put down to 0 in $`100`$ ms. The removed shell has been shown to be the outermost one which is in direct contact with air and therefore most subject to a rapid oxidation initiated by the current . We remove shells until no more current is detected. The final diameter has then decreased to typically $`24`$ nm as seen on the AFM picture in Fig. 1(b). The image shows that the peeling is uniform along the tube, except near the electrodes. Note that the failure is apparently different for the last shells since they are not evaporated. Failure may be local and situated at a small gap, too short to be observed with AFM.

We have designed an experiment where a series of electrodes is attached to a single long MWNT with various separations from $`200`$ to $`1200`$ nm. Fig. 2(a-c) show the result of the electrical-breakdown method applied to the different sections of this MWNT. Let us first consider the current-voltage characteristics ($`IV`$) of the pristine sample given by the upper $`IV`$ curve. It appears that the high-bias current is significantly larger in the shortest section. This result can be enlightened by analyzing the current carried by the different shells of the MWNT that are obtained by comparing the successive $`IV`$s. Inspection of Fig. 2(a-c) shows that the current increase in the 200 nm section is the combined effect of $`(i)`$ more jumps and $`(ii)`$ larger jumps.

We discuss first the jump number. This 15 nm diameter MWNT is expected to consist of $`19`$ shells. The jump number is however smaller and $`L`$ dependent. It is reduced from 17 to 13 when $`L`$ passes from 200 to 1200 nm. A smaller number of jumps might be due to a smaller number of shells and to a thinner diameter. However, this is unlikely since AFM imaging indicates a uniform diameter over the entire length. This observation suggests rather that some shells do not participate to the current, especially in the longer sections. Such an interpretation is supported by the following observation in Fig. 2(d), where the maximum current of each $`IV`$ is plotted as a function of the diameter $`d`$. The diameter has not been measured by AFM at each step but estimated assuming that the diameter decreases each time by twice the interlayer distance $`3.4`$ Å. The extrapolated curves should pass by $`25`$ $`\mu `$A at $`1.4`$ nm, the saturation of a SWNT. This is not the case. Rather, the extrapolations have negative current values at $`1.4`$ nm, which indicates that $`d`$ is overestimated and suggests that some shells do not participate to the current.

As mentioned before, the current jumps for the shortest length are larger in amplitude. The jump averages are 28.7, 19.1 and 17.4 $`\mu `$A for respectively the 200, 750 and 1200 nm lengths. Remarkably, the curves in Fig. 2d deviate from a linear behavior. The deviation is best seen for the $`200`$ nm length. It shows that large-diameter shells contribute more to current. This finding has been observed in more than $`20`$ MWNTs. Other examples are shown in Fig. 3(a) and (b), where the current jump $`\mathrm{\Delta }I`$ is plotted as a function of $`d`$. Both curves show that $`\mathrm{\Delta }I`$ is larger for larger $`d`$. However, the slope is much weaker for the $`650`$ nm length and is hardly observed due the measurement noise.

We note that an accurate determination of the current jump is difficult. Inset of Fig. 3(a) shows that the bias at which the shell fails increases irregularly and that the $`IV`$s are not smooth. We define the current-jump as the current difference at a bias $`V_0`$ given by the maximum voltage recorded for the two successive $`IV`$s. Another approach could have been to take the maximum of the difference or the difference of the maximums. However, the observation of a current increase with diameter is not altered by the choice of the estimation procedure. Importantly, the current jump estimation is precise enough to obtain close slopes of $`\mathrm{\Delta }I(d)`$ for different MWNTs with same length as seen in Fig. 3(b).

We now discuss the possible origins of these geometrical dependences in the shell current. We review first the electron-phonon scattering process that causes the current to saturate at $`2025\mu `$A in metal SWNTs . An electron that is accelarated by an electric field $``$ will gain energy until an optical phonon of energy $`\mathrm{}\mathrm{\Omega }`$ is emitted ($`\mathrm{}\mathrm{\Omega }`$ 0.16eV). At this point, the coupling to the phonon is so strong that the electron is backscattered. Using a Landauer-type argument the current saturation is shown to be proportional to the phonon energy, $`(4e/h)\mathrm{}\mathrm{\Omega }25\mu `$A. The mean free path $`l_{ph}`$ is equal to the distance an electron must travel to be accelerated to $`\mathrm{}\mathrm{\Omega }`$, $`l_{ph}=\mathrm{}\mathrm{\Omega }/e`$ . The total transmission trough the whole tube is obtained by dividing $`l_{ph}`$ by $`L`$

$$T_{ph}=\mathrm{}\mathrm{\Omega }/eV$$
(1)
where V is the voltage applied between the electrodes. Eq. (2) is valid for large $`V`$.

At first sight, the current jump variations might be explained from the energy variation of phonons as a function of diameter. However, $`ab`$ $`initio`$ calculations have shown that optical phonon energy changes barely with tube diameter . An alternative mechanism is thus needed to account for the current-jump variation.

We propose that the current sensitivity to diameter results from the contribution of Zener tunneling . Fig 4(a) shows a scheme that represents the potential variation in space of the band structure. The potential drop should be in principle determined by the self-consistent solution of the Poisson equation and the nonequilibrium electron density. We assume here that the voltage drop is linear, a good approximation when multiple electron-phonon scattering processes occur .

An electron that enters a non-crossing subband from the left electrode reaches at some point the top of the valence band and tunnels to the conduction band. The Zener-tunneling transmission  is given by

$$T_Z=\mathrm{exp}(\frac{4\sqrt{2m^{}}LE^{3/2}}{3e\mathrm{}V})$$
(2)
where we have assumed a linear voltage drop. $`m^{}`$ is the effective mass and $`E`$ the gap that goes like $`1/d`$.

Now, we compare the two transmissions as a function of diameter. The phonon transmission is independent on $`d`$ and $`T_{ph}=0.045`$ at $`3.5`$ V. For a large semiconducting shell of $`15`$ nm diameter and $`L=200`$ nm, we obtain $`T_Z0.7`$ for tunneling between the upper valence band and the lower conduction band. The effective mass has been estimated using the thigh-binding model for a zig-zag tube. The comparison between $`T_{ph}`$ and $`T_Z`$ shows that the current is limited by phonons. Zener reflexion is here unimportant and no difference in current is expected between metal and semiconducting shells in agreement with experiment. The situation changes for small-diameter shells. For $`d=3`$ nm we estimate $`T_Z10^4`$ so that the current is now limited by Zener tunneling .

To put the above analyze on a more quantitative basis, we use the Landauer formula to estimate the high-bias current of a shell. Since the considered transmissions are weak, the total transmission of the subband $`i`$ is given by $`T_i^1=T_Z^1(E_i)+T_{ph}^1`$. $`E_i`$ is the energy separation between the $`i`$ valence band to the $`i`$ conduction band, where $`i`$ is counted from the neutrality point. Note that alternative processes such as tunneling from the $`i`$ valence band to the first conduction band changes barely the results presented below. The current is given by $`I=(2e/h)_iT_idE`$. The integral is calculated between $`0`$ and $`eV`$. The use of the Landauer formula in presence of inelastic scattering in the non-equilibrium regime is rather crude since it does not account for the exclusion principle. But it allows to give some physical insight. A model based on Boltzmann equation as in  that incorporates Zener tunneling is probably needed for more accuracy.

Fig 4(b) and (c) show the calculated current as a function of diameter of metal and semiconducting shells. The current increases with diameter and decreases with length, in agreement with experiment. Indeed, larger diameters decrease the subband gap and enable better Zener transmission. Similarly, the tunneling transmission becomes weaker for longer nanotubes as the Zener barrier length increases.

Fig. 4(d) shows the calculated current as a function of the diameter $`D`$ of MWNTs. We assume that MWNTs consist of shells that carry current in parallel. The current through the MWNT is the sum of the current of a series of shells with diameter ranging from $`2`$ nm to $`D`$. For each shell an average is taken on the current of a semiconducting and a metal shell with the statistical weights $`\frac{2}{3}\frac{1}{3}`$. The resemblance to the experiment is remarkable. $`No`$ $`fitting`$ parameter has been used.

Here, we discuss the number $`N`$ of subbands within a shell that carry significantly current. The calculated curve $`\mathrm{\Delta }I(d)`$ in Fig. 4c has 3 plateaus indicated by arrows. They correspond to situations where 0, 1 and respectively 2 subbands carry current. These subbands are close to the neutrality point and the Zener barrier is negligible compared to phonon scattering, $`T_Z(E)>T_{ph}`$. The other subbands that have no contribution to the current correspond to deeper valence bands with large $`E`$. Zener reflection has thus a dramatic role in the high bias transport in MWNTs. About $`eV/E_{sub}60`$ subbands in a $`14`$ nm shell are counted for $`eV=3.5`$ V and have to be taken into account for the current estimation, but very few contribute effectively.

Fig. 4(c) shows the current of semiconducting shells as a function of diameter. It appears that thin diameter shells do not carry current. This happens when $`T_Z(E)<T_{ph}`$. Moreover, the diameter at which semiconducting shells start to conduct becomes larger for longer tubes. This may explain the discussed experimental observation that the number of jumps is reduced for longer section. Such an interpretation is further supported as follows. Inset of Fig. 3(a) shows that the $`IV`$ indicated by an arrow fails at a bias larger than the others. Such $`IV`$s have been observed several times but mostly for the last $`IV`$. At the failure the current is typically $`50100`$ $`\mu `$A and corresponds to several shells. We propose that the outermost shell is semiconducting in this situation and do not carry current. The current has to pass through inner shells that are protected from oxygen and that can sustain higher voltages. At the failure, the breakdown is mostly so violent that all the shells fail together.

We finish with a discussion on the performances of large-diameter MWNTs as interconnects. We have shown that the number of current-carrying subbands within a shell depends weakly on the diameter and the length and stays typically close to one. Thus, the current saturation of MWNTs is expected to vary roughly linearly with the shell number and thus the diameter. This is observed in Fig. 2(e), where the maximum current as a function of diameter is plotted for 10 samples. The functional form is not quadratic as it would be expected for a plain wire. This makes current-carrying capacities of MWNTs less efficient as interconnects when they become larger. However, the current carrying capacity of large MWNTs still surpasses Cu wires ($`10^6`$ A/cm2). For example, the current of a $`100`$ nm diameter MWNT is estimated to be $`10100`$ times larger.

The authors wish to thank N. Regnault, P. Morfin and C. Delalande. The research has been supported by the CNRS, DGA, ACN programs, Paris $`6`$ and Paris $`7`$.

corresponding author: bachtold@lpmc.ens.fr


format: text

text:
# Ratio of color and effective temperatures in X-ray burst spectra

## 1 Introduction

Discovery of the first X-ray burster was reported by Grindlay et al. (1976) and by Belian et al. (1976). We believe that these objects are low mass binary stars which contain the accreting neutron star. Since we do not observe X-ray pulsations, the magnetic field of the neutron star is lower than $`10^{10}`$ Gs. Weak magnetic field of such intensity allows for the existence of accretion disc. Basic features of X-ray burst are: fast rise $``$ 1 s and decay time $``$ 3 – 100 s, luminosity near the peak about $`10^{38}`$ erg/s and energy per burst $`10^{39}`$ ergs (cf. review article by Joss & Rappaport 1984). The idea that energy source of a burst is the thermonuclear flash was proposed by Woosley and Taam (1976) and Maraschi & Cavaliere (1977). This model quite well explains many global features of X-ray bursters.

Initially neutron star spectra were assumed to be blackbody spectra, but this assumption is not valid when the atmosphere was dominated by scattering on free electrons. In such case spectra of outgoing photons are shifted towards higher energy and the maximum is shifted in the same direction both for the coherent Thomson, and noncoherent Compton scattering. In electron scattering stellar atmosphere, photons are created by thermal processes below the photosphere, in so called thermalization layer. In this layer local temperature T is much higher than the effective temperature $`T_{eff}`$, therefore maximum of the Planck function $`B_\nu (T)`$ occurs in higher energy than the maximum of $`B_\nu (T_{eff})`$. Consequently, the blackbody temperature $`T_{BB}`$ as obtained from fitting of continuum spectra, can be much larger than the true $`T_{eff}`$. Temperature $`T_{BB}`$ is usually referred to as the color temperature $`T_c`$.

In principle, the ratio $`T_c/T_{eff}`$ can reach arbitrary large values in Thomson scattering stellar atmospheres, provided that true absorption in that atmospheres decreases to zero (Madej 1974 ). Compton scattering in stellar atmospheres reduces $`T_c/T_{eff}`$ by unknown factor, but still this ratio remains larger than unity. Model atmosphere computations with account of Compton scattering, and exact $`T_c/T_{eff}`$ determination, are very difficult.

Ebisuzaki, Hanawa & Sugimoto (1984) estimated $`T_c/T_{eff}1.4`$ by use the simple opacity formulae. Model atmospheres computations by London, Taam & Howard (1986) showed, that this ratio was in the range from 1.2 to 1.7 and depended on few parameters: $`T_{eff}`$, surface gravity log g and iron abundance. London, Taam & Howard (1986)  concluded that the solar iron abundance does not influence significantly the ratio $`T_c/T_{eff}`$. On the other hand, Ebisuzaki & Nomoto (1986) claimed that absorption by iron decrease this ratio.

## 2 Model equations

Our model equations describe transfer of radiation subject to constraints of hydrostatic and radiative equilibrium. We assume the equation of state of ideal gas in local thermodynamical equilibrium (LTE). The equation of transfer includes f-f absorption from all ions, b-f absorption from iron ions, and Compton scattering terms. The equation of transfer is as follows:

$`\mu {\displaystyle \frac{I_\nu }{\tau _\nu }}`$ $`=`$ $`I_\nu ϵ_\nu B_\nu (1ϵ_\nu )J_\nu `$
$`+`$ $`(1ϵ_\nu )J_\nu {\displaystyle _0^{\mathrm{}}}\mathrm{\Phi }_1(\nu ,\nu ^{})𝑑\nu ^{}`$
$``$ $`(1ϵ_\nu ){\displaystyle _0^{\mathrm{}}}J_\nu ^{}\mathrm{\Phi }_2(\nu ,\nu ^{})𝑑\nu ^{}`$
where: $`B_\nu `$ is the Planck function, $`J_\nu `$ the mean intensity, $`ϵ_\nu `$ is the ratio of true absorption to coefficient of total opacity, $`\tau _\nu `$ the monochromatic optical depth, and $`\mu `$ denotes cosine of zenithal angle.
Details of the above equation of transfer, and Compton scattering redistribution functions $`\mathrm{\Phi }_1(\nu ,\nu ^{})`$ and $`\mathrm{\Phi }_2(\nu ,\nu ^{})`$, and linearization procedure, are described in earlier papers (see Madej 1991 ; Madej & Różańska 2000a, b; Madej, Joss & Różańska 2003).

## 3 Computations

We calculated our models with the following assumptions: equation of state in LTE, hydrostatic and radiative equilibrium, plane-parallel geometry, non rotating star. Our equations allow for large photon-electron energy exchange at the time of Compton scattering. Details of the computer code and linearization scheme were extensively described in earlier papers quoted in previous section.

We calculated grid of models each of them computed on the 144 standard optical depth points (12 points per the decade) from $`10^4`$ to $`10^8`$ and 901 wavelength points. We computed few models of hot neutron star atmospheres with effective temperatures $`T_{eff}`$=1, 2, 3$`\times 10^7`$ K, and surface gravities $`g`$ from the critical gravity up to $`10^{15}`$ in cgs units. Models assume H and He in solar number abundance $`N_{He}/N_H=0.11`$ and iron with number abundance $`N_{Fe}/N_H=3.7\times 10^5`$ (approximately solar value) and 100 times greater.

## 4 Results

We claim that our computations can measure both the ratio $`T_c/T_{eff}`$ and its slight changes influenced by small differences of $`T_{eff}`$, $`\mathrm{log}g`$ and iron abundance. We obtain models of burster atmospheres where the constancy of the bolometric flux is almost everywhere better than 1%, and only in few optical depth points approches 2%.

Tables 1 and 2 show values of $`T_c/T_{eff}`$ computed over wide range of parameters $`T_{eff}`$, $`\mathrm{log}g`$ and iron abundance. We can see that for models with pure H-He atmospheres, these ratios are higher than for atmospheres with nonzero iron abundance, of course for the same $`T_{eff}`$ and $`\mathrm{log}g`$. For higher metal abundance these ratios decrease. This is because increase of metal abundance diminish the importance of Compton scattering, and $`T_c/T_{eff}`$ ratio is than determined by the redistribution of outgoing flux by nongray metal absorption.

What is interesting, for the lowest $`T_{eff}=10^7`$K these ratios exhibit local mimimum in some intermediate $`\mathrm{log}g`$. Tendencies, which we observe are due to iron opacity and the Compton scattering. For the lowest surface gravity Compton scattering dominates over thermal absorption, therefore shift of the peak spectra is the highest. In the presence of iron, thermal absorption is greater than in case of H-He atmospheres. Consequently, the $`T_c/T_{eff}`$ is lower for atmospheres containing iron.

In Figure 1 we show the outgoing photon spectra for various parameters. We compare blackbody spectra with the spectra of H-He atmosphere and with the atmosphere containing iron in Figure 1a. We can see that the blackbody spectra (solid line) have the peak in lower energy and the BB spectrum is softer than the spectrum of the model X-ray burster. In lower photon energies we obtain excess of the outgoing flux as compared to the blackbody spectra both for H-He models and atmospheres containing some amount of iron. Figures 1c and 1d present spectra of atmospheres with iron abundance 100 times greater than the solar abundance, and with the solar iron abundance, respectively. For the latter iron abundance we note, that the iron ionisation edge appear in emission for two lowest $`T_{eff}`$. Profiles of iron absorption lines (fundamental series of H-like, and He-like ions) can be very complicated (see Fig. 1b). In some models we obtain inverse emission in iron line cores, which strictly correspond to temperature inversion reproduced in our model atmospheres with Compton scattering.

We stress here, that we did not obtain values of $`T_c/T_{eff}`$ larger than 2.0 in any of our model atmospheres of X-ray burst sources.

## ACKNOWLEDGEMENTS

We acknowledge support of the grant No. 2 P03D 021 22 from the Polish Committee for Scientific Research.


text:
# Chapter 1 Whipple Telescope Observations of Potential TeV Gamma-Ray Sources Found by the Tibet Air Shower Array

## Chapter 1 Whipple Telescope Observations of Potential TeV Gamma-Ray Sources Found by the Tibet Air Shower Array

### Abstract

An all-sky survey performed with the Tibet Air Shower Array (Tibet AS) has found a number of potential point sources of TeV gamma rays. If they are steady sources, the implied Tibet AS fluxes should be visible with strong significance to the Whipple 10-m gamma-ray telescope (E $`>`$ 400 GeV) with only a short (5 hour) exposure. We have observed four candidate directions from the Tibet-II HD dataset for $``$5 hours each with the Whipple telescope. In addition, we observed a new candidate direction from the Tibet-III Phase 1 dataset for 7.5 hours. We have found no corresponding excesses at the flux levels implied, and we have set upper limits for each candidate.

### 1. Introduction

The Tibet AS is an array of particle detectors located in Tibet at an elevation of 4300 m and operated by the TibetAS$`\gamma `$ collaboration. The Tibet-II HD array operated from 1996 to 1999 with an energy threshold of 3 TeV and an angular resolution of $`0.9^{}`$. The Tibet-III array has operated since 1999 with the same threshold and resolution . The ability of these arrays to detect gamma-ray sources has been demonstrated through detection of the Crab Nebula  and Mrk 501 .
A wide angle survey conducted with the Tibet-II HD array found 19 directions with excesses greater than 4$`\sigma `$ over the average background. The TibetAS$`\gamma `$ collaboration noted that these may be explained as statistical fluctuations, but one direction corresponded to the Crab Nebula . We searched radio, optical, x-ray, and gamma-ray catalogues for corresponding objects within 1 of each of the Tibet directions, and we selected four promising candidate directions for observation with the Whipple telescope during the 2001-2002 season. After receiving an update from the Tibet-III all-sky survey, we selected one more candidate for observation during the 2002-2003 season. Table 1 gives a summary of the targets and observations. The candidates were selected as follows:

* Tibet1 had a high significance and showed steady increase through Tibet-II HD data.
* Tibet9 is $`0.3^{}`$ from a Seyfert 1 galaxy (RGB J1337+243).
* Tibet14 is $`0.7^{}`$ from an EGRET unidentified (3EG J2021+3716).
* Tibet16 had a high significance and is in the Cygnus star field.
* Tibet0554 showed steady increase through Tibet-II HD and Tibet-III and was second in significance to the Crab.

### 2. Observations and Analysis

Observations were made with the Whipple 10 m gamma-ray telescope with the 490 pixel camera . Only the inner 379 pixels ($`2.4^{}`$) were used in the analysis, and because of the large uncertainty in the Tibet AS source coordinates, 2-dimensional (2-D) analysis was required. In the analysis, potential gamma rays were selected by applying standard supercuts shape cuts. Then the distance and alpha cuts ($`\alpha <10^{}`$) were applied across a grid of points, with alpha and distance calculated with respect to each point. The background was estimated using either OFF-source data, or an average background comprised of many OFF runs. We applied this technique to data with the Crab off-axis and found it to be effective, as seen in Figure 1.

### 3. Results

The 2-D analysis of the 2001-2003 Whipple data has failed to detect a source in any of the target regions. A 2-D significance plot for Tibet0554 is shown in Figure 2 as an example. Assuming that each candidate is point-like and has a spectral shape similar to the Crab, we are able to set upper limits on the TeV gamma-ray flux from each candidate region. These flux limits are summarized in Table 2.

It should be noted that this analysis cannot rule out the possibility that the Tibet AS had observed an extended source ($`size>0.2^{}`$) or a line source with $`Energy>`$ 3 TeV. However, assuming the candidates are point-like with power-law energy spectra, we conclude that the five Tibet AS excess regions we studied were either statistical fluctuations in the Tibet AS data, or due to the non-contemporaneous observation periods, the emissions were episodic.

### 4. Acknowledgement

We gratefully acknowledge support from the University of Utah and the National Science Foundation under NSF Grant #0079704. This research is also supported by grants from the U.S. Department of Energy, by Enterprise Ireland and by PPARC in the UK. We also acknowledge the technical assistance of E. Roache and J. Melnick.

### 5. References

1. Amenomori et al. 2001, Proc. 27th Int. Cosmic Ray Conf (Hamburg), 2544

2. Amenomori et al. 2002, ApJ, 580, 887

3. Amenomori et al. 1999, ApJL 525, L93

4. Amenomori et al. 1999, ApJ, 532, 302

5. Finley, J.P. et al. 2001, Proc. 27th Int. Cosmic Ray Conf (Hamburg), 2827


text:
# High Performance Commodity Networking in a 512-CPU Teraflop Beowulf Cluster for Computational Astrophysics

## 1 Introduction

The Olympic motto ”Citius, altius, fortius” (”Swifter, higher, stronger”) succinctly describes the direction of 21st century parallel supercomputing - swifter processors, higher resolution and stronger fault-tolerant parallel algorithms. In computational astrophysics, the need for all of these qualities is perhaps greater than most supercomputing applications. While the bulk of the physics of the formation of the planets, stars, galaxies and the large-scale structure in the universe are now largely understood with the initial conditions well posed in some cases, the challenge of computing the formation of objects and structures in the universe is difficult. The standard methods of computational fluid dynamics (including magnetohydrodynamics) and gravitational N-body simulation are stressed by the large dynamic range in density, pressure, and temperature that exist in nature. Ever-finer computational meshes and greater numbers particles in N-body simulations are needed to capture the physics of the formation of things in the universe correctly.

In recent years, the relative low cost of commodity servers running the Linux operating system and more importantly switches and network interface cards has led to the growth in performance and competitiveness of Beowulf clusters versus vector machines and symmetric multi-processors. The recent list of the Top 500 supercomputers  last year reveals that 55 out of the top 100 are clusters. Most high performance Beowulfs now use low-latency, proprietary networking with costs several times that of current commodity gigabit networking. Large-scale 200+ port gigabit switches are a slightly cheaper option but still eat significantly into the cost of a cluster. Networking costs alone can dominate cluster costs.

The spirit of building a Beowulf cluster is to get the highest performance on a limited budget so the incorporation of expensive networking solutions defeats this fundamental motivation. We describe here a cheaper networking configuration based on off-the-shelf 24-port commodity gigabit switches in a 256-node/512-CPU cluster. We achieve a total cross-sectional bandwidth of 128 Gbit between 256 nodes. The innovation to achieving this performance is a networking topology involving both a stack of commodity 24-port gigabit SMC switches and the configuration of each Linux server as a network router through a second gigabit network port. Our economical networking strategy competes directly with both large-scale 200+ port switches and proprietary low-latency networking at a fraction of the cost. The high throughput provided by Linux routing capability makes this performance possible. We ran the HPL LINPACK standard benchmark and achieved 1.202 Tflops placing it in the top 50 computers according to the Top 500 list of November 2002.

We also describe the performance and results of 4 parallel applications for astrophysical problems that are currently running on the cluster spanning problems in large-scale structure formation, galactic dynamics, magnetohydrodynamic flows on blackholes and planet formation.

## 2 Cluster Hardware Specifications

In this section, we briefly describe the design and hardware specifications of the cluster. The hardware was assembled and installed by Mynix Technologies, Montreal, PQ at the University of Toronto . The installation of the Linux operation system was done using the OSCAR  cluster package by our group.

### 2.1 Compute Nodes

The CITA Beowulf Cluster dubbed McKenzie is comprised of 268 dual Xeon 2.4 Ghz rack-mounted format Linux servers (536 processors in total) distributed over 7 racks. All the nodes are 1U format with the exception of two 2U format head nodes. The hardware specifications for each server are given in Table 1. From the total of 268 nodes, 256 nodes are dedicated to parallel, message-passing supercomputing, 8 nodes for code development and smaller scale applications, 2 spare nodes also running to act as hot replacements in case of hardware failures on the compute nodes and 2 head nodes. The main head node “bob” contains the home disk space as well as the root space images and Linux kernels for the compute slave nodes. All slave nodes boot through the network from bob using PXElinux allowing easy kernel upgrades and maintenance for the cluster nodes. A secondary master node “doug” mirrors bob’s root and home directories on a daily basis and acts a backup in case of a bob crash. This allows quick recovery and continued operation if bob fails.

### 2.2 Switches and Cables

The cluster is networked using 19 SMC Tiger 24-port configurable gigabit switches. Seventeen switches are used for the main 256-node cluster network while the development cluster is connected to its own switch and the 19th switch is reserved as a spare and connected to the running spare nodes. We defer the discussion of the network topology to the section below.

### 2.3 Rack Configuration

The nodes and switches are mounted on 7 44U high racks. The main 256 compute nodes are mounted on the first 3 and last 3 racks. The central rack is reserved for the head nodes, development nodes and spares as well as the stack of switches. The switches are mounted on the back of the racks to simplify the cabling.

### 2.4 Cost

The cost for the compute nodes plus storage and incidentals was USD$513K. The cost for the gigabit switches plus cabling was USD$48K representing only 9% of the total cluster cost of $561K.

## 3 Networking

When considering options for networking the cluster, we decided against proprietary low-latency networking and large-scale 200+ port switches after realizing that they would consume a significant fraction of our modest budget. The prices of 24-port gigabit switches this past year had dropped significantly so we considered networking schemes that could take advantage of this inexpensive hardware.

Fat-tree networks that link switches together through a hierarchy are a straightforward configuration that provide connectivity but quite low cross-sectional bandwith. Each node comes with a second gigabit port, however, and could be exploited in some way using point to point connections with other machines. With this in mind we came up with the following scheme. Switches (or pairs of trunked switches) can be thought of as separate tightly coupled networks of up to a few dozen nodes. These nodes can be assigned to vertices in some more complex network topology. How can one connect these vertices to form a high-bandwidth global network linking all the compute nodes? After loading a switch with compute nodes there are few ports left to link to other vertices. However, each compute node provides a single connection through its second port that can in principle be connected to a compute node on an adjacent network vertex. If each linux node has routing enabled, network traffic can be relayed through the second port to other nodes. This was the strategy we ultimately employed in setting up a high-bandwidth global network.

### 3.1 The Fat-Tree Maintenance Network

Before setting up our high performance global network, we first set up a simple fat-tree that we treat as a robust but lower-performance maintenance network. This network was used to install the cluster originally but also acts as a bootstrap for configuring a high-performance network. In our final configuration, we trunk pairs of switches together making network vertices containing 48 ports with a total of 8 vertices using 16 switches. We connect 32 compute nodes to each vertex. For 6 vertices, we run a 4-port trunk to a 17th master switch filling it to capacity. The fat-tree networking is completed by connecting the last two vertices with 4-port trunks to 2 vertices trunked to the master switch. The head nodes bob and doug, plus the development and spare nodes are all plugged into this network as well allowing direct communication between all available nodes.

This fat-tree network configuration is robust and provides connectivity to all the 256 nodes and is adequate for maintenance, installation and embarassingly parallel applications. However, the cross-sectional bandwidth is not at all ideal for heavy duty MPI applications. We now describe how we use this fat-tree network to bootstrap to our high-bandwidth network which runs predominantly through point to point connections between the second network ports of all machines.

### 3.2 The Cross-Diagonally Connected Cube Network

The network vertices described above can be thought of as independent networks that we can connect to each other according to some topology. The use of a master switch to build a fat tree is one such topology but it has quite low cross-sectional bandwidth. How can we use the second port on each compute node to connect these vertices to increase the network bandwidth and minimize network hops to reduce latency? There are many possibilities but we finally settled on a cubic network topology (Figure 1).

In our chosen topology, the 8 network vertices are arranged on the corners of a cube. Communication between vertices can occur along the edges of the cube but we also connect opposite corners through diagonals that cross through the cube centroid. We call this topology a Cross-Diagonal Connected Cube (CDCC). Each vertex then has 4 outgoing lines of communication to its adjacent corners and opposite along the diagonal.

The CDCC topology requires hard-wired routing tables between nodes which we constructed using a variety of scripts. We have 32 gigabit lines on each vertex provided by the second network port on each compute node and 4 vertex connections. We therefore assign 8-gigabit lines to each vertex connection. These lines connect 8 compute nodes on one vertex to another 8 on a neighbouring vertex through direct port to port connections with a straight-thru cable. In this way, there is an 8-gigabit pipeline connecting each vertex which must be shared between 32 compute nodes. Also, the 4-port trunks connecting pairs of switches are 4- gigabit pipelines relaying traffic for 16 compute nodes per switch.

In this way, the CDCC topology with 256 nodes approximately represents a fully-switched network running at 250 Mbit, although it is slightly better since nodes on the same switch are fully-switched at 1 Gbit. The full-duplex cross-sectional bandwidth is 128 Gbit for 256 nodes. When constructing the routing tables, the maximum number of network hops to get from one node to another is 4 so latency is greater than a fully-switched system but we find in practice for our applications and standard benchmarks this is not a large hindrance. The performance of this configuration is about half the 500 Mbit performance that is expected for large gigabit switches that typically connect 16 ports to an 8- gigabit backplane internally. Latency is also about 4 times greater. However, the total cost for the 17 required 24-port switches plus cabling is much less than the cost of a single 256-port capable switch. We show below that our cluster provides competitive benchmark speeds for similar clusters with better networking and performs well for our applications.

### 3.3 Cabling

The task of cabling the cluster also proved to be formidable but we came up with some ways to minimize confusion and tangling. Connections from the first ethernet port to the switches is straightforward if the switches are centralized to one location at the centre of the rack array in a stack. Network cables are run over and under the cluster in small bundles and colour coded by rack number and then plugged into their appropriate ports in the central stack of switches. Although tedious, this stage of the cabling took a small team of people two days to complete. We point out that the same cabling effort is required even for more expensive multi-port switches so there is no significant overhead in manual labour.

In the second stage, it was necessary to wire the CDCC network using the point to point connections through the second network port on each node. Wiring these connections could potentially lead to a tangled mess of cables if the nodes are ordered sequentially according to their vertex numbers. A simple solution to this wiring problem is to re-order the nodes on the racks in logically connected pairs so that they are also physically connected to one another. In this way, only a short 1 foot length of cable is required for each contiguous pair and wiring can be done quickly. The logic of the connectivity of the CDCC network was used to create the correct ordered list of nodes on the rack for this arrangement as well as a wiring diagram for the connecting cables to the right switches in the fat-tree diagram. The only extra level of complexity for doing things this way in comparison to 200+ port monster switches is the node ordering and the need to make sure each cable is connected to the correct port in the switch stack. In practice, we have re-wired the network a couple of times to test the performance of different topologies and it generally only takes a few hours to do.

### 3.4 Routing

Another essential (and complex) part of getting the CDCC network to perform was the configuration of the routing tables for each Linux server. Network packets relayed through each server on this server have to know where to go according to our defined network topology.

A useful part of our design is that we have a back-up maintenance network (the fat-tree) that we can always use to gain access to each node in case of routing bugs or broken links in the CDCC network. This allowed us to experiment extensively with the scripts for generating the routing tables of the CDCC and led to an optimized system. Our strategy is to generate purely static routing tables on each server with a simple script at startup. Each of these scripts contains routing commands that establish about 50 static routes to individual nodes and sub-networks (vertices) in the CDCC network. Each node has its own unique routing table. We find this system works very well in practice and can be taken up and down quickly if necessary.

Our arrangement is also very fault-tolerant in case of node failures or broken links. If one or more nodes fails in the system, it will create holes in the CDCC network but as a contingency we can route around the broken nodes simply by using the fat-tree network. We have written a simple network repair script that pings all neighbours on the CDCC - if a node fails to answer we attempt to establish a direct route to it through the fat-tree - if that fails the node is assumed to be dead. In this way, we can continue to run jobs on the cluster even when main compute nodes go down. Either the hot spares, or development nodes can fill in while we replace the failed nodes. Since our networking is similarly fragmented into many switches, the hardware failure of one switch will only take out those nodes connected to it so most of the cluster can still run while the broken switch is being replaced.

Figure 1 shows some final details of our routing scheme that are worth pointing out. It turned out that the network performance of the trunks was not as good as we expected. We therefore modified our routing scheme to avoid using the trunks as much as possible. Network traffic from one node to another on different vertices used cross links that avoid the trunks. The only network traffic going through the trunks is between the nodes on bonded pairs of switches.

### 3.5 HPL Benchmark

As an initial test of performance, we ran the High-performance Linpack Benchmark (HPL), a portable implementation for distributed-memory computers to measure the speed of our cluster and compare with others on the Top 500 list . This code uses the Message Passing Interface (MPI) libraries. The actual benchmark involves the inversion of the largest matrix that can be stored on your cluster. In our case, the matrix size was $`160000\times 160000`$ elements. We compiled the code using the Intel C compiler icc version 7.0 and linked to the Intel MKL math libraries and Kazushige Goto’s optimized BLAS libraries . We used the LAM version 6.6b1 as the MPI library. We achieved a sustained speed of 1.202 Teraflops about 48% the theoretical peak speed of 2.46 Teraflops for 256 Xeon 2.4 GHz chips. This ratio of sustained to theoretical peak speed is comparable to many of the machines quoted on the November 2002 Top 500 list including similar sized clusters running with proprietary networking such as Myrinet. Our ranking on the list would currently be number 34. Of course, it remains to be seen how we will rank on the next list released in a few months since these trends tend to evolve very quickly!

## 4 Astrophysical Applications

The Beowulf cluster Mckenzie was constructed for problems in computational astrophysics. We describe here several applications that are being used to investigate areas on many astronomical scales from planets to the entire universe. We also highlight the performance of the various codes we have started to use.

### 4.1 Cosmology

Cosmological N-body simulation is the main numerical tool used to study the formation of structure in the universe. Dubinski et al.  have developed a new parallel algorithm that combines the standard techniques of particle-mesh methods (PM)  with the Barnes-Hut (BH) oct-tree algorithm . The PM method is used to measure the long-range gravitational forces from particles within a periodic cube representing a typical volume in the universe. These forces are refined on the sub-mesh scale using the BH algorithm by building a grid of oct-trees within the simulation domain. The new code is called GOTPM for Grid-(of)-Oct-Trees-Particle-Mesh. The code uses slab domain decomposition to distribute the particles among the CPUs. During the PM phase, Poisson’s equation for the gravitational potential is solved using Fourier methods. We use the MPI implementation of FFTW  to do the FFT’s during this phase which assumes slabs of equal width. A grid of BH oct-trees is then constructed within each slab and used to refine the forces at short range. Force accuracies are typically 0.5% using this method. To ensure load balance, the slab widths are allowed to vary during the tree phase with the width determined by the amount of computational work in the previous timestep. A fair amount of communication is required to move particle data at each step since they must move back and forth from an equal width slab decomposition for the PM phase and variable width slabs for tree phase. Despite this high demand for bandwidth, communication overhead is about 25% even with our inexpensive networking configuration.

Figure 2 shows one snapshot from a cosmological simulation using a standard model ($`\mathrm{\Omega }=0.3`$ and $`\mathrm{\Lambda }=0.7`$) in a 200 $`h^1`$ Mpc (650 million light years) box. This simulation used a $`1536^3`$ mesh to solve Poisson’s equation using Fourier methods in single precision. A total of $`N=768^3`$ (453M) particles are used and simulated for 5000 timesteps. An animation of a fly-through of this dark matter simulation as it structure develops is available at our website . In a 256 processor run, the wallclock time per timestep grows from 80 seconds at the beginning of the simulation to about 200 seconds once clustering develops with CPU time being dominated by the use of trees. The spatial resolution of these simulations is set by the gravitational softening radius used to shut down the $`1/r^2`$ Newtonian force at small separations. For this simulation, we set the softening radius to 1 kpc (3000 light years) so we achieve a spatial resolution about 65 times greater than what could be achieved using a PM method by itself.

### 4.2 Galaxy Dynamics

We have also applied a pure parallelized N-body treecode  to study galaxy interactions. The main mode by which spiral galaxies are transformed into ellipticals is through merging. When galaxies interact, the ordered energy of rotation and bulk motion of the two spiral galaxies is redistributed in the form of randomly oriented orbits. This process only requires a few dynamical times and the final structure of the merger remnant resembles and elliptical galaxy.

To study this process at high resolution, we built two models of spiral galaxies composed of a disk, a central bulge and surrounding dark halo using Kuijken & Dubinski’s method . To make the calculation more interesting, we used initial conditions chosen to represent a future configuration of the Milky Way and the Andromeda Galaxy which will likely merge in about 3 billion years. Each model is composed of 128M particles representing stars and 25.6 million particles representing the surrounding dark matter halos for a total of 307.2M particles. The two galaxies are on a collision course and we computed the trajectories for all particles for 5300 timesteps representing a physical timescale of 2.3 billion years. A sequence of 9 images showing the evolution of the system over this timescale is shown in Figure 3. An animation showing the galaxy collision is available here .

The code used for these calculations is PARTREE , a parallelized treecode based on Salmon’s  original algorithm and ideas of locally essential trees. The code has recently been modified to incorporate asynchronous message-passing in the construction of the locally essential trees to minimize communication overhead which becomes dominant for large numbers of CPUs. With current modifications, communication overhead only used 10% of the wallclock time of 180 seconds per step in the 307M particle simulation. The code should scale well to thousands of processors in its current form assuming the problem size grows as well. A several thousand CPU machine with 2GB/CPU is now large enough to follow the trajectory of about 10 billion particles, the same number of stars in some galaxies and only a factor of 40 away from the number of stars in the Milky Way. Galactic dynamics simulations will reach their ultimate resolution in terms of the number particles representing stars in less than a decade.

### 4.3 Blackhole Accretion

We have undertaken MHD simulations with 1400<sup>3</sup> zones arrayed in a uniform Cartesian grid, the largest MHD simulations to date (Figure 4). At this resolution, each full dimensional sweep corresponding to two timesteps took 40 seconds. A series of optimizations allows the code to exploit the hardware’s vector units, hyperthreading OpenMP parallelism, and uses the Message Passing Interface (MPI) to communicate between nodes. The code  is based on a 2nd order accurate in space and time high resolution Total-Variation-Diminishing (TVD) algorithm; it explicitly conserves the sum of kinetic, thermal and magnetic energy; hence magnetic dissipation (at the grid scale) heats gas directly. No explicit resistivity or viscosity is added, and reconnection and shocks occur through the solution of the flux conservation laws and the TVD constraints. Magnetic flux is conserved to machine precision by storing fluxes perpendicular to each cell face.

Each time step requires about 20 seconds, corresponding to 2.2 Tflop sustained in single precision. The operation count was obtained by adding the number of operations in the source code by hand. This is about 44% of theoretical peak speed achievable with the using the SSE2 vector units. The main parallelization overhead is the need for 16 buffer zones between computational domains. This increases the floating point count by 65%. Fatter nodes with more processors and more memory would decrease the surface to volume ratio. The domain decomposition was matched to the physical three dimensional connectivity of the network. Most traffic does not require any routing, and up both gigabit ports per node can be effectively exploited. All message passing latency is hidden by staggering the computation phase: the processes do not need to wait for data, and communication is asynchronous.

### 4.4 Dusty disks and the Formation of planets

The earliest stages of the planet formation process, in which micron sized particles of dust somehow come together to form meter and kilometer sized objects and eventually planets, is very poorly understood. Humble, Maddison and Murray  have developed a new parallel two-phase (dust and gas) smoothed particle hydrodynamics  and self-gravity code, so that for the first time we can model the non-linear evolution of two coupled fluids in a 3D protoplanetary dusty disk.

A typical calculation uses 125k gas and 125k dust particles and resolves the motion of the entire disk from $`2`$ to $`200`$AU. We investigate how dust evolves in a so-called minimum mass solar nebula in which the total mass of the dust is $`1`$% of the gas phase, which in turn is $`1`$% the mass of the protostar. Figure 5 shows the distribution of $`1`$cm grains after $`10^4`$ years.

Our simulations clearly capture the expected settling and density enhancement of the dust in the midplane of the disk where grain growth is most likely to occur, and where gravitational instabilities may emerge. The simulations also reveal the new result that the pressure supported differential rotation between gas and dust is maintained only in a relatively narrow region in radius. Rapid evolution of the dust occurs mostly in this narrow band, the location of which is a strong function of the dust grain size. This results in large variations of structure and dynamics between dusty disks of different characteristic grain radii.

The code is based on the hashed oct tree algorithms of Salmon and Warren , is written in C and MPI and uses HDF5  for i/o. The code has evolved to be a flexible parallel tree framework, being tunable for minimum bandwidth or maximum speed, and incorporating asynchronous MPI latency hiding techniques using multiple simultaneous tree traversals from each processor. The hash table allows O($`1`$) random access to all data in the tree and acts as a serialisation point for remote data retrieval, whilst more traditional pointers are used for rapid traversal of all locally cached tree data.

Most calculations performed to date utilise around $`32`$ processors with $`90`$% parallel efficiency. This number of processors is currently adequate as this is a new field and there are large regions of parameter space to explore. A relatively inexpensive machine like the McKenzie cluster is vital so that many runs with different configurations and disk physics can be explored simultaneously in a reasonable time. We are working on further scalability and speed enhancements of the code.

## 5 Conclusions

We have built Mckenzie, 512-CPU Beowulf cluster that incorporates commodity gigabit networking in a new topology that permits 128 Gbit cross-sectional bandwidth. Networking costs represent 9% the total cost of the cluster of $561K. The machine has been benchmarked with the Top 500 list HPL code and achieved a sustained speed of $`1.202`$ Teraflops or 48% of the theoretical peak speed making it competitive with clusters with more expensive networking solutions. By this measure, our computing costs are just $0.47/Megaflop. We note that the idea of using a Linux node as a router through a secondary network or even tertiary network is a general concept that can be applied to different network topologies with larger numbers of nodes. The main limitation of this method is not so much the bandwidth but the extra latency inherent to commodity switches and perhaps multiple hops in the chosen topology. Nevertheless, complex parallel algorithms with high bandwidth requirements can still run on systems like this based on our experience.

We have also described 4 astrophysical parallel applications running on the cluster that span the scale of the universe from solar system scale to the size of the universe. We have been able to run some of the largest simulations to date in cosmology, galaxy dynamics and MHD accretion flows using MPI parallel codes. N-body simulations of hundreds of millions of particles in cosmology and galaxy dynamics can be done routinely. The MHD code achieved a sustained speed of 2.2 teraflops in single precision or 44% of the theoretical peak speed in $`1400^3`$ zone calculation. Our low-cost networking also justifies using the machine in a parameter survey mode on smaller sub-units of the cluster while such strategies might be considered a waste of expensive resources on a large-scale SMP. The ability to run on 512 processors has allowed us to tune and optimize our codes to hide latency for large processor numbers and paves the way for porting codes to clusters of thousands of processors.

## Acknowledgments

We acknowledge the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Canadian Foundation for Innovation (CFI) for funding this project.


text:
# Electronic structure of nuclear-spin-polarization-induced quantum dots

## I INTRODUCTION

The theoretical and experimental researches of quantum dots have attracted much attention in recent years . Quantum dots are usually fabricated experimentally by applying lithographic and etching techniques to impose a lateral structure onto an otherwise two-dimensional electron system. Lateral structures introduce electrostatic potentials in the plane of the two-dimensional electron gas, which confines the electrons to a dot region. The energy levels of electrons in such quantum dots are fully quantized like in an atom. In such electrically confined quantum dots the confining potential can be well represented by a parabolic potential.

Another method of low-dimensional structure fabrication consists of the application of spatially inhomogeneous magnetic fields. There has been proposed several alternative magnetic structures subsequently realized experimentally. Among them: magnetic dots using a scanning tunneling microscope lithographic technique , magnetic superlattices by the patterning of ferromagnetic materials integrated by semiconductors , type-II superconducducting materials deposited on conventional heterostructures , and nonplanar two-dimensional electron gas (2DEG) systems grown by a molecular beam epitaxy . Such systems were studied theoretically in a series of papers by different authors .

In the present paper we study a quantum dot system which is different from the quantum dot systems discussed above: (1) the electrons are confined through local nuclear spin polarization, (2) the confinement potential is inherently nonparabolic and time-dependent, it is a solution of the diffusion equation when considering relaxation, and (3) the dot contains electrons with only one spin direction. Such system was proposed for the first time in Ref. . However, the properties of Nuclear-Spin-Polarization-Induced Quantum Dots (NSPIQD) have not been considered thus far and this is the motivation behind the present investigation. In our calculations we use some ideas from , where a nuclear-spin-polarization-induced quantum wire was proposed and investigated.

Electron and nuclear spins interact via the contact hyperfine interaction. Once the nuclear spins are polarized, the charge carrier spins feel the effective hyperfine field, $`𝐁_{hf}`$, which lifts the spin degeneracy. The maximum nuclear field in GaAs can be as high as $`B_{hf}=`$5.3T in the limit that all nuclear spins are fully polarized . This high level of nuclear spin polarization has been achieved experimentally. For example, the optical pumping of nuclear spins in 2DEG has demonstrated nuclear spin polarization on the order of $`90\%`$, . A similarly high polarization has been created by quantum hall edge states ($`85\%`$) . The spin splitting due to such a hyperfine magnetic field is comparable to the Fermi energy of 2DEG. It is important to note that the hyperfine field does not manifest itself magnetically due to the smallness of the nuclear magnetic moments. The electrons in the region where nuclear spins are polarized will preferably occupy the energetically more favorable states with the spins opposite to $`𝐁_{hf}`$. Furthermore, the nuclear polarization acts on the electrons as the effective confining potential. This effective confining potential can be used to create different nanostructures with polarized electrons in them. In this paper we consider a nuclear-spin-polarization-induced quantum dot (NSPIQD).

The proposed system is depicted in Fig. 1. The nuclear spins are polarized homogeneously along the z-axis perpendicular to the 2DEG in heterostructure by any other suitable experimental method. For example, the optical nuclear spin polarization  or the transport polarization  can be used. The region where the nuclear spins are polarized is indicated by the cylinder in Fig. 1. The NSPIQD is created in the region of intersection of the 2DEG with the region of local nuclear spin polarization. The gate electrode below the 2DEG is used to control the number of electrons in the NSPIQD. Moreover, the system is subjected to an external magnetic field along the z-axis. The magnetic field plays an important role in the nuclear spin polarization process and, under specific conditions, increases the nuclear spin relaxation time. Assuming a small magnetic field, we can neglect it in our calculations, focusing on the effects caused by the confining hyperfine field. Our paper is organized as follows. In Sec. II we discuss the properties of nonhomogeneous nuclear spin polarization and calculate the evolution of initially-created hyperfine-field profile which is taken, for simplicity, in the Gaussian form. Time dependence of the electron states in NSPIQD is studied in Sec. III. The conclusions of this investigation are presented in Sec. IV.

## II HYPERFINE-FIELD PROFILE

Let us assume that the method of optical nuclear spin polarization is used  to create a NSPIQD. To pattern a nanostructure it is proposed to illuminate the sample locally by, for example, putting a mask on it. The usual optical technique allows one to create the light beams of the width of the order of the wavelength ($`500nm`$). By using near fields optics the beam width can be sufficiently reduced ($`100nm`$). Hence a $`1\mu m`$-size NSPIQD can be easily created by the modern experimental technique.

There are two main mechanisms leading to the time dependence of the hyperfine field: the nuclear spin relaxation and the nuclear spin diffusion. We assume that the initial nuclear spin polarization is homogeneous in the $`z`$-direction. Then the hyperfine field evolution is described by the two-dimensional diffusion equation:

$$\frac{B_{hf}}{t}=D\mathrm{\Delta }B_{hf}\frac{1}{T_1}B_{hf}\text{,}$$
(1)
accounting for the relaxation processes. Here $`D`$ is the spin-diffusion coefficient, $`\mathrm{\Delta }`$ is two-dimensional Laplace operator, and $`T_1`$ is the nuclear spin relaxation time . The formal solution of Eq.(1) can be written as

$$B_{hf}=e^{\frac{t}{T_1}}G(𝐫𝐫^{},t)B_{hf}\left(𝐫^{},t=0\right)𝑑𝐫^{}.$$
(2)
Here $`G(𝐫,t)=\frac{e^{\frac{\left(𝐫𝐫^{}\right)^2}{4Dt}}}{4\pi Dt}`$ is the Green function of the diffusion equation and $`B_{hf}\left(𝐫^{},t=0\right)`$ is the initial hyperfine field profile.

In this paper we consider NSPIQD having the cylindrical symmetry; that is, the hyperfine field $`B_{hf}`$ is a function of $`r`$. In the simplest case, we can assume the initial condition to be of the Gaussian form: $`B_{hf}(r,0)=B_0\mathrm{exp}\left(\frac{r^2}{2d^2}\right)`$. The two parameters, $`d`$ and $`B_0`$, define the half-width and the amplitude of the initial distribution of the hyperfine field, respectively. Then the solution of Eq. (1) is:

$$B_{hf}(r,t)=B_0e^{\frac{t}{T_1}}\left(1+\frac{t}{t_0}\right)^1e^{\frac{r^2}{2d^2\left(1+\frac{t}{t_0}\right)}}\text{ ,}$$
(3)
where $`t_0=\frac{d^2}{2D}`$. The value of $`t_0`$ is the time it takes for the $`B_{hf}(0,t)`$ to reduce by factor of two from $`t=0`$ due to the nuclear spin diffusion. The nuclear-spin relaxation time, $`T_1`$, in semiconductors at sufficiently low temperatures is rather long. It varies from several hours to a few minutes . The available experimental values for the diffusion coefficient are $`D10^{13}`$ cm2s-1 for 75As in bulk GaAs  and $`D=10^{14}`$ cm2s-1 in Al0.35Ga0.65As . For $`d=1`$ and 5 $`\mu `$m taking $`D=10^{13}`$ cm2s-1 we have $`t_0=5\times 10^4`$, $`1.25\times 10^6`$ s.

## III ENERGY SPECTRUM

The microscopic description is based on the following Hamiltonian:

$$H=\frac{\mathrm{}^2}{2m^{}}\mathrm{\Delta }+\frac{1}{2}g^{}\mu _B\sigma 𝐁_{hf}(r,t)+U\left(z\right)$$
(4)
where $`m^{}`$ is the electron effective mass, $`g^{}`$ is the effective electron $`g`$-factor ($`g_{GaAs}^{}=0.44`$), $`\mu _B`$ is the Bohr magneton, $`\sigma `$ is the vector of Pauli matrices, $`𝐁_{hf}`$ is given by Eq. (3) and $`U(z)`$ is the 2DEG confining potential. We suppose, as is usually done for the 2DEG, that only the lowest sub-band, corresponding to the confinement in $`z`$-direction, is occupied and we can ignore the higher sub-bands. Thus, we omit $`z`$-dependence of the wave function in the following. The time scale introduced by a nuclear spin system is several orders of magnitude larger than the time scale of typical electron equilibration processes. In such a case the conduction electrons see a quasi-constant nuclear field. This simplifies calculation by avoiding the complications which would appear when solving the Schrödinger equation with the time dependence due to polarized nuclei. We take into account the electrons of only one spin direction (for which the effective potential is attractive).

The one-electron eigenvalue problem with the attractive Gaussian potential (Eq. (3)) does not admit analytical solutions. Different approximate methods  were implemented to solve this problem. In the present paper, an analytical solution of the Schrödinger equation is found within the parabolic approximation of the hyperfine field :

$$\stackrel{~}{B}_{hf}=abr^2$$
(5)
connected with Eq. (3) by the relations:

$$\stackrel{~}{B}_{hf}(0,t)=B_{hf}(0,t)$$
(6)
and

$$_0^{r_0}r\stackrel{~}{B}_{hf}(0,t)𝑑r=_0^{\mathrm{}}rB_{hf}(0,t)𝑑r.$$
(7)
Here $`r_0^2=a/b`$. Eq. (6) connects the depth of potentials, Eq. (7) provides equal nuclear-spin polarization for the two fields. From Eqs. (6,7) we obtain $`a=B_0\frac{e^{\frac{t}{T_1}}}{1+\frac{t}{t_0}}`$ and $`b=B_0\frac{e^{\frac{t}{T_1}}}{2d^2\left(1+\frac{t}{t_0}\right)}`$. The energy spectrum for the parabolic potential (5) in units of $`E_0=\frac{\mathrm{}^2}{2m^{}d^2}`$ is given by

$$\epsilon _{n,m}=\frac{g^{}\mu _BB_0}{2E_0}\frac{e^{\frac{t}{T_1}}}{1+\frac{t}{t_0}}+\sqrt{\frac{g^{}\mu _BB_0}{E_0}}\frac{e^{\frac{t}{2T_1}}}{1+\frac{t}{t_0}}\left(2n+|m|+1\right)$$
(8)
where $`n=0,1,\mathrm{}`$ and $`m=0,\pm 1,\mathrm{}`$ .

The exact solution of the Schrödinger equation with the Gaussian profile of the hyperfine field (Eq. (3)) was found numerically. Due to the cylindrical symmetry of the problem, the wave function can be written as $`\psi (\rho ,\varphi )=\frac{1}{\sqrt{2\pi }}e^{im\varphi }R(\rho )`$. The equation for the radial part $`R(r)`$ of wave function has a form

$$\left[\frac{1}{x}\frac{d}{dx}x\frac{d}{dx}\frac{m^2}{x^2}+\gamma \frac{B_{hf}(x,t)}{B_{hf}(0,0)}+\epsilon _{n,m}\right]R_{n,m}=0,$$
(9)
where $`x=r/d`$ is the dimensionless coordinate and $`\gamma =\frac{g^{}\mu _BB_{hf}(0,0)}{2E_0}`$. For $`d=1`$ and $`5`$ $`\mu `$m, taking $`m^{}=0.067m_e`$, we have $`E_0=0.57\times 10^3`$, $`0.023\times 10^3`$ meV; for $`B_{hf}(0,0)=2.65`$ (50% nuclear spin polarization) and $`5.3`$ T (100% nuclear spin polarization) corresponding energies are $`\frac{1}{2}g^{}\mu _BB_{hf}(0,0)=3.4\times 10^2`$ and $`6.8\times 10^2`$ meV. We have used the Shooting Method to solve Eq. (9), subjecting the solution to the following boundary conditions: $`R_{n,m}\left(\rho 0\right)=\rho ^{|m|}`$ and $`R_{n,m}\left(\rho \mathrm{}\right)=0`$. The results of the numerical calculations are presented below.

The time-dependence of the electron energy levels in the NSPIQD is determined by the time-dependence of the confining hyperfine field. There are two characteristic times in the problem: the diffusion characteristic time $`t_0`$ and the relaxation characteristic time $`T_1`$. We can distinguish the diffusive regime, when $`tt_0T_1`$, the intermediate regime, $`tt_0T_1`$ and the relaxation regime, $`tT_1t_0`$ . Here $`t`$ is the observation time.

Figure 2 shows the time dependence of the electron energy levels for the Gaussian and parabolic potentials in the diffusion regime. We emphasize that the parabolic potential can be regarded as a good approximation of the Gaussian potential only for the ground state. The excited-state energy levels for the parabolic potential reveal large deviations from those for the Gaussian potential, which manifest in the degeneracy of states and in the shift of levels. This result is qualitatively similar to those obtained for 3D Gaussian and parabolic potential . However, time-dependence of energy levels for both potentials show quite similar behavior. The number of energy levels in NSPIQD remains constant, whereas their depth decreases. From Eq. (8) it follows that in the diffusion regime the time-dependence of the energy levels in the parabolic potential is $`\epsilon _{n,m}(t)=\frac{\epsilon _{n,m}(0)}{1+\frac{t}{t_0}}`$. It can be shown that the energy levels in the Gaussian potential have the same time-dependence. Substituting Eq. (3) into Eq. (9) and introducing the variable $`\xi `$ as $`x=\xi \sqrt{1+t/t_0}`$ we obtain:

$$\left[\frac{1}{\xi }\frac{d}{d\xi }\xi \frac{d}{d\xi }\frac{m^2}{\xi ^2}+\gamma e^{\frac{\xi ^2}{2}}+\left(1+\frac{t}{t_0}\right)\epsilon _{n,m}\right]R_{n,m}=0.$$
(10)
The time-dependent factor, $`\left(1+\frac{t}{t_0}\right)`$, appears in Eq. (10) only as a product with $`\epsilon _{n,m}`$ thus proving the statement.

Figures 3 and 4 show the results obtained for the intermediate and relaxation regimes. On the contrary, the number of the energy levels in NSPIQD decreases in time in these regimes. This decrease occurs on the scale of $`T_1`$. We can not explicitly obtain time-dependence of energy levels for the Gaussian potential in these regimes. The parabolic approximation of the hyperfine field serves as a good approximation again only for the ground energy level. The evolution of excited-state energy levels in the Gaussian and in the parabolic potentials are different: the lifetimes of energy levels obtained in the case of the parabolic potential are shorter than in the case of the Gaussian potential.

It is important to know the lifetime of the NSPIQD. We can consider electron states in the NSPIQD up to the moment when the confining potential depth is more than the temperature. Consequently, the lifetime $`t_l`$ of the NSPIQD can be defined by the following condition: $`\frac{\left|g^{}\mu _BB_{hf}(0,t_l)\right|}{2}=k_BT`$, where $`k_B`$ is the Boltzman constant and $`T`$ is the temperature. Using Eq.(3), we calculate time $`t_l`$ for two limiting cases: $`T_1t_0`$ and $`T_1t_0`$. In the first case (the relaxation regime), $`t_lT_1\mathrm{ln}\frac{\left|g^{}\mu _BB_0\right|}{2k_BT}`$. In the second case (the diffusion regime), $`t_lt_0\left|\frac{g^{}\mu _BB_0}{2k_BT}1\right|`$. Time-dependence of the half-width of NSPIQD is $`d(t)=d\sqrt{1+t/t_0}`$. Let us estimate it at $`t=t_l`$. For $`T^{}=30mK`$ and $`B_0=2.65`$T we have $`d(t)=d`$ in the relaxation regime and $`d(t)=3.6d`$ in the diffusion regime.

## IV CONCLUSIONS

We have studied the electron energy levels of a NSPIQD created in the region of the intersection of a local nuclear spin polarization with a 2DEG. The properties of the NSPIQD are time-dependent because of the nuclear spin diffusion and relaxation. There are two characteristic time and three corresponding regimes: the diffusion regime, the intermediate regime and the relaxation regime. In the diffusion regime, the number of electron energy levels remains constant with time. In the relaxation and intermediate regime, the number of electron energy levels decreases with time. Time-dependence of the electron energy levels in the diffusion regime has a simple form. Since the characteristic relaxation time is proportional to the square of the NSPIQD radius at $`t=0`$, it is possible to create NSPIQDs operating in different regimes using the same experimental setup.

The numerical estimations allow us to conclude that the system under study can be realized experimentally. For a hyperfine field of just a few teslas, the experiment could be made at a temperature of the order of $`10`$mK. The modern experimental technique allows one to create a region with local nuclear spin polarization of characteristic sizes $`100`$nm, making the NSPIQD having a small size. The spectroscopy of the NSPIQD could be used to obtain some information about nuclear spin interactions in solids, for example, the nuclear spin relaxation time and the nuclear spin diffusion coefficient.

It should be pointed out that a simplified model was used in this paper to describe the single electron states in the NSPIQD. We considered the influence of a nuclear spin-related hyperfine field on the electron states, whereas the electrons could also alter the nuclear spin dynamics. The well-known examples of such phenomena are the indirect long-range nuclear-spin interaction, electron-assisted mechanisms of nuclear spin relaxation and nuclear spin precession in an effective field created by the electrons . Investigation of these effects is beyond the scope of this paper. Results of such investigations will be published elsewhere.

## ACKNOWLEDGMENTS

We acknowledge useful discussions with V. Privman, S. Shevchenko, and I. Vagner. This research was supported by the National Security Agency and the Advanced Research and Development Activity under Army Research Office contracts DAAD-19-99-1-0342 and DAAD 19-02-1-0035, and by the National Science Foundation, grants DMR-0121146 and ECS-0102500.


text:
# Chapter 1 Radio Emission From EAS - Coherent Geosynchrotron Radiation

## Chapter 1 Radio Emission From EAS - Coherent Geosynchrotron Radiation

### Abstract

Extensive air showers (EAS) have been known for over 30 years to emit pulses of radio emission at frequencies from a few to a few hundred MHz, an effect that offers great opportunities for the study of EAS with the next generation of “software radio interferometers” such as LOFAR and LOPES. The details of the emission mechanism, however, remain rather uncertain to date. Following past suggestions that the bulk of the emission is of geomagnetic origin, we model the radio pulses as “coherent geosynchrotron radiation” arising from the deflection of electrons and positrons in the earth’s magnetic field. We analytically develop our model in a step-by-step procedure to disentangle the coherence effects arising from different scales present in the shower structure and infer which shower characteristics govern the frequency spectrum and radial dependence of the emission. The effect is unavoidable and our predictions are in good agreement with the available experimental data within their large margins of error.

### 1. Introduction

Research concerning radio emission from EAS almost ceased completely in the 1970s. With the advent of fully digital “software radio interferometers” such as LOFAR, however, it once more becomes very attractive to measure EAS through their radio emission. The technique is also ideally suited to be combined with classical particle detector arrays such as KASCADE Grande (as evident from the LOPES experiment described in ) or AUGER and offers the advantage of a very high duty cycle and cost-effectiveness over air fluorescence methods.

Past experiments have established that the dominant emission mechanism is related to the earth’s magnetic field. The strength of the emission is, however, still rather unclear and past modeling efforts have not been developed to sufficient depth for application to concrete experiments. We therefore take the new approach of modeling the radio emission as “coherent geosynchrotron emission” as first proposed by  and developed in depth in  and . Compared with recent Monte-Carlo simulations by  we take a more analytical approach in hope to better understand the various coherence effects shaping the emission and additionally take into account the longitudinal shower evolution.

### 2. The Model

Starting with the works of  we base our model on the frequency spectra of synchrotron pulses emitted by highly relativistic electron-positron pairs:

$$\stackrel{}{E}(\stackrel{}{R},\omega )=\left(\frac{4\pi }{c}\right)^{1/2}\frac{1}{R}\frac{\omega e}{\sqrt{8c}\pi }\mathrm{e}^{\mathrm{i}(\omega \frac{R}{c}\frac{\pi }{2})}\left(\stackrel{}{\widehat{e}}_{}\right)A_{}(\omega ),$$
(1)
where

$$A_{}(\omega )=\mathrm{i}\frac{2\rho }{\sqrt{3}c}\left(\frac{1}{\gamma ^2}+\theta ^2\right)K_{2/3}(\xi )$$
(2)
with

$$\xi =\frac{\omega \rho }{3c}\left(\frac{1}{\gamma ^2}+\theta ^2\right)^{3/2}\mathrm{and}\rho =\frac{\gamma m_ec^2}{eB\mathrm{sin}\alpha }.$$
(3)
Our gradual development of the integration over the shower particles helps in understanding the coherence effects arising from the different scales present in the shower. Examination of a point source demonstrates the importance of the intrinsic beaming of the synchrotron pulses for the radial dependence of the emission. Going from a point source to a longitudinal line charge of a few metres length (corresponding to the thickness of the shower “pancake”) shows that the associated coherence losses limit the frequency spectrum to the regime $`<100`$ MHz. Extending the model to a spherical shell and later on a “flaring disk” reveals coherence effects arising from a more realistic shower geometry. We take the lateral particle distribution into account via NKG-functions and incorporate the longitudinal spread as a function of distance from the shower axis through empirical data. An enhancement of the emission at high distances in the E-W compared to the N-S direction arises from the dependence on the geomagnetic field. In the last step we integrate over the longitudinal shower evolution as a whole, which has a strong impact on the total emission strength and the radial dependence of the emission. We also consider the energy distribution of the particles.

### 3. Comparison with Data

At the moment, the available historic data on radio emission from EAS is subject to large uncertainties. This is mainly due to discrepancies in the emission strength of an order of magnitude between the data sets from different groups, probably arising from radio calibration issues, and a lack of precision in the documentation of the more than 30 years old data sets. Additionally, there are some subtleties involved in the conversion of theoretical values ($`E_\omega `$) to measured quantities ($`ϵ_\nu `$). We have chosen the well documented data of  and compared them with our predictions for the radial dependence of the emission of a $`10^{17}`$ eV vertical air shower in Figure 1. For an assessment of the spectral dependence as shown in Figure 2 we have taken data from  and  and rescaled them so that they are consistent with the data of .

### 4. Discussion

As can be seen in Figures 1 and 2, our calculations successfully reproduce the trends visible in the data sets. As for the absolute strength of the emission, our model somewhat over-predicts the data of , which is already at the upper end of the available data sets. Given the very simplified integration over the longitudinal shower evolution this is, however, not too surprising. We also find a significant E-W to N-S anisotropy in the emission strength (not shown here). In general, it is thus clear that a significant contribution by the geomagnetic emission mechanism is unavoidable. At the moment, the acquisition of reliable experimental data is of paramount importance. As our signal-to-noise calculations indicate, this will be accomplished by LOPES in the near future, see .

### 5. Conclusions

Our model of “coherent geosynchrotron radiation” reproduces the available experimental data within their large margins of error. The step-by-step modeling procedure illustrates that the high-frequency cutoff in the emission spectra is mainly governed by the longitudinal spread, i.e. thickness, of the shower “pancake”. The radial dependence is governed by the intrinsic beaming of the pulses as well as the integration over the shower evolution as a whole. The predicted asymmetry of emission strength in N-S and E-W directions will be verifiable with LOPES and would help in the detection of EAS out to high distances in the E-W direction. In a next step, we will use this model as a basis for the development of a sophisticated numerical simulation that includes additional aspects such as near-field effects, polarisation and Askaryan-type Čerenkov radiation.

### 6. List of Symbols/Nomenclature

| $`\stackrel{}{R}`$ = distance observer-particle | $`\stackrel{}{E}`$ = electric field |
| --- | --- |
| $`c`$ = speed of light | $`e`$/$`\mathrm{m}_e`$ = electron unit charge/mass |
| $`\omega `$ = 2$`\pi `$ observing frequency | $`\theta `$ = minimum angle to line of sight |
| $`\gamma `$ = particle Lorentz factor | $`K_\nu `$ = modified Bessel-function of order $`\nu `$ |
| $`B`$ = earth’s magnetic field | $`\alpha `$ = pitch angle of particle trajectory |

### 7. References

1. Allan H. R., Clay R. W., Jones J. K. 1970, Nature 227, 1116

2. Falcke H., Gorham P. W. 2003, Astropart. Phys. in press, astro-ph/0207226

3. Horneffer A., Falcke H. 2003, Proc. of the 28th ICRC, these proceedings

4. Huege T., Falcke H. 2003, A&A submitted

5. Jackson J. D. 1975, Classical Electrodynamics (J. Wiley & Sons, New York)

6. Prah, J. H. 1971, M.Phil. thesis, University of London

7. Spencer R. E. 1969, Nature 222, 460

8. Suprun D. A., Gorham P. W., Rosner J. L. 2003, astro-ph/0211273