format: text

text:
# Some equivalences for Martin’s Axiom in asymmetric topology

## Abstract.

We find some statements in the language of asymmetric topology and continuous partial orders which are equivalent to the statements $`\kappa <𝔪`$ or $`\kappa <𝔭`$.

###### Key words and phrases:

Martin’s Axiom, locally compact, core compact, sober, ccc, continuous lattice  Bruce S. Burdick, Some equivalences for Martin’s Axiom in asymmetric topology, Proceedings of the Ninth Prague Topological Symposium (Prague, 2001), pp. 9–13, Topology Atlas, Toronto, 2002

We think of asymmetric topology as those parts of topology in which the specialization order, $`xy`$ if and only if $`xc\{y\}`$, need not be symmetric. (See  for some of the motivations.)

Martin’s Axiom has many equivalent statements, consequences, and variations in the literature which can be stated in topological terms. Most of the treatments we have seen so far from set-theoretic topologists assume that spaces are Hausdorff. In view of recent interest in asymmetric topology, in which even $`T_1`$ spaces are a highly symmetric special case, we give some equivalences for Martin’s Axiom which utilize the terms of this field.

Our reference for properties related to Martin’s Axiom is , and for properties related to continuous lattices we referred to .

###### Definition 1.

A partially ordered set $`(𝒫,)`$ is upwards-ccc if any uncountable subset of $`𝒫`$ must have two distinct members which have a common upper bound in $`𝒫`$. The cardinal $`𝔪`$ is the least cardinal such that there exists a non-empty upwards-ccc partially ordered set $`(𝒫,)`$ and a collection $`\{D_\alpha |\alpha <𝔪\}`$ of cofinal subsets of $`𝒫`$ such that no upwards-directed subset of $`𝒫`$ meets each $`D_\alpha `$.

It can be shown that $`\omega _1𝔪𝔠`$. The Martin’s Axiom of the title is the statement $`𝔪=𝔠`$.

###### Definition 2.

A topological space is ccc if any uncountable collection of open sets has two distinct members which are not disjoint. A space is locally compact if every open set contains a compact neighborhood of each of its points. For two sets $`U`$ and $`V`$, we say $`U`$ is compact in $`V`$, denoted $`UV`$, if for every open cover of $`V`$ there is a finite subcollection which covers $`U`$. A space is core compact if every open set $`O`$ is the union of open sets $`O^{}`$ with $`O^{}O`$.

###### Definition 3.

A subset $`C`$ of a space is irreducible if it is closed and the non-empty open sets in the subspace topology on $`C`$ form a filterbase. A space is sober if it is $`T_0`$ and every irreducible subset is the closure of a point. A space is supersober if it is $`T_0`$ and every ultrafilter on the space has either the empty set or the closure of a point as its set of limits. A filter $``$ of open sets on a space is Scott open if whenever $`𝒟`$ is a collection of open sets which is directed upward under inclusion, $`𝒟`$ implies that some element of $`𝒟`$ is a member of $``$.

We will make use below of the Hofmann-Mislove Theorem, that on a sober space the intersection of a Scott open filter of open sets is non-empty. (See  for the original reference, and see  for a new proof of the form of the theorem that we have just stated.)

###### Definition 4.

In a partially ordered set $`(𝒫,)`$, for two elements $`x`$ and $`y`$ we say $`x`$ is way below $`y`$, denoted $`xy`$, if for any directed subset $`D`$ of $`𝒫`$ with supremum $`supD`$, if $`ysupD`$ then there is some $`dD`$ with $`xd`$. $`(𝒫,)`$ is upwards-continuous if for all $`y𝒫`$ the set $`\{x𝒫|xy\}`$ is a directed set with supremum $`y`$.

###### Definition 5.

A lattice is distributive if it satisfies the distributive laws,

$$a(bc)=(ab)(ac),$$

$$a(bc)=(ab)(ac).$$

A complete lattice is downwards-ccc if in any uncountable subset there are two distinct members $`x`$ and $`y`$ with $`xy0`$, where 0 is the least element. An element $`x`$ in a complete lattice is a non-zero divisor if $`xy0`$ for any $`y0`$. An element $`x`$ is irreducible if whenever $`x=ab`$ then $`x=a`$ or $`x=b`$.

###### Proposition 1.

Let $`\kappa `$ be a cardinal. The following are equivalent.

1. $`\kappa <𝔪`$
2. In a ccc locally compact sober space the intersection of $`\kappa `$ or fewer open dense sets is non-empty.
3. In a ccc core compact space if $`𝒟`$ is a collection of $`\kappa `$ or fewer open dense sets then there is an irreducible set meeting every member of $`𝒟`$.
4. Let $``$ be a downwards-ccc upwards-continuous distributive complete lattice. If $`D`$ is a set of $`\kappa `$ or fewer non-zero divisors of $``$ then there is an irreducible $`p`$ such that no element of $`D`$ is below $`p`$.

###### Proof.

(1) implies (2). Given a ccc locally compact sober space $`(X,𝒯)`$ let $`𝒫=𝒯\{\varphi \}`$ and define a partial order on $`𝒫`$ by $`O_1O_2`$ if either $`O_1=O_2`$ or $`O_2`$ is compact in $`O_1`$. By local compactness, two open sets $`O_1`$ and $`O_2`$ have a common upper bound in $`(𝒫,)`$ if and only if $`O_1O_2\varphi `$. So $`(𝒫,)`$ is upwards-ccc. Given a collection $`\{O_\alpha |\alpha <\kappa \}`$ of open dense subsets of $`X`$, for each $`\alpha `$ define $`D_\alpha =\{O𝒫|OO_\alpha \}`$. Then each $`D_\alpha `$ is cofinal in $`(𝒫,)`$, again by local compactness. Since $`\kappa <𝔪`$ there is a directed set $`S𝒫`$ that meets each $`D_\alpha `$. We may assume that $`S`$ is closed upwards under $``$, i.e., $`S`$ is a filter of open sets on $`(X,𝒯)`$.

If $`S\varphi `$ we are done since an element of $`S`$ must be contained in some member of each $`D_\alpha `$ and so is in $`_{\alpha <\kappa }O_\alpha `$. If $`S`$ has a maximal element under $``$ then it is $`S`$ and we are done. Otherwise $`S`$ is Scott open, since if $`𝒟`$ is a directed collection of open sets with $`𝒟S`$ then some element $`OS`$ is strictly greater than $`𝒟`$, which means it is compact in $`𝒟`$, and so it is a subset of some element of $`𝒟`$. By the Hofmann-Mislove Theorem $`S\varphi `$.

(2) implies (4). Given a downwards-ccc upwards-continuous distributive complete lattice $``$, by Theorem 1.2 in  there exists a locally compact sober space $`(X,𝒯)`$ such that $``$ is order isomorphic to $`𝒯`$ with inclusion. We note that this $`(X,𝒯)`$ must be a ccc space. Given $`\{d_\alpha |\alpha <\kappa \}`$, a set of non-zero divisors, the image of each $`d_\alpha `$ under the isomorphism must be an open dense subset $`O_\alpha `$ of $`X`$. So there is a point $`x_{\alpha <\kappa }O_\alpha `$. Let $`O=Xc\{x\}`$. Then $`O`$ corresponds in $``$ to an irreducible $`p`$, and since $`x`$ is in each $`O_\alpha `$ we can’t have $`d_\alpha p`$ for any $`\alpha `$.

(4) implies (3). Given a ccc core compact space $`(X,𝒯)`$, by Theorem 1.1 in  $`(𝒯,)`$ is an upwards-continuous lattice. It is clear that it is also downwards-ccc, distributive, and complete. Given a collection $`\{O_\alpha |\alpha <\kappa \}`$ of open dense subsets of $`X`$, we see as above that each $`O_\alpha `$ is a non-zero divisor. So there is an irreducible $`O`$ in the lattice which does not contain any $`O_\alpha `$. $`C=XO`$ is an irreducible set for $`(X,𝒯)`$, and $`C`$ meets every $`O_\alpha `$.

(3) implies (1). Given an upwards-ccc partially ordered set $`(𝒫,)`$ consider the partial order topology $`𝒯`$ generated by all sets

$$x=\{y𝒫|xy\}$$

for $`x𝒫`$. $`(𝒫,𝒯)`$ is a ccc space. It is locally compact, hence core compact. Given a collection $`\{C_\alpha |\alpha <\kappa \}`$ of cofinal subsets of $`𝒫`$, each $`C_\alpha `$ is dense, so there is an irreducible set $`C`$ meeting each $`C_\alpha `$. But in the partial order topology a set is irreducible only if it is directed, and so we are done. ∎

We note that each of properties (2) or (3) in Proposition 1 also implies the Baire category formulation of $`\kappa <𝔪`$, to wit, no ccc compact Hausdorff space is the union of $`\kappa `$ or fewer nowhere dense sets. We gave in the proof an argument that (3) implies (1) so that we may now suggest that (3) is a simultaneous generalization of the partial order and Baire category forms of $`\kappa <𝔪`$.

We also wish to point out that the equivalence of properties (2) and (3) may be established directly using the sobrification of the space $`(X,𝒯)`$. The sobrification of $`(X,𝒯)`$ is the collection of irreducible subsets of $`X`$, topologized by the lower Vietoris topology, and $`(X,𝒯)`$ can be mapped continuously to its sobrification by sending each $`x`$ to $`c\{x\}`$ (and thus the image of $`(X,𝒯)`$ is the $`T_0`$-identification of $`(X,𝒯)`$). A space is core compact if and only if its sobrification is locally compact .

In view of the Baire category form of $`\kappa <𝔪`$ mentioned above, the core compactness in property (3) may be replaced by any stronger property that is weaker than compact Hausdorff. Likewise, locally compact sober in property (2) may be replaced by anything stronger which is implied by compact Hausdorff, including locally compact supersober. Sober may also be weakened to quasisober, which is sober without the $`T_0`$ assumption, since that is sufficient for the Hofmann-Mislove Theorem.

###### Corollary 1.

The following are equivalent.

1. Martin’s Axiom.
2. In a ccc locally compact sober space the intersection of fewer than $`𝕔`$ open dense sets is non-empty.
3. In a ccc core compact space if $`𝒟`$ is a collection of fewer than $`𝔠`$ open dense sets then there is an irreducible set meeting every member of $`𝒟`$.
4. Let $``$ be a downwards-ccc upwards-continuous distributive complete lattice. If $`D`$ is a set of fewer than $`𝔠`$ non-zero divisors of $``$ then there is an irreducible $`p`$ such that no element of $`D`$ is below $`p`$.

There is another cardinal, $`𝔭`$, for which we can get some similar results.

###### Definition 6.

In a partially ordered set $`(𝒫,)`$ a set $`R`$ is upwards-centered if every finite subset of $`R`$ has an upper bound in $`𝒫`$. $`(𝒫,)`$ is upwards-$`\sigma `$-centered if $`𝒫`$ is the union of countably many centered subsets. The cardinal $`𝔭`$ is the least cardinal such that there exists a non-empty upwards-$`\sigma `$-centered partially ordered set $`(𝒫,)`$ and a collection $`\{D_\alpha |\alpha <𝔭\}`$ of cofinal subsets of $`𝒫`$ such that no upwards-directed subset of $`𝒫`$ meets each $`D_\alpha `$.111This is not the original definition. if $`𝒜`$ is a collection of infinite sets we say that an infinite set $`B`$ is a pseudo-intersection of $`𝒜`$ if $`BA`$ is finite for every $`A𝒜`$. Then $`𝔭`$ is the least cardinal such that there exists a collection $`𝒜`$ of subsets of $`\omega `$, with the cardinaltiy of $`𝒜`$ equal to $`𝔭`$, and although every finite subset of $`𝒜`$ has an infinite intersection, $`𝒜`$ has no infinite pseudo-intersection. It was Murray Bell  who proved that this definition is equivalent to the $`\sigma `$-centered partial order definition given above, and Fremlin  refers on page 25 to this result as “Bell’s Theorem.” A space $`(X,𝒯)`$ is $`\sigma `$-centered if $`𝒯\{\varphi \}`$ with reverse inclusion is upwards-$`\sigma `$-centered. A complete lattice $``$ is downwards-$`\sigma `$-centered if $`\{0\}`$ with order reversed is upwards-$`\sigma `$-centered.

The Baire category formulation of $`\kappa <𝔭`$ is that no separable compact Hausdorff space is the union of fewer than $`𝔭`$ nowhere dense sets. It can be shown that $`𝔪𝔭𝔠`$.

###### Proposition 2.

Let $`\kappa `$ be a cardinal. The following are equivalent.

1. $`\kappa <𝔭`$
2. In a separable locally compact sober space the intersection of $`\kappa `$ or fewer open dense sets is non-empty.
3. In a separable core compact space if $`𝒟`$ is a collection of $`\kappa `$ or fewer open dense sets then there is an irreducible set meeting every member of $`𝒟`$.
4. In a $`\sigma `$-centered locally compact sober space the intersection of $`\kappa `$ or fewer open dense sets is non-empty.
5. In a $`\sigma `$-centered core compact space if $`𝒟`$ is a collection of $`\kappa `$ or fewer open dense sets then there is an irreducible set meeting every member of $`𝒟`$.
6. Let $``$ be a downwards-$`\sigma `$-centered upwards-continuous distributive complete lattice. If $`D`$ is a set of $`\kappa `$ or fewer non-zero divisors of $``$ then there is an irreducible $`p`$ such that no element of $`D`$ is below $`p`$.

###### Proof.

The equivalence of properties (1), (4), (5), and (6) may be established as in the proof of Proposition 1. Properties (2) and (3) are implied by (4) and (5), respectively, and each of them implies the Baire category form of property (1). ∎

We do not know if $`\sigma `$-centered is equivalent to separable for locally compact sober spaces, but this is true for locally compact supersober spaces.


text:
# Deterministic diffusion in flower shape billiards

## I Introduction

One of the central themes in the theory of nonequilibrium statistical mechanics is to assess the importance of deterministic chaos for understanding transport processes such as diffusion . Simple model systems appear to be most suited for studying the detailed relation between microscopic chaos and macroscopic transport. Along this line of research, the parameter dependent diffusion coefficients of strongly chaotic dynamical systems have been investigated for one- and two-dimensional mappings , periodic Lorentz gases , and billiards in an external field . That diffusion coefficients can be fractal functions of control parameters was first observed in a simple one-dimensional mapping generalizing a random walk on the line . The origin of this fractality may be attributed to the topological instability of orbits under parameter variation that affects the parameter dependence of the diffusion coefficient in a nontrivial way. Based on the analysis of such simple systems, it was conjectured that fractal diffusion coefficients are rather generic for low-dimensional fully chaotic dynamical systems exhibiting some spatial periodicity . Indeed, recently it was found that in case of billiards in an external field the diffusion coefficient again exhibits a highly irregular structure .

The standard periodic Lorentz gas is one of the typical models for studying deterministic normal diffusion (see, e.g., Refs.  and further references therein). That it is strongly chaotic and exhibits normal diffusion was proven by Bunimovich and Sinai . Machta and Zwanzig have calculated the diffusion coefficient of this model from computer simulations at some parameter values, and they have matched their results to a simple analytical random walk approximation . That the diffusion coefficient in the standard periodic Lorentz gas is indeed a non-trivial function of the parameter was first reported in Ref. . Here the analysis by Machta and Zwanzig was refined by suggesting two methods for systematically correcting their random walk approximation. However, whether the numerically detected irregularities in the diffusion coefficient were of a fractal nature remains an open question. More recently, a third approximation scheme was proposed by deriving a Green-Kubo formula that exactly generalizes the Machta-Zwanzig approximation . Applying all these methods led to the conclusion that including long-term correlations, or memory effects, was indispensable to reproduce the precise functional form of the parameter dependent diffusion coefficient for the standard periodic Lorentz gas.

One of the essential problems in the analysis of diffusion in the standard periodic Lorentz gas is that the parameter range of normal diffusion is very limited. In this small region, the irregular behavior of the parameter dependent diffusion coefficient shows up on very fine scales and appears to be rather smooth within the range of precision available from computer simulations . Consequently, the question about the existence of a fractal diffusion coefficient is very difficult to answer for this model. As the main reason for this behavior it might be suspected that the topological instability of the standard periodic Lorentz gas under parameter variation is not strong enough to generate more pronounced irregularities in this region. The main purpose of this paper is therefore to propose a billiard without external field which is very similar to the standard periodic Lorentz gas, but which has a geometry, and an associated range of control parameters exhibiting normal diffusion, with stronger topological instabilities. This way, we intend to learn more about the emergence of possible fractal structures for diffusion coefficients in billiards. As we will show, our model indeed generates a considerably stronger irregular parameter dependence of the diffusion coefficient than in the standard Lorentz gas. By applying the set of approximation methods mentioned above we argue that long-range dynamical correlations, or memory effects of orbits, are again at the origin of this irregularity, as in case of simple one- and two-dimensional maps.

Our paper is composed of seven sections: In section II, we introduce the flower shape billiard. Numerical results depicting the non-trivial parameter dependence of the diffusion coefficient are shown in section III. In section IV,V, and VI, we briefly review the different approaches to understand the parameter dependence of diffusion coefficients in deterministic dynamical systems, i.e., the Machta-Zwanzig approximation, Klages-Dellago correction methods, as well as the approach based on a suitable Green-Kubo formula for diffusion, and we apply them to the flower shape billiard. Summary and conclusions are contained in section VII.

## II The flower shape billiard

The two-dimensional class of billiards we consider here consists of a point particle of mass $`m`$ moving in a plane such that its Hamiltonian is

$$H=\frac{1}{2m}p_x^2+\frac{1}{2m}p_y^2,$$
(1)
where $`x`$ and $`y`$ denote the Cartesian coordinates of the position in the plane while $`p_x`$ and $`p_y`$ are the corresponding momenta. The point particle undergoes elastic collisions with obstacles that are fixed in the plane. All the obstacles have the same shape, and their centers are situated on a triangular lattice according to

$$𝒒_c=m_c\mathbf{}_\mathrm{𝟏}+n_c\mathbf{}_\mathrm{𝟐},$$
(2)
as defined in terms of the fundamental translation vectors of the triangular lattice

$$\mathbf{}_\mathrm{𝟏}=(0,1)$$
(3)
and

$$\mathbf{}_\mathrm{𝟐}=(\frac{\sqrt{3}}{2},\frac{1}{2}),$$
(4)
where $`m_c`$ and $`n_c`$ are integers.

If all the pairs of integers are selected, we fill the whole triangular lattice with hard wall obstacles, and the billiard is invariant under the group of spatial translations generated by the vectors Eq. (2). Accordingly, the whole lattice can be mapped onto a so-called Wigner-Seitz cell with periodic boundary conditions. The elementary Wigner-Seitz cell of the triangular lattice is a hexagon of area

$$A_{WS}=|\mathbf{}_\mathrm{𝟏}\times \mathbf{}_\mathrm{𝟐}|=\frac{\sqrt{3}}{2}.$$
(5)
In this paper, we propose an open billiard consisting of flower shape obstacles instead of disks, which belongs to the general class of periodic Lorentz gases whose normal diffusion has been proven by Bunimovich and Sinai . The mixing property and the extension of such billiards to higher-dimensional gases have been studied by Chernov . As shown in Fig. 1, the space between the obstacles forms the two-dimensional domain of the billiard where the point particle moves freely and collides with the obstacles obeying the law of elastic reflection.

A single scatterer of our billiard is defined as follows: First, we consider the inner hexagon whose vertices are on the middle points of the sides of the hexagon of the elementary Wigner-Seitz cell, as depicted by the dotted lines in Fig. 2. Next, we join six arcs which have the same radii and touch the inner hexagon. Then we obtain the flower-shaped obstacle shown in Fig. 2. Note that the radius $`r`$ of one arc which consists in a petal of the flower-shaped obstacle can be changed from $`1/(4\sqrt{3})`$ to infinity. According to this construction, the position space forms a two-dimensional torus. The motion of the point particle in the infinite lattice is unbounded so that transport by diffusion is a priori possible. Indeed, we will show that the diffusion of point particles in the billiard of the flower-shaped obstacles is normal. When the dynamics is reduced to the Wigner-Seitz cell, the position of the particle inside this cell must be supplemented by a lattice vector of the type of Eq. (2) in order to determine the actual position of the particle in the infinite lattice. This lattice vector changes in discrete steps at each crossing of the border of the elementary Wigner-Seitz cell.

A billiard whose obstacles are disks, or, in higher dimensions, spheres, is called a periodic Lorentz gas, and this model serves as a typical example for studying deterministic diffusion . The diffusion coefficient of this standard periodic Lorentz gas has been studied in various ways both analytically and numerically, where recent work focused particularly onto its density dependence (see Refs.  and further references therein). However, the density of this model cannot be varied much because of the condition of finite horizon, which prohibits collision-free ballistic motion and thus guarantees the existence of normal diffusion . Consequently, the diffusion coefficient exists in a very limited range of parameters only, and whether the diffusion coefficient of the standard periodic Lorentz gas is a fractal function of the density of scatterers appears to be an open question.

Let us introduce the Liouville equilibrium invariant measure given by

$$d\mu _e=I(x,y)\delta (HE)dxdydp_xdp_y,$$
(6)
where $`I(x,y)`$ is the indicator function of the billiard domain, and $`E`$ is the energy of the point particle. Averages over this invariant measure are denoted by $`<>`$. This measure is normalizable for the reduced dynamics in an elementary Wigner-Seitz cell where the area of the billiard domain takes a finite value. In this finite case, the Liouville invariant measure is a probability measure which defines the microcanonical ensemble of equilibrium statistical mechanics. The flower shape billiard belongs to the class of dispersing billiards whose hyperbolicity has been proven by Sinai . Consequently, it is known that the motion of the point particle in the elementary Wigner-Seitz cell of our billiard is hyperbolic in the sense that all orbits are unstable of saddle type with nonvanishing Lyapunov exponents, and time averages are equal to averages over the Liouville equilibrium invariant measure.

## III  Curvature dependence of the diffusion coefficient

Since the system of flower-shaped obstacles is fully chaotic, and by working in the regime of finite horizon, we may expect that diffusion is normal in the sense that the position is asymptotically a Gaussian random variable with a variance growing linearly in time. Consequently, the diffusion coefficient exists and is finite . Indeed, we checked numerically that the variance is proportional to time after sufficiently long time evolution.

The diffusion coefficient $`D`$ is given by the Einstein formula,

$$D=\underset{t\mathrm{}}{lim}\frac{1}{4t}<\{𝒒(t)𝒒(0)\}^2>,$$
(7)
and according to this formula the diffusion coefficient was calculated from computer simulations in the flower shape billiard where the curvature $`\kappa `$ of the petals is varied from 0 to its maximum $`4\sqrt{3}`$. The results are depicted in Fig. 3. In this figure, we observe a non-trivial structure depending on the curvature $`\kappa `$ of the arc defining the petal of the flower-shaped obstacles.

The gross features of the curvature dependence for the diffusion coefficient can qualitatively be explained as follows: When the curvature of the petal of the flower-shaped obstacle is zero, the inner hexagon shown by the dotted lines in Fig. 2 connects to the six hexagons surrounding it. In this case, the point particle remains forever localized in compact domains bounded by the three neighbouring hexagons. For this specific value of the control parameter, the motion of the point particle is completely predictable because the compact domain is an equilateral triangle, and the system is integrable.

When the curvature becomes positive, the point particle can run away from the compact domain, and diffusion occurs. As already explained, at all positive curvatures of the petal, even if they are very small, the motion of the point particle is fully chaotic and the horizon is finite, hence diffusion is expected to be normal. The diffusion coefficient starts to increase from zero according to the linear increase of the curvature of the petal, and related to the fact that the space between petals also increases.

When the radius of the petal is equal to $`R_L=\sqrt{3}/40.433`$, which is the distance between the center of the hexagon and the tangent point to the hexagon, the obstacle becomes a disk of radius $`R_L`$, that is, for this parameter value our billiard is precisely the same as the conventional periodic Lorentz gas. This point corresponds to the curvature $`\kappa 2.309`$ in Fig. 2.

When the radius $`r`$ of the curvature of the petal decreases below $`R_L`$, the point particle is much more likely to be trapped in the space between two obstacles. This appears to be due to the formation of wedges between any two petals of a flower shape obstacle.

The inset of Fig. 3 depicts a zoom on the curve showing the fine structure on smaller scales with respect to curvature. We remark that the apparently continuous fluctuations therein are within the numerical errors, that is, we confirmed the convergence of our results within a precision of order $`10^4`$ by taking an average over $`10^{10}`$ initial conditions. Unfortunately, with our computational power it is impossible to check whether this oscillatory behavior persists on even finer scales.

## IV  Machta-Zwanzig approximation for diffusion coefficients

In Ref. , Machta and Zwanzig have obtained a simple analytical approximation for the diffusion coefficient of the periodic Lorentz gas which yields asymptotically correct results in the limit of small gaps between disks. In this case, the particle is for a long time somewhat trapped in the triangular regions between three adjacent scatterers. Hence, the particle is supposed to loose the memory of its past itinerary due to the multiple scattering in the trap region, and the transition probabilities to the neighbouring triangular cells are assumed to be equivalent. As was shown in Ref. , the average rate $`\tau ^1`$ of such transitions can be calculated from the fraction of phase space available for leaving the trap divided by the total phase space volume of the trap leading to

$$\tau =\pi A/(3W),$$
(8)
where $`A`$ is the area of the trap and $`W`$ the width of the gap between the disks.

The flower shape billiard has similar types of traps as the periodic Lorentz gas. Accordingly, the Machta-Zwanzig approximation can be applied to the flower shape billiard as well, and Eq. (8) holds again for the average trapping time. Hence, we only need to calculate the area of the trap and the gap between the petals from simple geometrical considerations yielding

$$A=\frac{3\sqrt{3}}{4}3h\left[\sqrt{3}h+\sqrt{r^2h^2}\right]$$
(9)
and

$$W=\frac{1}{2}\left[\sqrt{3}h+\sqrt{r^2h^2}\right],$$
(10)
where

$$h=\frac{1}{2}\left(\frac{\sqrt{3}}{4}r\right).$$
(11)
In the above, $`r`$ denotes the radius of the curvature of the petal.

The distance $`l`$ between the centers of the flower shape obstacles is $`1/\sqrt{3}`$. Assuming that the gap size $`W`$ is very narrow leads to the Machta-Zwanzig random walk approximation for the diffusion coefficient

$$D_{MZ}=\frac{l^2}{4\tau }$$
(12)
with $`\tau `$ being given by Eq. (8) and supplemented by Eqs. (9)-(11) for the flower shape case. As is shown in Fig. 4, the Machta-Zwanzig approximation works very well in the vicinity of zero curvature of the petal only.

## V  Klages-Dellago corrections of the Machta-Zwanzig approximation

In Ref. , Klages and Dellago have generalized the Machta-Zwanzig approximation for the standard periodic Lorentz gas by taking memory effects of orbits into account. Their generalization is based on the observation that, except in the asymptotic limit of narrow gap sizes, the diffusive dynamics is not a simple Markov process, in the sense that there exist non-vanishing dynamical correlations. By mapping the orbit of a particle onto a suitable symbolic dynamics they numerically calculated the probabilities to obtain certain symbol sequences of finite length. Increasing the length of these symbol sequences yielded systematic corrections of the Machta-Zwanzig approximation. In Ref. , two schemes directly emerging from this approach were discussed, one suggesting simple heuristic corrections to the simple random walk model of diffusion Eq. (7), and another one employing lattice gas computer simulations defined by these probabilities. In this section, we apply these two methods to the flower shape billiard in order to systematically correct the Machta-Zwanzig approximation. A third scheme starting from a Green-Kubo formula for diffusion will be discussed in section VI.

The Machta-Zwanzig approximation assumes that a particle jumps from one trap to a neighbouring trap situated on the hexagonal lattice of traps. However, there exist non-vanishing probabilities that a particle can jump to next nearest neighbours, or even farther, without collisions. Accordingly, we should correct the Machta-Zwanzig approximation for the flower shape billiard by using the probabilities $`p_{cf1}`$ and $`p_{cf2}`$ of those collisionless flights which lead from one cell directly to its second nearest neighbours, or to the third nearest neighbours, respectively. The distances $`l_1`$ and $`l_2`$ between the centers of a trap to respective second and third neighbours are

$$l_1=\sqrt{3}l,l_2=\sqrt{7}l.$$
(13)
The diffusion coefficient $`D_{cf}`$ with corrections due to these collisionless flights then reads

$`D_{cf}`$ $`=`$ $`(1p_{cf1}p_{cf2})D_{MZ}+p_{cf1}{\displaystyle \frac{l_1^2}{4\tau }}+p_{cf2}{\displaystyle \frac{l_2^2}{4\tau }}`$ (14)
$`=`$ $`(1+2p_{cf1}+6p_{cf2})D_{MZ}.`$ (15)
Next we take memory effects of orbits due to backscattering into account. For this purpose, orbits are coded by labeling the entrance through which a particle enters a trap with $`z`$, the exit to the left of this entrance with $`l`$, and the one to the right with $`r`$. Thus, an orbit can be mapped onto a sequence of symbols $`z,l`$, and $`r`$. For example, $`p(z)`$ is the backscattering probability $`p_{bs}`$, which is the probability of the moving particle to leave the trap through the same gate where it entered. The Machta-Zwanzig approximation assumes that $`p(z)=p(l)=p(r)=1/3`$. However, in general $`p(z)`$ is not close to $`1/3`$ as shown in Fig. 5, because the actual orbits do not loose their memory during their itineraries.

A more profound explanation for the complicated functional form of $`p(z)`$ may be provided in terms of the theory of chaotic scattering: Chaotic scattering systems with multiple exit modes typically have fractal phase space boundaries separating the sets of initial conditions (basins) going to the various exits. However, open systems such as a three-disk scatterer of the periodic Lorentz gas possess the even stronger property of being Wada, that is, any initial condition which is on the boundary of one exit basin is also simultaneously on the boundary of all the other exit basins . Changing the curvature $`\kappa `$ sensitively affects the highly irregular structure of these basin boundaries. Consequently, Fig. 5 may be understood as reflecting the topological instability of Wada basins under parameter variation, and as we will now show this is reflected in the parameter dependence of the diffusion coefficient.

Modifying the Machta-Zwanzig random walk by including the backscattering probability $`p(z)`$ we obtain the diffusion coefficient

$$D_{BS}=\frac{(1p(z))l_2^2}{4(2\tau )}=(1p(z))\frac{3}{2}D_{MZ}.$$
(16)
Combining the effects of collisionless flights and backscattering yields as a first order approximation

$$D_1=\frac{3}{2}(1p(z))(1+2p_{cf1}+6p_{cf2})D_{MZ}.$$
(17)
Higher-order approximations of the diffusion coefficient, as related to longer symbol sequences and respective probabilities such as $`p(lrz\mathrm{})`$, can be derived in the same way . For the flower shape billiard, respective results are shown in Fig. 6.

The above correction methods assume that all orbits follow higher-order Markov processes, where correlations are present in form of initial transient times before the variance becomes linear in time. This dynamics appears to be more suitably represented in form of lattice gas simulations on a honeycomb lattice, where the sites of the lattice represent the traps. Indeed, for the periodic Lorentz gas such lattice gas simulations were performed in Ref.  confirming fast convergence to the numerically exact results. Compared to that scheme, the convergence of the intuitive correction method described above is firstly slower, and secondly it is not everywhere converging to the numerically exact results, which is due to the fact that this approach was purely of a heuristic nature.

We also performed lattice gas simulation in case of the flower shape billiard according to the following prescription: Particles hop from site to site with frequency $`\tau ^1`$, which is identical to the hopping frequency used in the Machta-Zwanzig approximation. The hopping probabilities are given by the backscattering probability $`p(z)`$ and by those corresponding to respective longer symbol sequences. The diffusion coefficient is then obtained from the Einstein formula Eq. (7) in the limit when the variance is getting proportional to time. The correlations in the actual orbits are thus systematically and exactly filtered out according to the length of the symbol sequences.

In Fig. 7, the results of such higher-order approximations according to lattice gas simulations are shown. One can see that the convergence to the numerically exact results is not only much better than in Fig. 6, but even exact. Strong memory effects are clearly visible especially after the diffusion coefficient curve takes its maximum. In the previous heuristic modifications to the simple random walk model, the dynamics was only modeled for a limited number of time steps as a Markov process. Fig. 6 suggests that correlations as contained in the symbol sequences are more suitably represented by higher-order iterations in form of lattice gas simulations. However, a disadvantage is that the lattice gas scheme requires a second round of computations which is put on top of the previous simulations, by again looking at the time evolution of an initial ensemble of points.

## VI  The Green-Kubo formula approach

The main drawbacks of the two methods described above were, firstly, that the heuristic corrections of the Einstein formula were not converging exactly to the numerically exact results, and, secondly, that the lattice gas simulations were merely a numerical scheme without being represented in form of analytical approximations. These deficiencies were essentially resolved in Ref.  by the derivation of a Green-Kubo formula which employs the symbolic dynamics on the hexagonal lattice of traps introduced in section V. The result reads

$$D=\frac{1}{4\tau }C_0+\frac{1}{2\tau }\underset{n=1}{\overset{\mathrm{}}{}}C_n,$$
(18)
with

$$C_n:=𝐣(𝐱_0)𝐣(𝐱_n)$$
(19)
being the velocity autocorrelation function related to jumps $`𝐣(𝐱_n)`$ on the hexagonal lattice at time step $`n`$. These jumps are suitably defined in terms of the lattice vectors Eqs. (3),(4). That is, any symbol sequence of an orbit on the hexagonal lattice of traps defines a respective chain of lattice vectors. The averages indicated by the brackets in Eq. (19) are calculated by weighting the respective scalar products of lattice vectors with the corresponding conditional probablities $`p(\alpha \beta \gamma \mathrm{}),\alpha ,\beta ,\gamma \{l,r,z\}`$. In Eq. (19), $`\tau `$ is the mean time of free flight between symbol changes, and it is given by Eq. (8). Eq. (18) is thus the honeycomb lattice analogue to the Green-Kubo formula derived by Gaspard for the Poincaré-Birkhoff map of the periodic Lorentz gas .

It is easy to see that the first term in Eq. (18) yields the Machta-Zwanzig approximation Eq. (12). Higher-order corrections can then be calculated by defining the hierarchy of approximations

$$D_n=\frac{l^2}{4\tau }+\frac{1}{2\tau }\underset{\alpha \beta \gamma \mathrm{}}{}p(\alpha \beta \gamma \mathrm{})\mathbf{}\mathbf{}(\alpha \beta \gamma \mathrm{})$$
(20)
with $`n>0`$ being the number of symbols and $`D_0(w)`$ given by Eq. (12), where, again, $`\mathbf{}(\alpha \beta \gamma \mathrm{})`$ are suitable lattice vectors.

The impact of dynamical correlations on the diffusion coefficient can now be understood by analyzing the single contributions in terms of the correlation function $`C_n`$ as contained in the Green-Kubo formula Eq. (18). In fully chaotic systems such as the periodic Lorentz gas and the flower shape billiard, the velocity correlation function decays exponentially, which is in agreement to the results depicted in Fig. 8. By comparing this figure to Fig. 9 one can learn how the irregularities of the correlation function determine the parameter dependent diffusion coefficient: Let us start with the first-order approximation of Eq. (20) which reads $`D_1=D_0+D_0(13p(z))`$. The functional form of $`p(z)`$ in Fig. 5 thus qualitatively explains the position of the global maximum of the diffusion coefficient, because at this value of the curvature the probability of backscattering is minimal. Adding up the three-jump contributions coming from $`C_2`$ furthermore yields the most important quantitative contributions in this region of the curvature. In the region of large curvature the diffusion coefficient decays monotonically according to the effect of two-hop correlations covered by $`C_1`$. However, note the large fluctuations of the correlation function $`C_n`$ as well as of the diffusion coefficient approximations $`D_n`$ in this regime both indicating the dominant effect of long-range higher-order correlations. Studying the detailed convergence of the approximations depicted in Fig. 9 shows that correlations due to orbits with longer symbol sequences yield irregularities in the parameter dependence of the diffusion coefficients on finer and finer scales.

## VII  Summary and conclusion

In this paper, we have introduced a novel variant of the periodic Lorentz gas by assigning a flower-shaped geometry to the scatterers. Although both systems are rather similar in the sense that they are both fully chaotic and exhibit normal diffusion in a certain parameter range, we have found that the diffusion coefficient in the flower shape geometry is considerably more irregular under parameter variation than the one obtained from circular disks as scatterers. We have analyzed these irregularities by three different methods which all start from correcting the Machta-Zwanzig random walk approximation for the diffusion coefficient. All these improved approximation schemes use a symbolic dynamics which maps the orbits of moving particles to symbol sequences according to traps situated on a hexagonal lattice. We have discussed the convergence of these different approximation schemes, and we have shown how they enable a detailed understanding of the precise shape of the parameter dependent diffusion coefficient in the flower shape billiard in terms of long-range dynamical correlations.

The Green-Kubo formula introduced in Ref.  appears to be most suitable for understanding the irregular behavior of the parameter dependent diffusion coefficient, because it conveniently transforms the diffusive dynamics into a sum over the velocity correlation function, whose specific parameter dependence can in turn be analyzed step by step. In particular, this approach yields an exact convergence to the parameter dependent diffusion coefficient as obtained from simulations.

Interestingly, when the correlation function decays in time, the frequency of oscillations as a function of the control parameter increases. The relation between this decay in time and the increase of the frequency of these oscillations determines the strength of the irregularities on fine scales of the resulting parameter dependent diffusion coefficient. The question of the existence of fractal diffusion coefficients in billiards such as periodic Lorentz gases with circular or flower-shaped scatterers might thus be answered by using Green-Kubo formulas if the respective correlation functions could be evaluated more precisely for large enough times. Indeed, in Ref.  the highly irregular diffusion coefficient of an open billiard in an external field has already been investigated along these lines by relating the Poincaré-Birkhoff version of the Green-Kubo formula to fractal Weierstrass functions. The joint efforts compiled in Refs.  may therefore be considered as first steps towards answering the conjecture of Refs. , which suggested a possible universality of fractal diffusion coefficients in low-dimensional fully chaotic dynamical systems exhibiting some spatial periodicity, for the case of chaotic Hamiltonian dynamical systems such as particle billiards.

###### Acknowledgements.

T. H. thanks Dr. B. Komiyama for support and encouragement in this research. R. K. thanks Mr. J. Aguirre for pointing out Refs.  to him. P. G. thanks the National Fund for Scientific Research (F.N.R.S. Belgium) for financial support.


text:
# Inclusive Spectra and Quantum Stochastic Processes

## Abstract

We construct a simple, exactly solvable model illustrating the conjecture that the evolution of the wave function of a system of very many interacting particles is likely to drive the one-particle spectrum towards a ”thermal” shape. We speculate that this pseudo-thermalization, operating at the level of multiparton wave functions could perhaps explain the apparently thermal spectra observed in those high-energy processes, where the conditions for a genuine thermalization are not met.

PACS numbers: 02.50Ey, 13.85Hd

Key words: stochastic, hadron, inclusive

LPT Orsay 02/31 (April 2002)

Consider a multiparticle state in a finite box. If the state were belonging to a classical statistical ensemble, it would thermalize in the presence of an interaction between particles. In particular, the one-particle momentum distribution would tend to the Boltzmann form with time going on. Suppose now that the state is a pure quantum state. We wish to address the following questions: is it likely that the one-particle spectrum is approaching, for time $`+\mathrm{}`$, a limit independent of the initial state of the system, like in a generic Markov process? Is it plausible that this limit distribution looks thermal, reflecting the features of the corresponding classical ensemble? Although these questions may sound academic, they are in fact motivated by phenomenology. However, since the answers are, we believe, of intrinsic theoretical interest, we postpone the discussion of the phenomenological matters to the final paragraphs and focus first on what appears to be, hopefully, more than an amusing speculation.

It is difficult to approach the problem in its full generality. Consequently, we adopt the one-particle approximation: we follow a single test particle moving in an external mean field, whose action represents the multiparticle environement. We know, that in a classical gas this motion generically leads to thermalization. A time arrow is attached to the motion of the test particle. This time arrow is also present when one goes over to the quantum description. This means that the evolution generator describing the motion of the test particle is non-hermitian. Hence, it will, in general, have complex eigenvalues. Consequently, a single quantum state of the test particle, independent of its initial state, will be selected for time $`+\mathrm{}`$. This state, to be regarded as the stationary one, is reached by a process, which corresponds to the decay of a quantum fluctuation and which is likely to resemble thermalization. That is the general idea of this paper. It will be illustrated by a simple example.

Consider a gas of classical particles first. How does the thermalization occur? A given particle of the gas is slowed down (accelerated) when it collides with a particle less (more) energetic than itself. Eventually, its momentum fluctuates near the average momentum measured by the gas temperature. To mock this behavior assume that our test particle gets random kicks, such that the transition $`pp\pm \delta `$ occurs with probability

$$Prob(pp\pm \delta )e^{\mathrm{\Delta }^\pm (p)}$$
(1)
where

$$\mathrm{\Delta }^\pm (p)=\frac{1}{2T}[E(p)E(p\pm \delta )]$$
(2)
Here $`E(p)`$ is the energy of a particle. The transition probability (1) satisfies the detailed balance equation

$$Prob(p)Prob(pp\pm \delta )=Prob(p\pm \delta )Prob(p\pm \delta p)$$
(3)
with $`Prob(p)\mathrm{exp}[E(p)/T]`$, which is therefore the limit distribution of the process defined by (1) foot3 . Equation (1) violates the time reversal. It must be so, the process has the time arrow built into it, since it mocks the effect of the interaction of the test particle with its multiparticle environement.

In order to simplify writing we assume first that $`\delta `$ is a constant vector. Thus the motion is one-dimensional. The extension of the discussion to a spectrum of momentum transfers is trivial.

Eq. (1) implies the following master equation for the momentum probability distribution $`w(p,t)`$:

$`w(p,t)/t=w(p+\delta ,t)e^{\mathrm{\Delta }^{}(p+\delta )}`$
$`+w(p\delta ,t)e^{\mathrm{\Delta }^+(p\delta )}w(p,t)[e^{\mathrm{\Delta }^{}(p)}+e^{\mathrm{\Delta }^+(p)}]`$ (4)
where the time $`t`$ is measured in appropriate units, so that the interaction strength does not appear explicitly. It can be immediately seen that $`w(p)=e^{E(p)/T}`$ is a stationary solution of (4), as expected.

We now go over to the quantum world and replace the particle by a wave. We assume that the kicks are due to the interaction

$$V=\underset{p}{}[\eta ^{}e^{\mathrm{\Delta }^{}(p+\delta )}a_p^{}a_{p+\delta }+\eta ^+e^{\mathrm{\Delta }^+(p\delta )}a_p^{}a_{p\delta }]$$
(5)
where $`\eta ^\pm `$ are some phase factors, scalar functions of $`p`$ and $`\delta `$, while $`a_p^{}`$ ($`a_p`$) are time independent creation (annihilation) operators. The interaction operator is not hermitian, the time reversal invariance is broken like in the classical model.

The number of particles is - by assumption - conserved foot4 . Hence the state at time $`t`$ has necessarily the form

$$t=\underset{p}{}\psi (p,t)a_p^{}0$$
(6)
The norm of this state is time dependent, the evolution being non-unitary. This does not matter for our purposes, since we are interested in this paper in the shape of the momentum distribution only. The Schroedinger equation describing the evolution of the wave packet (6) reads

$`i\psi (p,t)/t=\eta ^{}e^{\mathrm{\Delta }^{}(p+\delta )}\psi (p+\delta ,t)`$
$`+\eta ^+e^{\mathrm{\Delta }^+(p\delta )}\psi (p\delta ,t)`$ (7)
In order to determine the behavior of the solution we have to solve first the eigenvalue equation:

$$\eta ^{}e^{\mathrm{\Delta }^{}(p+\delta )}\psi (p+\delta )+\eta ^+e^{\mathrm{\Delta }^+(p\delta )}\psi (p\delta )=\lambda \psi (p)$$
(8)
The transition probabilities are determined by the values of $`\mathrm{\Delta }^\pm `$. In the absence of a better insight into the problem it is reasonable to assume that $`\eta ^\pm `$ also depends primarily on this parameter. In order to obtain an analytically solvable model we postulate a linear behavior

$$\text{Arg}(\eta ^\pm )=\xi _0+\xi _1\mathrm{\Delta }^\pm (p\delta )$$
(9)
Set $`\psi (p)=e^{E(p)(1+2i\xi _1)/2T}z(p)`$ to get

$$e^{i\xi _0}[z(p+\delta )+z(p\delta )]=\lambda z(p)$$
(10)
From here on the calculation becomes elementary. Assuming periodic boundary conditions one has

$$2e^{i\xi _0}\mathrm{cos}(x\delta )\stackrel{~}{z}(x)=\lambda \stackrel{~}{z}(x)$$
(11)
where

$$\stackrel{~}{z}(x)=\underset{p}{}z(p)e^{ipx}$$
(12)
Hence the sought eigenvalues and eigenvectors are

$$\lambda _x=2e^{i\xi _0}\mathrm{cos}(x\delta )$$
(13)
and

$`\psi _x^{(s)}(p)=e^{E(p)(1+2i\xi _1)/2T}\mathrm{cos}(px)`$ (14)
$`\psi _x^{(a)}(p)=e^{E(p)(1+2i\xi _1)/2T}\mathrm{sin}(px)`$ (15)
When one considers a spectrum of momentum transfers, one gets instead of (13):

$$\lambda _x=e^{i\xi _0}\underset{\delta }{}\mu (\delta )\mathrm{cos}(x\delta )$$
(16)
with some $`\mu (\delta )`$, which we assume to be a positive weight so that, generically, $`x=0`$ corresponds to the eigenvalue with the largest modulus. The general solution of the Schroedinger equation is

$$\psi (p,t)=\underset{x}{}[A_x\psi _x^{(a)}(p)+S_x\psi _x^{(s)}(p)]e^{it\lambda _x}$$
(17)
Let $`0<\xi _0<\pi `$. Then only the contribution corresponding to $`x=0`$ survives in the limit $`t+\mathrm{}`$, because this is the one whose norm is ”explosing” the most rapidly (remember that the state is not normalized; the $`x=0`$ eigenstate should be regarded as the stationary solution). Consequently, the asymptotic spectrum is $`w(p)e^{E(p)/T}`$. This is the generic behavior, which can be avoided by choosing very unnatural initial conditions only. At this point the goal we have set for ourselves at the beginning of this paper has been achieved. A rather robust quantum model has been constructed illustrating the following conjecture: the evolution of the wave function of a system of very many interacting particles is likely to drive the one-particle spectrum towards a ”thermal” shape. This behavior, which we would like to call pseudo-thermalization, is likely but not mandatory, in the model it occurs for a (wide) range of model parameters only.

Notice, that the physical contents of the classical and of the quantum model are different. In the classical model a particle keeps changing its momentum. In the quantum model the components of a wave keep splitting. All these components are a priori present at any time, only their amplitudes get modified. However, in both models the observable one-particle spectrum has the same ”thermal” shape. Moreover, both models have the Markov property: the generic limiting behavior is independent of the initial state of the system.

As long as one stays within the one-particle framework, several specific assumptions of our quantum model can be relaxed. Still, even in this framework, it would be interesting to check how generic is the result obtained in this paper, viz. does it hold for processes which are significantly different from (1). The difficulty is mainly technical, at least if one looks for an analytic solution. However, going beyond the one-particle picture and checking the conjecture in a genuine multiparticle context requires a significant conceptual advance. On the theory side, we made just a first step and much remains to be done.

As mentioned at the beginning of this text, our research has actually a phenomenological motivation. Indeed, in high energy processes involving hadrons the inclusive spectra often look thermal, in circumstances where genuine thermalization seems very unlikely. For example, in soft hadronic collisions the inclusive transverse momentum spectrum at not too large $`p_{}`$ has an approximately boltzmannian shape,

$$\frac{d\sigma }{d^2p_{}}e^{E(p_{})/T},$$
(18)
as if it resulted from a thermalization process occuring during the collision. Here, $`E(p_{})`$ denotes the ”transverse energy” of a secondary foot1  and $`T`$ is a ”temperature” parameter, roughly independent of the secondary hadron species. However, the thermalization in the final state, debatable but plausible in a heavy-ion collision, is nearly impossible in a minimum bias nucleon-nucleon collision, where (18) is the best observed.

As is well known, in the infinite momentum frame $`P\mathrm{}`$ a hadron can be regarded as a 2D nonrelativistic, coherent multiparton state enclosed in a box, with dimensions fixed by the confinement scale. The kinetic energy of a parton is $`p_{}^2/2xP`$, $`x`$ denoting its longitudinal momentum fraction. It is plausible that pseudo-thermalization holds for the multiparton wave function. If so, then one can argue, invoking the parton-hadron duality, that after the hadron is broken by the interaction the inclusive spectrum reflects the ”thermal” distribution holding at the parton level foot2 .

Notice, that since $`xP`$ plays here the role of ”mass”, one expects that the average transverse momentum is larger for larger longitudinal momentum (in cms). Such a correlation is indeed observed. On the other hand it is not expected when the Schwinger mechanism is operating.

Another example: in ref. bhal  a surprisingly good and economic description of structure functions is presented, using a statistical ansatz. This might be another example of pseudo-thermalization, this time oberved directly at the parton level.

On the phenomenology side, our exercise supports the idea that the ”thermal” spectra observed in high-energy processes involving hadrons do actually reflect pseudo-thermalization affecting the structure of the multiparton wave functions of the incoming hadrons. If true, this is to large extent independent of the details of the microscopic dynamics - a common feature of many multiparticle production phenomena - except that the validity of the parton picture requires an asymptotically free theory, of course foot5 .

Acknowledgements: I am indebted to Roger Balian, Andrzej Bialas, Ian Kogan and Jean-Pierre Leroy for helpful comments and suggestions. This work was partially supported by EC IHP grant HPRN-CT-1999-000161. Laboratoire de Physique Théorique is Unité Mixte du CNRS UMR 8627.

Note added: Preparing the original version of this note I have completedy overlooked results obtained in the context of so-called quantum chaos theory, and in particular the early paper by M. Srednicki sred , much more profound that the heuristic discussion presented above. Working with a gas of hard spheres, he has given strong arguments in favor of the conjecture that one-particle distribution in a stationary pure quantum state is Boltzmann, provided the corresponding classical state is chaotic.


text:
# Orbital Magnetism in Three-Dimensional Quantum Dots

## Abstract

We study orbital magnetism in a three-dimensional (3D) quantum dot with a parabolic confining potential. We calculate the free energy of the system as a function of the magnetic field and the temperature. By this, we show that the temperature-field plane can be classified into three regions in terms of the characteristic behavior of the magnetization: the Landau diamagnetism, de Haas-van Alphen oscillation and mesoscopic fluctuation of magnetization. We also calculate numerically the magnetization of the system and then the current density distribution. As for the oscillation of the magnetization when the field is varied, the 3D quantum dot shows a longer period than a 2D quantum dot which contains the same number of electrons. A large paramagnetism appears at low temperatures when the magnetic field is very weak.

Orbital magnetism of a bulk system appears as the Landau diamagnetism (LD). In the past, Landau’s work has provoked discussion regarding his treatment of boundary rf:1 ; rf:2  and the effects of boundary on the orbital magnetism have been studied by many people. rf:3 ; rf:4  Recently, two-dimensional (2D) quantum dots have become experimentally available with the capability of controlling the dots’ size, shape and number of electrons. rf:5 ; rf:6 ; rf:7  Stimulated by experimental development, many theorists have been studying on 2D quantum dots with parabolic confining potentials. In particular, Ishikawa and Fukuyama (IF) rf:3  have claimed that the temperature-field ($`T`$-$`B`$) plane can be classified into three regions in terms of the characteristic behavior of the magnetization $`M(T,B)`$: the LD, de Haas-van Alphen (dHvA) oscillation and mesoscopic fluctuation (MF) of magnetization. Such magnetic behavior is mainly governed by the energy levels near the Fermi level of the system. In the case of 3D quantum dots, the additional degree of freedom makes the density of states larger than that of 2D quantum dots. Hence, it is expected that the temperature affects the magnetization more easily, which changes the magnetic behavior to some extent and possibly invalidates the classification.

In this letter, we present our study on the orbital magnetism of 3D quantum dots. We consider noninteracting electrons confined in a parabolic potential and calculate analytically the free energy $`\mathrm{\Omega }`$ of the system. Thereby, we show that the temperature-field plane can still be classified into the above three regions in terms of the characteristic behavior of the magnetization, $`M(T,B)=\mathrm{\Omega }/B`$. In order to confirm this classification, we calculate numerically the magnetization of the system and then the current density distribution. In addition to the consistency between the analytical and numerical calculations, we report some remarkable findings of our numerical calculation. Hereafter, we set the fundamental constants as $`\mathrm{}=c=k_\mathrm{B}=1`$.

We consider spinless electrons of effective mass $`m`$ and charge $`e`$ confined in a spherical quantum dot with a parabolic confining potential in a uniform magnetic field $`𝐁`$. The total Hamiltonian $`H`$ of the system is given by the sum of single electron Hamiltonians as

$`H`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{\frac{1}{2m}\left[\mathrm{i}_i+e𝐀(𝐫_i)\right]^2+\frac{m}{2}\omega _3^2𝐫_i^2\right\},`$ (1)
where $`N`$ is the total number of electrons, and $`𝐀`$ is a vector potential satisfying $`𝐁=\mathrm{rot}𝐀`$. The last term represents the confining potential, where $`\omega _3`$ parametrizes the strength of the isotropic confinement.

The eigenstates of the single electron Hamiltonian comprise the Fock-Darwin states in the $`x`$-$`y`$ plane rf:3  and the eigenstates of a harmonic oscillator in the $`z`$-direction parallel to $`𝐁`$. In the cylindrical coordinate $`(\rho ,\varphi ,z)`$, the eigenfunction $`\psi `$ is given by

$$\psi _{n\alpha \nu }(𝐫)=\left(2\pi \right)^{\frac{1}{2}}\mathrm{exp}(\mathrm{i}\alpha \varphi )R_{n\alpha }(\rho )Z_\nu (z),$$
(2)
where $`n`$, $`\alpha `$, $`\nu `$ are integers and $`n`$, $`\nu 0`$. $`(n,\alpha ,\nu )`$ consists of a set of quantum numbers. $`R_{n\alpha }(\rho )`$ is the radial wave function of the Fock-Darwin state specified by $`(n,\alpha )`$ and $`Z_\nu (z)`$ is the $`\nu `$-th eigenfunction of the harmonic oscillator. The eigenenergy $`E_{n\alpha \nu }`$ is given by

$`E_{n\alpha \nu }=\left(n+\frac{1}{2}\right)\omega _1+|\alpha |\omega _2+\left(\nu +\frac{1}{2}\right)\omega _3,`$ (3)
$`\omega _1=\sqrt{\omega _B^2+4\omega _3^2},`$ (4)
$`\omega _2(\pm )=\omega _\pm =\frac{1}{2}(\omega _1\pm \omega _B),`$ (5)
where $`\omega _B`$ denotes the cyclotron frequency, $`\omega _B=eB/m`$. In eq. (5), the characteristic frequency $`\omega _2=\omega _2(\sigma )`$ depends on the sign $`\sigma (=\pm )`$ of the angular momentum $`\alpha `$. Hereafter, we set $`\omega _3`$ as the energy unit, though we keep the symbol $`\omega _3`$ in the several equations below for clarity.

The magnetization $`M`$ of the entire system is given by $`M=(\mathrm{\Omega }/B)_{T\mu }`$, rf:8  where $`T`$ is the temperature, $`\mu `$ is the chemical potential, and $`\mathrm{\Omega }`$ is the free energy of the system, $`\mathrm{\Omega }=T_{n\alpha \nu }\mathrm{ln}\{1+\mathrm{exp}[(\mu E_{n\alpha \nu })/T]\}`$. By means of the Poisson summation formula, rf:9  the triple sum with respect to $`n,\alpha ,\nu `$ is transformed into another triple sum, $`\mathrm{\Omega }=_{p,q,r=\mathrm{}}^{\mathrm{}}\mathrm{\Omega }_{pqr},`$ where $`\mathrm{\Omega }_{pqr}`$ is expressed as the Fourier integral which we can calculate analytically by assuming $`\mu T`$ and $`\mu \omega _B`$. With these calculations carried out, the free energy $`\mathrm{\Omega }`$ is given as follows,

$`\mathrm{\Omega }\mathrm{\Omega }_0+\stackrel{~}{\mathrm{\Omega }},`$ (6)
$`\mathrm{\Omega }_0=\frac{1}{24}\mu ^4+\left(\frac{1}{48}\omega _B^2\frac{\pi ^2}{12}T^2+\frac{1}{16}\right)\mu ^2,`$ (7)
$`\stackrel{~}{\mathrm{\Omega }}={\displaystyle \underset{i=1}{\overset{3}{}}}\mathrm{\Omega }_i={\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \underset{\sigma =\pm }{}}{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\mathrm{\Omega }_i^\sigma (p),`$ (8)
$`\mathrm{\Omega }_i^\sigma (p)=Q_i^\sigma (p)\mathrm{\Psi }\left(2p\pi ^2T/\omega _i\right)\mathrm{cos}\left(2p\pi \mu /\omega _i\right),`$ (9)
$`Q_i^\sigma (p)={\displaystyle \frac{(1)^{ip}\omega _i^3}{8\pi ^4p^4\omega _j\omega _k}}`$
$`{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^{ip}}{4\pi ^4p^2}}\left[{\displaystyle \frac{(1)^{jq}\omega _i^3\omega _j}{\omega _k(p^2\omega _j^2q^2\omega _i^2)}}+{\displaystyle \frac{(1)^{kq}\omega _i^3\omega _k}{\omega _j(p^2\omega _k^2q^2\omega _i^2)}}\right]`$
$`{\displaystyle \underset{q=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{r=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^{ip+jq+kr}\omega _i^3\omega _j\omega _k}{2\pi ^4(p^2\omega _j^2q^2\omega _i^2)(p^2\omega _k^2r^2\omega _i^2)}},`$ (10)
where $`\mathrm{\Psi }(x)=x/\mathrm{sinh}(x)`$, and in eq.(10), the coefficients $`Q_1^\sigma ,Q_2^\sigma ,Q_3^\sigma `$ take the set of indices $`(i,j,k)=(1,2,3),(2,3,1),(3,1,2)`$, respectively. The summation in eq.(10) is executed further and we obtain

$`Q_1^\sigma (p)`$ $`=`$ $`{\displaystyle \frac{(1)^p\omega _1}{8\pi ^2p^2\mathrm{tan}\left(\pi p\omega _2/\omega _1\right)\mathrm{sin}\left(\pi p\omega _3/\omega _1\right)}},`$ (11)
$`Q_2^\sigma (p)`$ $`=`$ $`{\displaystyle \frac{\omega _2}{8\pi ^2p^2\mathrm{sin}\left(\pi p\omega _3/\omega _2\right)\mathrm{sin}\left(\pi p\omega _1/\omega _2\right)}},`$ (12)
$`Q_3^\sigma (p)`$ $`=`$ $`{\displaystyle \frac{(1)^p\omega _3}{8\pi ^2p^2\mathrm{sin}\left(\pi p\omega _1/\omega _3\right)\mathrm{tan}\left(\pi p\omega _2/\omega _3\right)}}.`$ (13)
The free energy $`\mathrm{\Omega }_0`$ in eq. (7) is the nonoscillating part with respect to the magnetic field, while $`\stackrel{~}{\mathrm{\Omega }}=\mathrm{\Omega }_1+\mathrm{\Omega }_2+\mathrm{\Omega }_3`$ in eq. (8) represents the oscillating part of the total free energy $`\mathrm{\Omega }`$. For the oscillating part, the free energy $`\mathrm{\Omega }_1=_\sigma _p\mathrm{\Omega }_1^\sigma (p)`$ vanishes because $`_\sigma \mathrm{\Omega }_1^\sigma (p)_\sigma Q_1^\sigma (p)`$ is identically zero due to eq. (5). As for $`\mathrm{\Omega }_2`$, $`\mathrm{\Omega }_2^\pm (p)`$ has the oscillating factor, $`\mathrm{cos}(2p\pi \mu /\omega _\pm )`$, in eq. (9) as a function of the field parameter $`\omega _B`$ since $`\omega _\pm `$ varies monotonically with respect to $`\omega _B`$. Although the remaining factor $`Q_2^\pm (p)`$ in $`\mathrm{\Omega }_2^\pm (p)`$ also oscillates by $`\mathrm{sin}(\pi p\omega _3/\omega _\pm )\mathrm{sin}(\pi p\omega _1/\omega _\pm )`$ in eq. (12), we can neglect this oscillation, compared with $`\mathrm{cos}(2p\pi \mu /\omega _\pm )`$, since we can assume $`\omega _3\mu `$ for the sufficiently large size of the system and also $`\omega _1\mu `$ due to the assumption, $`\mu \omega _B`$. As for $`\mathrm{\Omega }_3`$, $`\mathrm{\Omega }_3^\sigma (p)`$ in eq. (9) has only a negligible oscillation caused by $`\mathrm{sin}(\pi p\omega _1/\omega _3)\mathrm{tan}(\pi p\omega _2/\omega _3)`$ in eq. (13), because $`\mathrm{cos}(2p\pi \mu /\omega _3)`$ is independent of the magnetic field. Consequently, $`\mathrm{\Omega }_2`$ is dominant in eq. (8) in terms of the oscillation of magnetization. Hereafter, we write the effective part of the free energy $`\stackrel{~}{\mathrm{\Omega }}`$ as

$`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}={\displaystyle \underset{\sigma =\pm }{}}{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}Q_2^\sigma (p)\mathrm{\Psi }\left(2p\pi ^2T/\omega _\sigma \right)\mathrm{cos}\left(2p\pi \mu /\omega _\sigma \right),`$ (14)
where $`\omega _\sigma =\omega _2(\sigma )`$.

It should be noted, however, that eqs. (11)-(13) have spurious divergent singularities when $`\omega _i/\omega _j`$ is a rational number. To eliminate such divergences, we rewrite $`\stackrel{~}{\mathrm{\Omega }}`$ as defined in eqs. (8)-(10), by rearranging terms, into

$`\stackrel{~}{\mathrm{\Omega }}={\displaystyle \underset{\sigma =\pm }{}}\left({\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\widehat{\mathrm{\Omega }}_p^\sigma +{\displaystyle \underset{p,q=1}{\overset{\mathrm{}}{}}}\widehat{\mathrm{\Omega }}_{p,q}^\sigma +{\displaystyle \underset{p,q,r=1}{\overset{\mathrm{}}{}}}\widehat{\mathrm{\Omega }}_{p,q,r}^\sigma \right),`$ (15)
$`\widehat{\mathrm{\Omega }}_p^\sigma ={\displaystyle \frac{(8\pi ^4)^1}{\omega _1\omega _2\omega _3}}{\displaystyle \underset{i=1}{\overset{3}{}}}(1)^{ip}{\displaystyle \frac{Y\left(p/\omega _i\right)}{(p/\omega _i)^4}},`$ (16)
$`\widehat{\mathrm{\Omega }}_{p,q}^\sigma ={\displaystyle \frac{(4\pi ^4)^1}{\omega _1\omega _2\omega _3}}{\displaystyle \underset{(i,j)}{}}(1)^{ip+jq}{\displaystyle \frac{\frac{Y\left(p/\omega _i\right)}{(p/\omega _i)^2}\frac{Y\left(q/\omega _j\right)}{(q/\omega _j)^2}}{(p/\omega _i)^2(q/\omega _j)^2}},`$ (17)
$`\widehat{\mathrm{\Omega }}_{p,q,r}^\sigma ={\displaystyle \frac{(2\pi ^4)^1}{\omega _1\omega _2\omega _3}}(1)^{p+r}`$
$`\times {\displaystyle \frac{\frac{Y\left(p/\omega _1\right)Y\left(r/\omega _3\right)}{(p/\omega _1)^2(r/\omega _3)^2}\frac{Y\left(q/\omega _2\right)Y\left(r/\omega _3\right)}{(q/\omega _2)^2(r/\omega _3)^2}}{(p/\omega _1)^2(q/\omega _2)^2}},`$ (18)
where $`Y(x)=\mathrm{\Psi }(2\pi ^2Tx)\mathrm{cos}(2\pi \mu x)`$, and in eq.(17), the index $`(i,j)`$ takes $`(1,2),(2,3)`$ and $`(3,1)`$. Equations (16)-(18) are apparently free from divergences. It is expected that as the temperature $`T`$ rises, the total amplitude of oscillation is rapidly suppressed by the exponentially decreasing factor, $`\mathrm{\Psi }(2p\pi ^2T/\omega _i)(4p\pi ^2T/\omega _i)\mathrm{exp}(2p\pi T/\omega _i)`$.

At sufficiently high temperatures, it is expected that $`\mathrm{\Omega }_0`$ in eq. (7) makes a dominant contribution to the total free energy $`\mathrm{\Omega }`$, and the magnetization $`M_0`$ is given by

$$M_0/\mu _\mathrm{B}=2\left(\mathrm{\Omega }_0/\omega _B\right)_{T\mu }=\frac{1}{12}\mu ^2\omega _B,$$
(19)
where $`\mu _\mathrm{B}=e/(2m)`$ is the Bohr magneton. The magnetization $`M_0`$ is negative and its absolute value increases linearly with increasing $`\omega _B`$, which is characteristic of the LD in a bulk system. rf:9  On the other hand, at low temperatures, by using the free energy $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}`$ in eq. (14), we obtain an oscillating part of the magnetization, $`\stackrel{~}{M}/\mu _\mathrm{B}=2(\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}/\omega _B)_{T\mu }`$, which contains the trigonometric functions of $`2\pi p\mu /\omega _+`$ and $`2\pi p\mu /\omega _{}`$. It should be noted that as the temperature $`T`$ decreases, the amplitude of the oscillation of $`\stackrel{~}{M}`$ becomes larger.

In the following, we discuss how the three regions introduced by IF rf:3  are modified in 3D. We regard the LD region as the region where the amplitudes of the oscillations specified by $`p`$ are much smaller than the confining energy $`\mathrm{}\omega _3`$, which is unity in our notation. Hence, the criterion for the LD region is given by $`|Q_2^\pm (p)|\mathrm{\Psi }(2p\pi ^2T/\omega _\pm )1`$ for all $`p`$. Introducing the cutoff $`ϵ(1)`$, we defined the LD region as the region where the following two inequalities hold simultaneously,

$$\mathrm{\Psi }\left(2p\pi ^2T/\omega _+\right)ϵ,\mathrm{\Psi }\left(2p\pi ^2T/\omega _{}\right)ϵ.$$
(20)
In this letter, we set $`ϵ=0.1\%`$, but it will be necessary, however, to reexamine the value of $`ϵ`$ when our analytic classification is compared to any experiment. Noting $`\mathrm{\Psi }^{}(x)<0`$ and $`\omega _+>\omega _{}`$, the region determined by eq. (20) is equivalent to $`T\gamma \omega _+`$, where $`\gamma =\mathrm{\Psi }^1(ϵ)/(2\pi ^2)0.5`$. As a result, the boundary between the LD and the oscillatory region is given by $`T=0.5\omega _+`$, as plotted by the solid line in Fig.1.

The oscillatory region $`T<\gamma \omega _+`$ is divided into two subregions as shown in Fig.1; (A) $`\gamma \omega _+>T\gamma \omega _{}`$, and (B) $`\gamma \omega _{}>T`$. In subregion (A), the first inequality in eq. (20) fails while the second one still holds. Hence, only $`Q_2^+(p)\mathrm{\Psi }(2\pi ^2T/\omega _+)\mathrm{cos}(2\pi \mu /\omega _+)`$ contributes to $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}`$. The effective free energy $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}^\mathrm{A}`$ is given by $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}^\mathrm{A}=_{p=1}^{\mathrm{}}Q_2^+(p)\mathrm{\Psi }\left(2p\pi ^2T/\omega _+\right)\mathrm{cos}\left(2p\pi \mu /\omega _+\right)`$. On the other hand, both inequalities fail in subregion (B), where two types of oscillations regarding $`\omega _\pm `$ contribute to $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}`$, and the free energy $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}^\mathrm{B}`$ is given by eq. (14). As the temperature $`T`$ is lowered, the terms with larger $`p`$ start contributing to $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}`$.

For subregion (A), magnetization $`\stackrel{~}{M}`$ derived from $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}^\mathrm{A}`$ contains a trigonometric function having the period $`\omega _+/\mu `$. The dHvA region corresponds to the region $`\gamma \omega _{}T<\gamma \omega _+`$, characterized by the oscillation having the period $`\omega _+/\mu `$. Particularly in a strong field region, since $`\omega _+\omega _B`$, magnetization $`\stackrel{~}{M}`$ behaves as the periodic function of variable $`\omega _B^1`$. This periodicity is typical of the dHvA oscillation in a bulk system. rf:8 ; rf:9

For subregion (B), magnetization $`\stackrel{~}{M}`$ derived from $`\stackrel{~}{\mathrm{\Omega }}_{\mathrm{eff}}^\mathrm{B}`$ contains a trigonometric function having the period $`\omega _{}/\mu `$ in addition to $`\omega _+/\mu `$. In this region, $`\omega _+`$ and $`\omega _{}`$ are not so close to each other and the two modes fluctuate separately. Consequently, the coexistence of the above two modes causes the rapid and irregular oscillation of magnetization. This magnetic behavior corresponds to the MF. Hence, subregion (B) is regarded as the MF region.

However, when the magnetic field is extremely weak in the MF region, the periods of $`\omega _\pm /\mu `$ become much closer since $`\omega _\pm 1\pm \frac{1}{2}\omega _B`$. It is impossible to regard $`\omega _\pm `$ as distinguishable and therefore the separation of the two fluctuations breaks down in the above argument. In this case, the largest contribution to $`\stackrel{~}{\mathrm{\Omega }}`$ in eq. (15) comes from the modes satisfying $`p:q:r=\omega _1:\omega _2:\omega _3`$, namely, $`p/\omega _1q/\omega _2=q/\omega _2r/\omega _3=r/\omega _3p/\omega _1=0`$ for $`\widehat{\mathrm{\Omega }}_{p,q,r}^\sigma `$ in eq. (18). From the behavior of characteristic frequency $`\omega _12+\frac{1}{4}\omega _B^2`$, we can see that the modes of $`p=2j`$, $`q=j`$, $`r=j`$ ($`j`$ is a natural number) make a dominant contribution to magnetization $`\stackrel{~}{M}`$. The largest contribution to $`\stackrel{~}{M}`$ is evaluated as

$`{\displaystyle \frac{\stackrel{~}{M}}{\mu _\mathrm{B}}}{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\widehat{M}_{2j,j,j}}{\mu _\mathrm{B}}}=2{\displaystyle \underset{\sigma =\pm }{}}{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\widehat{\mathrm{\Omega }}_{2j,j,j}^\sigma }{\omega _B}}\right)_{T\mu }.`$ (21)
Furthermore, we can estimate the upper bound $`U(p,q,r)`$ of $`|\widehat{\mathrm{\Omega }}_{p,q,r}^\sigma |`$ by using the triangle inequality and $`|Y(x)|1`$ for eq.(18), and it is found that the upper bound $`U(2j,j,j)`$ decreases as $`U(2j,j,j)=j^4U(2,1,1)`$ for larger $`j`$. Hence, it is expected that a few modes with small $`j`$ make a major contribution. The mode with $`j=1`$, $`\widehat{M}_{2,1,1}`$, is plotted in Fig.2(a), where we take $`T=0.01`$ and $`N=5000`$. The magnetization curve in Fig.2(a) shows a remarkable paramagnetism at weak magnetic fields below $`\omega _B0.02`$. The region corresponding to the paramagnetism is indicated by P in Fig.1.

Next, we describe the procedure of the numerical calculation for magnetization. We consider magnetization as the total dipole moment of electrons with the confining potential. From this standpoint, magnetization $`𝐌`$ is given by $`𝐌=_{\mathrm{}}dV\frac{1}{2}\left[𝐫\times 𝐉(𝐫)\right]`$, where $`𝐉(𝐫)`$ is the current density distribution of the system. The calculation of $`𝐌`$ reduces to that of $`𝐉(𝐫)`$ which is given by $`𝐉(𝐫)=\mathrm{Re}_{s=(n,\alpha ,\nu )}f(E_s)\psi _s^{}(𝐫)\frac{e}{m}\left(\mathrm{i}+e𝐀\right)\psi _s(𝐫)`$. Here, $`f(E)`$ denotes the Fermi distribution function, $`f(E)=\{\mathrm{exp}[(E\mu )/T]+1\}^1`$. By using eq.(2), the current density distribution leads to $`𝐉(𝐫)=J(\rho ,z)𝐞_\varphi `$, and

$`J(\rho ,z)=\sqrt{{\displaystyle \frac{e}{m}}}{\displaystyle \underset{s=(n,\alpha ,\nu )}{}}f(E_s)\left|\psi _s(𝐫)\right|^2\left(\alpha {\displaystyle \frac{\xi }{\rho }}+{\displaystyle \frac{\omega _B}{2}}{\displaystyle \frac{\rho }{\xi }}\right),`$ (22)
where $`\xi =m^{\frac{1}{2}}`$. By using eq.(22), the total dipole moment $`𝐌`$ is calculated analytically rf:3  and then $`𝐌=M𝐞_z`$. The $`z`$-component $`M`$ is given by

$`M/\mu _\mathrm{B}={\displaystyle \underset{s=(n,\alpha ,\nu )}{}}f(E_s)\left[\alpha +\left(2n+|\alpha |+1\right)\omega _B/\omega _1\right],`$ (23)
which is also derived from $`M/\mu _\mathrm{B}=2(\mathrm{\Omega }/\omega _B)_{T\mu }=2_sf(E_s)(E_s/\omega _B)`$ due to eq.(3). This coincidence of two derivations, of course, occurs regardless of the type of confining potentials. In the numerical calculation, we calculate the chemical potential first by solving $`N=_sf(E_s)`$ for $`\mu `$, and then compute magnetization $`M`$ by using eq.(23).

Magnetization $`M`$ calculated numerically is shown in Fig.2 for $`N=5000`$. At low temperatures ($`T=0.01`$, Figs.2(b) and 2(c)), the magnetization oscillates over the entire range of magnetic fields. At extremely weak fields below $`\omega _B0.1`$ (Fig.2(b)), the magnetization oscillates smoothly with significantly large amplitude. The smoothness implies that not so many modes are dominant in the oscillation. Moreover, Fig.2(b) agrees qualitatively with the behavior of $`\widehat{M}_{2,1,1}(T,\mu ,B)`$ (Fig.2(a)) with the chemical potential $`\mu (T,N,B)`$ calculated numerically for $`T=0.01`$ and $`N=5000`$. The large smooth oscillation in Fig.2(b) is consistent with our analysis in terms of free energy. At weak magnetic fields above $`\omega _B0.1`$, the large oscillation is suppressed and the magnetization irregularly takes a positive or negative value. This fluctuation continues until $`\omega _B5`$ with the linear magnetization superposed (Fig.2(c)). At high magnetic fields above $`\omega _B5`$, the large smooth oscillation is clear in Fig.2(c). It is considered that this magnetization curve reflects the dHvA oscillation. On the other hand, at high temperatures ($`T=2`$, Fig.2(d)), the magnetization is linear in $`\omega _B`$ until $`\omega _B10`$ and then the dHvA oscillation increases gradually.

Compared with the magnetization in 2D reported by IF, rf:3  it is found that the period of oscillation is longer in 3D. For the 2D system, they also calculated the free energy analytically and its oscillating part is expressed as $`\stackrel{~}{\mathrm{\Omega }}=_{i=1}^2_\sigma _p\mathrm{\Omega }_i^\sigma (p)`$, and $`\mathrm{\Omega }_i^\sigma (p)=Q_i^\sigma (p)\mathrm{\Psi }(2p\pi ^2T/\omega _i)\mathrm{sin}(2p\pi \mu /\omega _i)`$, in our notation. Here, we find that for 2D and 3D, these oscillating factors have the common phase $`2p\pi \mu /\omega _i`$. Noting the common phase, when $`\omega _B`$ is varied, the period with respect to $`\omega _i^1`$, $`\mathrm{\Delta }(1/\omega _i)`$, is given by $`\mathrm{\Delta }(1/\omega _i)=p^1\mu ^1`$. Hence, $`\mu ^1`$ serves as the scale for the period of oscillation and the difference between the periods can be explained by that between the chemical potentials. The chemical potentials $`\mu _{3\mathrm{D}},\mu _{2\mathrm{D}}`$ are estimated by $`N(\mathrm{\Omega }_0/\mu )_{TB}`$ and then $`\mu _{3\mathrm{D}}\sqrt[3]{6N}<\mu _{2\mathrm{D}}\sqrt{2N}`$ for the same number of electrons. Consequently, the period of oscillation in 3D is longer than that in 2D.

The current density distribution $`J(\rho ,z)`$ is computed by using eq.(22) and is shown in Fig.3 for $`N=5000`$. The positive value of the vertical axis corresponds to diamagnetism. When the remarkable paramagnetism appears (Fig.3(a)), a paramagnetic region for $`J(\rho ,z)`$ extends throughout the spherical quantum dot except for a small region near the $`z`$-axis ($`\rho =0`$). In the MF region (Fig.3(b)), the paramagnetic and diamagnetic regions are distributed irregularly. In the LD region (Fig.3(c)), the current density distribution is localized near the surface of the quantum dot. This surface current corresponds to the edge current obtained by IF rf:3  for the 2D system. In the dHvA region (Fig.3(d)), the current density distribution oscillates smoothly with a large amplitude.

In conclusion, we have studied orbital magnetism in a 3D quantum dot with a parabolic confining potential. From the analysis of the free energy, we have determined the boundaries between the LD, dHvA, and MF regions. We have also computed the magnetization of the system and then the current density distribution, which indeed reflects the character of each region. The results of numerical calculation show a large paramagnetism when the magnetic fields are extremely weak in the MF region. The dHvA oscillation has a longer period in 3D than that in a 2D quantum dot. These results are explained in terms of the free energy. Furthermore, we have found surface current in the LD region, corresponding to edge current for 2D quantum dots.

The authors are grateful to H. Fukuyama and T. Sasaki for their useful discussions. H. I. is supported by MEXT, Grant-in-Aid for Encouragement of Young Scientists, 13740197.


text:
# Untitled Document

A Percolative Model of Soft Breakdown in Ultrathin Oxides

1C. Pennetta, 1L. Reggiani and 2Gy. Trefán

1INFM - National Nanotechnology Laboratory, Dipartimento di Ingegneria dell’Innovazione

Università di Lecce, Via Arnesano, I-73100, Italy

2Eindhoven University of Technology, Dept. of Electrical Engineering

5600 MB Eindhoven, The Netherlands

Abstract

The degradation of ultrathin oxide layers in the presence of a stress voltage is modeled in terms of two antagonist percolation processes taking place in a random resistor network. The resistance and leakage current fluctuations are studied by Monte Carlo simulations for voltages below the breakdown threshold. An increase of excess noise together with a noticeable non-Gaussian behavior is found in the pre-breakdown regime in agreement with experimental results.

Keywords: breakdown, dielectrics, Monte Carlo, percolation

1. Introduction and model

Soft breakdown occurs in ultrathin oxide layers during constant voltage stress. This process manifests itself by a large increase of the stress-induced leakage currents (SILC) and by the occurrence of giant current fluctuations associated with resistance fluctuations. These fluctuations can be related to trapping and detrapping of electrons in the percolation cluster formed between the electrodes of the capacitor at the soft breakdown . In this paper, we investigate the resistance in the pre-breakdown region and the associated fluctuations using two percolations which evolve in competition on a random resistor network , as briefly summarized below. We describe a thin film as a two-dimensional square-lattice network of resistors of resistance $`r_n`$, laying on an insulating substrate at temperature $`T_0`$ acting as a thermal reservoir. Geometrical effects arising from different shapes are neglected. Initially all the resistors are identical, $`r_n=r_0`$. We take a square geometry $`N\times N`$ where $`N`$ determines the linear sizes of the lattice and $`N_{tot}=2N^2`$ is the total number of resistors. Electrical contacts, realized by perfectly conducting bars, are placed at the left and right hand sides of the network where a constant voltage $`V`$ is applied. By increasing the applied voltage, defects are introduced by replacing single resistors with short circuits. Thus, the network degrades by undergoing an insulator-conductor transition when the defects gradually become so dense that finally they form a continuous short circuit path between the contacts . The temperature $`T_n`$ of the $`n`$-th resistor is computed by :

$$T_n=T_0+A\left[r_ni_n^2+\frac{B}{N_{neig}}\underset{l=1}{\overset{N_{neig}}{}}\left(r_li_l^2r_ni_n^2\right)\right]$$
$`(1)`$
where the term which adds to $`T_0`$ in the r.h.s. of Eq. (1) accounts for Joule heating associated with the current flowing in the resistor and the coupling with the nearest neighbours. Here, $`N_{neig}`$ is the number of first neighbours around the $`n`$-th resistor with current $`i_n`$, $`A`$ the heat coupling parameter and $`B=3/4`$ provides a uniform heating everywhere in the perfect network. We note that Eq. (1) assumes an instantaneous thermalisation between neighbor resistors. Probabilities of creating and recovering a defect are taken as $`W_D=exp(E_D/K_BT_n)`$ and $`W_R=exp(E_R/K_BT_n)`$ with $`E_D`$ and $`E_R`$ characteristic activation energies, respectively, and $`K_B`$ the Boltzmann constant. Monte Carlo simulations are carried out using the following procedure. (i) Starting from the perfect lattice, we calculate the total network resistance $`R`$, the external current $`I`$, and the local currents by solving Kirchhoff’s loop equations by the Gauss elimination method. The local temperatures are then calculated and used for the successive update of the network. (ii) The defects are generated with probability $`W_D`$. The local currents and the local temperatures are then recalculated. (iii) The defects are recovered with probability $`W_R`$, and the total network resistance, the local currents, and the temperatures are finally calculated. This procedure is iterated from (ii). Each iteration step can be associated with an elementary time step. The iteration runs until either the defect percolation threshold is reached (i.e. the network resistance drops below $`10^3`$ times to its initial value), or steady-state is achieved. In this latter case, the iteration is let to run long enough for a fluctuation analysis to be carried out.

2. Results and discussion

We compare the results of numerical simulations to the corresponding experimental observations and validate the present approach by showing that it reproduces the key experimental features. Since soft breakdown experiments report leakage current evolutions, noise spectra and probability density of fluctuations, we investigate the same quantities by numerical simulations. In the simulations we applied a constant voltage stress with the following realistic parameters $`N=75`$, $`r_0=10^7\mathrm{\Omega }`$, $`T_0=300K`$, $`A=5\times 10^7{}_{}{}^{o}C/W`$, $`E_D=0.19eV`$, $`E_R=0.13eV`$. Figure 1 reports the leakage current evolutions for stress voltages, respectively, of $`0.1`$, $`0.3`$, $`0.5`$, $`0.6`$, $`0.7`$, $`0.8`$, $`0.85`$ and $`0.9V`$. We note that the values of the activation energies $`E_D`$ and $`E_R`$ control the speed of the evolution toward steady state and/or breakdown and as such should be calibrated from a quantitative comparison with experiments. However, the qualitative features of the breakdown process we are presently interested in do not depend significantly from their absolute values. We found that for voltages lower than $`0.9V`$ no breakdown occurs but a dynamical balance takes place between defect generation and defect recovery. The steady state so achieved exhibits current fluctuations whose statistical and spectral features characterize the electrical quality of the dielectric. However, for stress voltages equal or larger than $`0.9V`$ after a few thousand iteration steps the leakage current grows more than a thousand fold which indicates a breakdown situation. We also observe that the breakdown does not occur as a result of a systematic current growth but more as a result of current bursts. In fact, a set of bursts generally occurs that lead to the breakdown of the insulator. Thus an interesting pre-breakdown region appears like in experiments . From a detailed analysis of a set of simulations we also conclude that a threshold voltage exists below which no breakdown occurs. Figure 2 shows how the breakdown voltage depends on the resistance of the virgin dielectric, modeled as an initial perfect lattice with increasing value of the elemental resistor. This breakdown voltage is found to increase proportional to the square root of the perfect network resistance. Figure 3 shows the steady state current noise spectrum calculated by FFT for stress voltages ranging from $`0.1`$ to $`0.85V`$. Within numerical uncertainity, we detect Lorentzian spectra of the same corner frequency which indicates that the characteristic times of fluctuations are independent of stress voltages. Interestingly, the strong increase of the value of the plateaux at increasing stress voltages shows a super-quadratic increase of the noise with the applied stress. The current-noise spectra in Fig. 3 are found to be in qualitative agreement with experiments . Figure 4 reports the normalized variance of resistance fluctuations as a function of the applied voltage. Starting from the intrinsic value of the network , the variance is found to increase significantly as a net effect of the stressing voltage which is ultimately responsible of the breakdown at about $`0.9V`$. Recent observations of SILC measurements reported non-Gaussian current fluctuations in the soft breakdown region of ultra-thin dielectrics . Therefore, we computed the distribution function of the current fluctuations $`P(I)`$ for stress voltages of $`0.5V`$ and $`0.85V`$ and plotted them in Fig. 5 as a function of $`I<I>`$, where $`<I>`$ is the average current value. On this linear-log plot the corresponding Gaussian distribution would be a downward parabola. The distribution for an applied voltage near breakdown of $`0.85V`$ clearly exhibits a non-Gaussian tail at high currents over the average due to the large upward bursts which are not balanced by downward bursts in the current evolution. Since a further increase of the stress voltage causes the film to break down, the emerging non-Gaussian distribution can be thought of as a precursor of failure as observed in experiments .

3. Conclusions

We have developed a percolative approach to study the soft breakdown and the associated current fluctuations of ultrathin insulating films.The major results we found are summarized as follows. (i) The current evolution exhibits a pre-breakdown region. (ii) A threshold voltage depending on the resistance of the otherwise perfect film separates breakdown and steady state conditions. (iii) The spectrum of the current fluctuations is basically of Lorentzian type with a correlation time which is practically independent of the stress voltage. (iv) The variance of resistance fluctuations exhibits a giant super-quadratic enhancement at increasing the stress voltage. (v) At high enough voltage stress the current fluctuations exhibit a non-Gaussian behavior. These features, particularly the super-quadratic and the non-Gaussian behaviors of the current noise at high stress, are of relevant interest as failure precursors. These results are in satisfactory qualitative agreement with experimental findings.

Acknowledgements

This research is performed within the STATE project of INFM. Partial support is also provided by CNR MADESS II project.

References

 N. Vandewalle et al., Appl. Phys. Lett. 74 (1999) 1579.

 M. Houssa et al., Appl. Phys. Lett. 73 (1998) 514.

 C. Pennetta et al., IEEE Trans. on Electron. Dev. 47 (2000) 1986.

 G. Alers et al., Appl. Phys. Lett. 69 (1996) 2885.

 D.J. DiMaria and J.H. Stathis, Appl. Phys. Lett. 70 (1997) 2708

 C. Pennetta et al., Phys. Rev. Lett. 85, (2000) 5238.

Figure Captions

Figure 1 - Leakage current evolutions in a degrading dielectric kept under constant voltage for stress voltages (from bottom to top) of: $`0.1`$, $`0.3`$, $`0.5`$, $`0.6`$, $`0.7`$, $`0.8`$, $`0.85`$ and $`0.9V`$.

Figure 2 - Breakdown voltage $`V_b`$ as a function of the perfect network resistance $`R_{per}`$. The graph shows $`V_b\sqrt{R_{per}}`$.

Figure 3 - Power spectral density of current fluctuations of the steady-state dielectric film kept under constant voltage of (from bottom to top): $`0.1`$, $`0.3`$, $`0.5`$, $`0.6`$, $`0.7`$ and $`0.85V`$.

Figure 4 - Normalized variance of the network resistance fluctuations as a function of the applied voltage.

Figure 5 - Gaussianity check of current fluctuations for a stress voltage of $`0.5V`$ (open squares) and of $`0.85V`$ (filled circles).