text:
# References

$`\theta `$-Vacuum Systems Via Real Action Simulations

V. Azcoitia, G. Di Carlob, A. Galantec,b and V. Lalienaa

a Departamento de Física Teórica, Facultad de Ciencias, Universidad de Zaragoza,

50009 Zaragoza (Spain).

b Istituto Nazionale di Fisica Nucleare, Laboratori Nazionali del Gran Sasso,

67010 Assergi (L’Aquila), (Italy).

c Dipartimento di Fisica dell’Università di L’Aquila, 67100 L’Aquila, (Italy).

ABSTRACT

Inspired by the results of the Ising model within an imaginary external magnetic field, we introduce a transformation in quantum systems with a $`\theta `$-vacuum term that amounts to a rescaling of $`z=\mathrm{cos}(\frac{\theta }{2})`$. Making use of this transformation we are able to determine the order parameter as a function of $`\theta `$. The approach is successfully tested in models with both broken and unbroken $`CP`$ symmetry at $`\theta =\pi `$.

Quantum systems with complex euclidean actions have received an important amount of attention in recent time. Indeed there are two very relevant subjects in quantum field theory the analysis of which requires the use of complex actions: the study of topological structures and the behavior of matter at high temperature and density. Notwithstanding the big interest in these subjects, it is currently not possible to write an efficient code to perform numerical simulations of systems with complex actions. All efficient algorithms to simulate statistical systems are based in the use of a true probability distribution function as Boltzmann weight. This is the reason why all attempts to analyze systems with complex actions are based $`(i)`$ on the use of specific properties of the particular model under study , $`(ii)`$ uncontrolled analytical continuations, or $`(iii)`$ in a very delicate process of summing up all the contributions which appear in the definition of the partition function , , , , .

In a recent paper  we proposed a new approach to analyze systems with a $`\theta `$-vacuum term, tested it in two analytically solvable models, and applied to the study of $`\text{CP}^3`$ with a topological term in the strong coupling region. This approach is based on the use of a high accuracy method to measure the probability distribution function of the topological charge density at $`\theta =0`$ from numerical simulations at imaginary $`\theta `$. It is combined with a multiprecision algorithm which allows to sum up all the terms which appear in the partition function, that differ by many orders of magnitude, running from 1 to $`e^V`$, where $`V`$ is the lattice volume.

Even if the results reported in  are quite promising, the method is not free from systematics and it would be therefore important to have alternative approaches available in order to compare physical predictions. This is indeed the purpose of present paper.

We should say from the beginning that contrary to the approach reported in  where no assumption on the phase structure of the system was made, here we will make the assumption that the system has no phase transition at real $`\theta `$ except at most at $`\theta =\pi `$. This means that our proposal should be used in practice as complementary to other approaches and one should look for consistency between their results.

The method is based in extrapolating a suitably defined function to the origin. The function turns out to be very smooth in all the cases we have considered up to now and this makes us confident on the whole procedure. At the very end it is equivalent to an analytical continuation from imaginary to real $`\theta `$. However it has not to be confused with analytical continuation of the free energy density and/or the order parameter which are extremely sensitive to statistical errors and to the choice of fitting functions and therefore with small or no practical utility.

To fix the ideas and to keep this paper selfcontained, let us write the general form of the partition function of a system with a topological $`\theta `$-term as a sum over all topological sectors of the partition functions for each sector weighted with the proper topological weight,

$$𝒵_V(\theta )=\underset{n}{}p_V(n)\mathrm{e}^{\mathrm{i}\theta n},$$
(1)
where $`p_V(n)`$ is, up to a normalization constant, the probability of the topological sector $`n`$ at $`\theta =0`$, and the sum runs over all integers $`n`$. Writing equation (1) one needs to assume that $`CP`$ symmetry is realized at $`\theta =0`$ since otherwise $`Z_V(\theta )`$ would be ill-defined .

In almost all practical cases the sum in (1) has a number of terms of order $`V`$ since the maximum value of the topological charge at finite volume is of this order. The partition function (1) can also be written as a sum over the density of topological charge $`x_n=n/V`$, and setting $`p_V(n)=\mathrm{exp}[Vf_V(x_n)]`$ one gets:

$$𝒵_V(\theta )=\underset{x_n}{}e^{Vf_V(x_n)}\mathrm{e}^{\mathrm{i}\theta Vx_n}.$$
(2)
Equation (2) defines a $`2\pi `$ periodic function of $`\theta `$. For models $`CP`$-invariant at $`\theta =0`$ $`f_V(x_n)`$ is an even function; this property allow us to write:

$$𝒵_V(\theta )=\underset{y_n}{}G_V(y_n)\left(\mathrm{cos}^2\frac{\theta }{2}\right)^{Vy_n}$$
(3)
where $`y_n`$, which is non negative by construction, takes the same non negative values as $`x_n`$ and $`G_V(y_n)`$ is a functional of $`f_V(x_n)`$. The mean value $`y(\mathrm{cos}\frac{\theta }{2})`$ of $`y_n`$ computed from (3) is related to the mean value $`x(\theta )`$ of $`i`$ times the density of topological charge as follows:

$$x(\theta )=iy\left(\mathrm{cos}\frac{\theta }{2}\right)\mathrm{tan}\frac{\theta }{2}$$
(4)
which shows explicitly the $`2\pi `$ periodicity. At imaginary $`\theta `$, $`\theta =ih`$, the previous equation becomes

$$x(ih)=y\left(\mathrm{cosh}\frac{h}{2}\right)\mathrm{tanh}\frac{h}{2}$$
(5)
Furthermore, at real $`h`$, the system described by (1) can easily be simulated by standard numerical algorithms since in this case the Boltzmann weight is well defined. From numerical simulations at real $`h`$ one can get $`y(z)`$ (5) for values of the argument $`z=\mathrm{cosh}\frac{h}{2}`$ between 1 and $`\mathrm{}`$. In order to get the topological charge density at real $`\theta `$ one should extend this function to the $`(0,1)`$ interval in $`z`$. But as stated before, an extension of $`y(z)`$ is sensitive to statistical errors and fitting functions so that it makes this procedure not reliable.

The function $`y(z)`$ depends also on the couplings of the model under consideration. For the one-dimensional Ising model within an external real or imaginary magnetic field $`h`$, the function $`y(z)`$ can be analytically computed and the final result is

$$y(z)=\frac{\left(e^{4F}1\right)^{1/2}z}{\left(1+\left(e^{4F}1\right)^1z^2\right)^{1/2}}$$
(6)
which shows explicitly how this model, at imaginary values of the magnetic field $`(0<z<1)`$ and antiferromagnetic couplings $`\left(F=\frac{J}{KT}<0\right)`$ breaks spontaneously the $`Z(2)`$ symmetry at $`h=i\pi `$ (see (4)). But this model shows also a very interesting property: a change in the coupling $`F`$, or what is the same, a change in the physical temperature $`T`$, is equivalent to a rescaling of the variable $`z`$ since $`y`$ depends on $`z`$ through the combination $`e^{\frac{\lambda }{2}}z`$, with $`\lambda =\mathrm{log}(e^{4F}1)`$. In other words, if we define the transformation:

$$y_\lambda (z)=y\left(e^{\frac{\lambda }{2}}z\right)$$
(7)
which would allow us, for negative values of $`\lambda `$, to extend the order parameter from imaginary $`\theta `$ to real $`\theta `$, the $`1d`$ Ising model shows the exceptional feature that the transformation (7) is equivalent to a change of the coupling. This property, if true for other models, would be extremely interesting since it would reduce the problem of complex actions to simulations with real actions. Unfortunately this seems not to be the case since it is easy to check that already compact $`U(1)`$ in two dimensions with a $`\theta `$-term does not exhibit such a property. However this is not the end of the story and we will try in what follows to make use of transformation (7) in a much less ambitious way.

Many of the physically relevant models, as $`\text{CP}^{N1}`$ or $`QCD`$, are asymptotically free. Asymptotic freedom implies that the continuum limit is reached in the weak coupling region, but in this region the order parameter (and $`y(z)`$) will be very small for values of $`z`$ near 1 since the density of topological structures in lattice units is very suppressed at weak coupling. This means that, from simulations at real $`h`$, we can plot $`y_\lambda /y`$ against $`y`$ up to values of $`y`$ very close to the origin, provided that $`y(z)`$ does not vanish for $`z>0`$. We cannot exclude this possibility completely, but it does not happen in any of the exactly solvable models we know. Indeed, it is easy to show that $`y(z)0`$ for $`z>0`$ if all the coefficients $`G_V(y_n)`$ entering the right hand side of equation (3) are positive, as it happens in the models studied in the present work. In such a case, the fluctuation-dissipation theorem ensures that $`y(z)`$ is a non decreasing function of $`z`$. Furthermore, if the order parameter does not diverge at $`\theta =\pi `$, $`y(z)`$ should vanish at least linearly with $`z`$ at $`z=0`$ (that corresponds to $`\theta =\pi `$).

Moreover, if $`y_\lambda /y`$ as a function of $`y`$ shows a smooth behavior near $`y=0`$, we can expect a simple extrapolation to $`y=0`$ being reliable. This point is crucial in what follows and obviously can not be taken for granted. Apart from being reasonable it has been verified in all the models whose analytical solution is known to us. In particular it is true for the one-dimensional Ising model, compact $`U(1)`$ in two dimensions, $`\text{CP}^3`$, and a toy model that exhibits no phase transition at $`\theta =\pi `$, thus covering all the possible different qualitative possibilities.

In the Ising case we have done simulations of the model within an external real magnetic field $`h`$ at $`F=2.0`$. This coupling was chosen in order to be in the ”weak coupling region” where $`y(z)`$ is small near $`z=1`$. The model breaks spontaneously the $`Z(2)`$ symmetry at $`\theta =ih=\pi `$ and therefore $`y(z)`$ vanishes linearly with $`z`$ (see (4), (5)).

If we define an effective exponent $`\gamma _\lambda `$ by the following equation:

$$\gamma _\lambda =\frac{2}{\lambda }\mathrm{log}\left(\frac{y_\lambda }{y}\right)$$
(8)
this exponent, in the limit $`y(z)0`$ ($`z0`$), will give the dominant power of $`y(z)`$ as a function of $`z`$ near $`z=0`$. A value for this exponent that goes to 1 when $`y0`$ will imply spontaneous symmetry breaking, and values between 1 and 2 will signal a second order phase transition at $`\theta =\pi `$ with a divergent susceptibility. The reason is that the order parameter, $`x=\mathrm{tan}(\theta /2)y[cos(\theta /2)]`$, will behave as $`(\pi \theta )^{\gamma _\lambda 1}`$ for $`\theta \pi `$.

Fig. 1 reports $`y_\lambda /y`$ against $`y`$ for the Ising model at $`\lambda =0.5`$. The solid line is a quadratic fit of the points. This figure contains also the results for the effective exponent $`\gamma _\lambda `$ against $`y`$ (lower points). A simple quadratic fit predicts $`\gamma _\lambda (0)=0.997(3)`$ i.e., spontaneous symmetry breaking. We have checked that, as expected, we get compatible results using different values of $`\lambda =O(1)`$, and the same is true for the rest of the models considered.

The fit of $`y_\lambda /y`$ against $`y`$ together with the dependence of $`y(z)`$ on $`z`$ for $`z1`$ (computed from numerical simulations at real $`h`$) allows one to reconstruct the order parameter (4) at real $`\theta `$. To simplify the discussion let us assume we have computed $`y(z_0)=y^0`$. Next using the fit of Fig. 1 we plot $`y_\lambda `$ against $`y`$ and look in this plot for the value of $`y=y^1`$ which corresponds to $`y_\lambda =y^0`$ (see Fig. 2). It is simple to verify that $`y^1=y(e^{\lambda /2}z_0)`$. The next step is to iterate this procedure: we replace $`y^0`$ by $`y^1`$ and get $`y^2=y(e^\lambda z_0)`$, and so on. Once we have computed in this way $`y(z)`$ for as many values of $`z1`$ as desired, the order parameter can be computed using equation (4).

Fig. 3 contains the results for the Ising order parameter. We report also in this figure the exact analytical results in order to put in evidence the agreement.

For two-dimensional compact $`U(1)`$ with topological charge we chose to work at $`\beta =1.0`$ and $`\lambda =0.5`$. Figures similar to 1 and 3 were obtained for this model. The agreement between numerical and analytical results for the order parameter was even better than for the Ising model. Because of space limitations the corresponding figures are not included.

The next model we want to analyze here is $`\text{CP}^3`$ with topological charge. Contrary to the Ising and two-dimensional compact $`U(1)`$ cases, in which we can compare with the exact analytical results, $`\text{CP}^{N1}`$ models, which share many interesting properties with $`SU(N)`$, are not analytically solvable.

We have adopted for the action the “auxiliary U(1) field” formulation:

$$S_g=N\beta \underset{n,\mu }{}(\overline{z}_{n+\mu }z_nU_{n,\mu }+\overline{z}_nz_{n+\mu }\overline{U}_{n,\mu }2)$$
(9)
where $`z_n`$ is a $`N`$-component complex scalar field that satisfies $`\overline{z}_nz_n=1`$ and $`U_{n,\mu }`$is a $`U(1)`$ “gauge field”. The topological charge operator can be constructed directly from the $`U(1)`$ field:

$$S_\theta =i\frac{\theta }{2\pi }\underset{p}{}\mathrm{log}(U_p)$$
(10)
where $`U_p`$ is the product of the $`U(1)`$ field around the plaquette and $`\pi <\mathrm{log}(U_p)\pi `$.

Simulations were done at $`\beta =0.4`$. This value of the coupling is not in the weak-coupling scaling region, but the values of $`y(z)`$ near $`z=1`$ are smooth enough to allow a simple extension of $`y_\lambda /y`$ until $`y=0`$. Fig. 4 contains the effective exponent $`\gamma _\lambda `$ against $`y`$ at $`\lambda =0.5`$. A simple quadratic fit predicts $`\gamma (0)=1.03(5)`$ i.e., spontaneous symmetry breaking at $`\theta =\pi `$. By imposing, as suggested by the previous fit, $`\gamma (0)=1.0`$ in the fit of $`y_\lambda /y`$ against $`y`$, we get for the order parameter the results reported in Fig. 5. The continuous line in this figure represents the numerical results obtained within the approach reported in . The very good agreement between the two approaches makes quite reliable the numerical results and $`CP`$ structure at $`\theta =\pi `$ of Fig. 5.

The models previously analyzed have a common feature: $`CP`$ symmetry is always spontaneously broken at $`\theta =\pi `$. In order to check that our approach works also in the symmetric case we have analyzed a simple toy model, with only one effective degree of freedom, where $`CP`$ symmetry is enforced at $`\theta =\pi `$. This model, which resembles very much the free instanton gas model, is described at finite volume $`V`$ by the following partition function:

$$𝒵_V(\theta )=\left(1+A\mathrm{cos}\theta \right)^V.$$
(11)
The order parameter at real $`\theta `$

$$x(\theta )=i\frac{A\mathrm{sin}\theta }{1+A\mathrm{cos}\theta }$$
(12)
shows explicitly $`CP`$ symmetry at $`\theta =\pi `$, as it should be since this is a model of one degree of freedom. We have not done for obvious reasons numerical simulations of this model at real $`h=i\theta `$ but in order to work in realistic conditions, we modified the exact value of the order parameter at real $`h`$ by adding to it a random gaussian error of 1 percent. Fig. 6 contains the exact analytical results for the order parameter at $`A=0.005`$ and the points obtained within our approach, using always for the parameter $`\lambda `$ defined in (7) the value 0.5. The agreement is quite good.

A final comment on the possibility of applying this approach to the inverse problem i.e., finite density $`QCD`$, is worthwhile. $`QCD`$ at imaginary chemical potential and staggered fermions can be simulated numerically because the fermion determinant is positive definite. Its partition function belongs to the general class of partition functions described by (2), where $`x_n`$ now stands for the density of baryonic charge and $`\mathrm{exp}[Vf_V(x_n)]`$ is the probability distribution function of the baryon density at vanishing chemical potential. From numerical simulation of QCD at imaginary chemical potential one can measure the mean value of the density of baryonic charge and from it we can get $`y(z)`$ for $`z`$-values between 0 and 1. In order to get the density of baryonic charge at real chemical potential we should extend $`y(z)`$ to $`z`$-values larger than 1. Summarizing we know $`y_\lambda /y`$ as a function of $`y`$ around $`y=0`$ and we need to extend this quantity to larger values of $`y`$. But since a phase transition is expected in QCD at finite chemical potential, a singular behavior of $`y_\lambda /y`$ as function of $`y`$ should also be expected, and any analytical extension of this quantity should fail when trying to reproduce the correct behavior.

This work has been partially supported by an INFN-CICyT collaboration and MCYT (Spain), grant FPA2000-1252. Victor Laliena has been supported by Ministerio de Ciencia y Tecnología (Spain) under the Ramón y Cajal program.


format: text

text:
# The post-AGB evolution of AGB mass loss variations

## 1 Introduction

High resolution imaging with the Hubble Space Telescope HST have revealed the presence of so-called ‘rings’ around Planetary Nebulae (PNe) and proto-PNe. Examples can be found around the proto-PNe CRL 2688 (Sahai et al. 1998, Egg Nebula), IRAS 17441-2411 (Su et al. 1998), IRAS 17150-3224 (Kwok et al. 1998), IRAS 16594-4656, IRAS 20028+3910 (Hrivnak et al. 2001), and the PNe NGC 7027, (where they show up in scattered light), NGC 6543 and Hb 5 (Terzian & Hajian 2000). Inspection of HST archival images of NGC 3918 also reveal ring-like structures in its halo, confirmed by ground-based data (Corradi et al. 2003). Similar rings are observed around one AGB star, the extreme Carbon star IRC+10216 (Mauron & Huggins 2000).

The name ‘rings’ is somewhat of a misnomer since these are three-dimensional structures. However, the term is now commonly used and we will adher to this convention. Although never completely closed, the rings show an overall spherical appearance, and seem to be associated with a time when the mass loss from the AGB star was roughly spherical. Interestingly, they tend to sit right around the aspherical nebular structures, suggesting they are connected to the phase just before the onset of asphericity in the system.

Their origin has been attributed to three different processes. Mastrodemos & Morris (1999) showed that binary interaction can produce ring-like structures, both in and perpendicular to the orbital plane. In this case the rings are actually the spiral waves produced by the presence of the companion star. García-Segura et al. (2001) reproduce the ring features by invoking a magnetic cycle, similar to that in the Sun, for the AGB star. The period is then set by the period of this magnetic activity cycle, which is a parameter in their model. Frank (private comm) has pointed out that any magnetic dynamo cycle in AGB stars would have a much shorter period than the one corresponding to the rings. The third model comes out of a study of the dust-driven winds of AGB stars. Simis et al. (2001) showed that intricate coupling between dust and gas, leads to cycles in which predominantly either large or small grains are produced, which then produce periodic variations in the wind density and velocity. The period depends on the sound travel time of the subsonic wind region, and is found to lie close to the observed values for the rings.

In case the rings are ionized it is possible to use plasma diagnostics to derive the temperature and density in the rings. To date this has not been done in detail. Balick et al. (2001) estimated some properties of the rings around NGC 6543 on the basis of HST images and limited spectroscopic data. Hyung et al. (2001) used HST images taken with WFPC2 in narrow band filters near the lines of \[O III\] 5007Å and 4363Å to find the electron temperature in the rings of the same object. Although Balick et al. (2001) found no indication for unusual electron temperatures in the rings, Hyung et al. (2001) claim to find substantially higher electron temperatures, indicating the presence of shocks. Both methods have their limitations, and the result remains inconclusive. Balick et al. (2001) derived density variations from H$`\alpha `$ emission variations, and found them to be consistent with individual shells with negligible amounts of material between them.

Both Bryce et al. (1992) and Balick et al. (2001) report that emission lines in the rings area of NGC 6543 are unusually broad, with FWHM of around 30 km s-1, for which they could offer no explanation. Hyung et al. (2001) presented a simple model for explaining the putative higher electron temperatures, as well as the broadened emission lines. Following the evolution of sinusoidal density and velocity variations using a numerical hydrodynamic model including the effects of ionization, they could reproduce temperatures and broadened lines. However, the model only included the rings, not the core nebula, nor did it follow the full evolution of the rings, and in order to keep the velocity field under control, unreasonably high velocities had to be used.

Studying the evolution of the rings from the AGB to the PN phase is interesting for several reasons. Firstly, one may be able to find differences between the three proposed models. Secondly, it can help us understand whether the limited number of PNe found to have rings is intrinsic, e.g. because ionization changes the properties of the rings, or rather due to lack of high resolution observations. Also, it will provide insight in the connection between the AGB mass loss variations and the ionized PN rings.

In this paper we study the evolution of mass loss variations on the AGB through the proto-PN phase to the PN phase. For this we will use an improved version of the model of Simis et al. (2001) to produce the AGB mass loss variations, and an improved version of the model used in Mellema (1994) to follow their evolution during the post-AGB phase. Although we pick one model for the origin of the mass loss variations, we will show that most of the results have relevance for the other models too.

The layout of the paper is as follows. In Sect. 2 we describe the numerical methods used to produce the mass loss variations and to study their fate from the AGB through the post-AGB phase. Section 3 presents the results of a typical simulation, for which we derive observables in Sect. 4. In Sect. 5 we discuss the results and how they increase our understanding of the origin and evolution of the rings. Our conclusions are summed up in Sect. 6.

## 2 Two models

In order to model the evolution of AGB mass loss variations we first need to produce them. To keep the model as physical as possible, we produce the mass loss fluctuations using a slightly improved version of the two-fluid dust-gas code described by Simis et al. (2001). These are then used as input for a gas-only hydrodynamics code which can treat the detailed ionization, heating and cooling processes in step with the evolution of the star. We will now describe these codes in some more detail.

### 2.1 A dust-driven AGB wind model

We provide here a short description of the two-fluid hydrodynamics code which we use to model the AGB. For a more detailed description of the code and the implementation of the two-fluid flow in particular, the reader is referred to Simis et al. (2001, 2003). The new aspect of this code is the fact that it combines time-dependent hydrodynamics with two-fluid flow and a self consistent description of the dust (including nucleation and growth). It turns out that the interaction of these physical mechanisms is critical for the formation of the shells, which are a result of alternating low- and high-mass loss episodes of outflow. The mass loss rate is low when the momentum transfer from grains to gas is inefficient. This is the case when the gas has a relatively low velocity while the dust grains move relatively fast (i.e. they have a high drift velocity). This results in a rapid passage of the grains through the zone in which grain growth is efficient, and hence in somewhat smaller grains. As a consequence, the drag force exerted on the gas is relatively weak. This leads to a decrease of the gas density and hence an even larger drift velocity, because the frequency of gas-grain collisions decreases. So, during the low mass-loss phase the average grain size gradually decreases whereas its drift speed increases. In a phase of high mass loss, the rate of momentum transfer from the radiatively accelerated grains to the gas is high. This results in a large velocity of the gas and a low drift velocity. Because the grains move relatively slow through the gas they have sufficient time to grow. The average grain size during this phase is relatively large, which causes an efficient momentum transfer. This enables a high gas density, a gradual growth of the grains and a further decrease of their drift speed. The transition between the two phases is triggered by the sudden increase of the drag force at the end of the low mass loss phase, which is a result of the quadratic appearance of the (increasing) drift velocity in the drag force. The sudden increase of the drag force leads to an increasing momentum transfer per collision and it prevents the ongoing decrease of the collision rate that was a result of the gradually decreasing grain size. This starts the high mass loss phase. Curiously, it also triggers the onset of the next low mass loss phase: due to the suddenness of the transition, an outward shock develops. At the same time, a rarefaction wave moves towards the stellar surface. This wave will dilute the gas, thereby again creating a condition for inefficient grain growth, small grains, low collisional cross sections and hence a low mass loss rate. Hence, the characteristic time scale for the mass loss fluctuations is the time used by the rarefaction wave to cross the dust forming part of the subsonic region of the outflow.

The code is a time dependent, Eulerian, explicit hydrodynamics code, based on the FCT/LCD algorithm assuming spherical symmetry (Boris 1976; Icke 1991). Physical processes taken into account in the code involve (equilibrium) gas chemistry, nucleation and growth of grains (Gail et al. 1984; Gail & Sedlmayr 1988; Dorfi & Höfner 1991), semi-analytical radiative transfer (Fleischer et al. 1992) and the behaviour of grains in collision with gas particles. The temperature is calculated assuming local radiative equilibrium, following Lucy (1976). Grain-gas collision is implemented to provide a realistic description of the momentum transfer between the two fluids. This enables a full two-fluid description, i.e. no assumptions about the motion of the grains relative to the gas need to be made. The extra degree of freedom introduced thereby turns out to be essential for the formation of quasi periodic shells: when comparing two otherwise identical model calculations, the (two-fluid) model in which grains can drift with respect to the gas does produce the shells, whereas the single fluid model, i.e. the model in which grains are always assumed to have the same velocity as the gas, has a more or less stationary outflow.

### 2.2 A post-AGB model

The numerical method we use to follow the subsequent evolution in the post-AGB phase, descends from the method used by Mellema (1994). This code solves the hydrodynamic equations in one-dimension, using the spherical radial coordinate. The solver of choice is the ‘Roe solver’ (Roe 1981; Eulderink & Mellema 1995), a second-order accurate, approximate Riemann solver.

The atomic physics is treated according to the approach described in Mellema & Lundqvist (2002) (DORIC), but with the addition of photo-ionization from an evolving source. This means that we follow the time-dependent ionization of all ions of H, He, C, N, O, and Ne, taking into account collisional ionization, photo-ionization, radiative and dielectronic recombination and charge exchange with hydrogen. The local density and temperature of the ions is used to calculate the local non-equilibrium cooling rate. The heating rate is given by photo-ionization heating of H and He.

The transfer of ionizing radiation is treated according to the method described in Frank & Mellema (1994) and Mellema et al. (1998), i.e. using the three frequencies corresponding to the ionization thresholds of H0, He0, and He+, but taking the frequency dependence of the cross sections into account when calculating the integrals over the spectral energy distribution. At all times, the star is assumed to radiate as a black body with the effective temperature given by the stellar evolution model used (see below).

The outer boundary is an outflow condition, the inner boundary condition is provided by the evolving stellar wind in the post-AGB phase.

#### 2.2.1 The post-AGB wind

During the post-AGB phase a star can reach temperatures over 105 K. The bulk of the radiation is emitted in the UV, where the outer layers have an abundance of absorption lines. Castor et al. (1975) found that for hot stars, absorption lines dominate in transferring momentum from the radiation to the gas, and the consequent stellar wind is therefore called ‘line-driven’. The theory of line-driven winds is well developed and tested. Here we use the analytical approximation derived by Kudritzki et al. (1989), which for a given effective temperature and luminosity gives the stellar wind mass loss rate and terminal velocity.

The mass loss rate of a line-driven wind is proportional to the luminosity and the number of absorbing lines. The acceleration of the atoms is proportional to a so called force multiplier, which is approximated by a simple function:

$$M(t)=kt^\alpha (10^{11}n_e/W)^\delta ,$$
(1)
where $`t`$ is a dimensionless optical depth parameter, $`k`$, $`\alpha `$ and $`\delta `$ force multiplier parameters and $`W`$ the dilution factor. The acceleration also depends on the velocity field in the wind as a function of radius. This velocity field is given by:

$$v(r)=v_{\mathrm{}}(1\frac{R_{}}{r})^\beta $$
(2)
Apart from the stellar luminosity and effective temperature, the model also needs values for these force multiplier parameters, $`k`$, $`\alpha `$, $`\beta `$ and $`\delta `$, as well as for the stellar mass, the helium abundance, and the average number of electrons produced by helium.

The luminosity, mass and effective temperature are taken from Blöcker (1995), who presented a compilation of evolutionary tracks for post-AGB masses in the range 0.53 to 0.94 M. Three of these tracks were originally calculated by Schönberner (1983). The force multiplier parameters we take from the compilation in Lamers & Cassinelli (1999, Table 8.2).

With these parameters the wind depends strongly on the stellar effective temperature. When the star just enters the post-AGB phase, its effective temperature is low, and the mass loss rate and wind velocity are low too. As the effective temperature goes up, the mass loss rate increases, levelling off to a roughly constant rate from about 10000 K. At this point the number of effective lines reaches a maximum. When nuclear burning stops, the luminosity of the star starts to drop and the mass loss rate goes down, but the velocity of the wind keeps rising. Figure 1 shows the evolution of the wind as a function of stellar effective temperature for four different evolutionary tracks.

## 3 Simulations

We performed five simulations for different initial density distributions and stellar masses. In this section we describe the results of one and refer to the results of the other simulations only when they help us understand certain features of behaviour.

### 3.1 The initial density distribution

In Sect. 2.1 the mechanism of the mass loss variations during the AGB was discussed. From the two-fluid dust-gas simulations of Simis et al. (2001) we took one result as the input for our simulation. These mass loss and velocities variations were calculated using the properties of IRC+10216 (Winters et al. 1994): a stellar mass of 0.7 $`M_{}`$, an effective temperature of $`1.98\times 10^3`$ K, luminosity $`2.4\times 10^4`$ $`L_{}`$, and a carbon/oxygen ratio of 1.4.

Our main aim is to follow the post-AGB evolution of the type of variations in density and velocity which are being produced by the two-fluid calculations. We therefore scale the outcome of the two-fluid simulation to obtain a more typical AGB wind with an average mass loss rate of $`2.510^5M_{}`$ yr-1, and an average velocity of 10 km s-1. This means that the AGB and post-AGB simulations are no longer physically connected and the parameters used for the AGB star listed above are no longer relevant for the results. What we are effectively doing is using the two-fluid calculations to generate a spectrum of density and velocity variations which at least has some physical basis (rather than assuming sinusoidal variations), and apply this to our AGB wind. This is also motivated by the fact that the quantitative outcome of the two-fluid simulation is quantitatively uncertain since the model does not include effects such as pulsation, evolution of the AGB star, etc.

The computational grid of the dust-gas simulations extended to a radius of only $`1.810^{16}`$ cm (1200 AU) and to be able to follow the long term (post-AGB) evolution we first had to fill a larger space with mass loss and velocitiy fluctuations. For this we extracted the density and velocity as a function of time at a radius of $`5.010^{15}`$ cm (333 AU), and used this as an inflow condition for a period of 31000 years, by which time the fluctuations fill most of the grid. The temperature at the inner boundary is taken from the two-fluid calculation, and we extrapolate outward using a power law, $`Tr^{2(\gamma 1)}`$, with an imposed minimum value of 10 K. The actual choice for the temperature profile is not very important for the dynamics, since the gas is flowing out supersonically.

As the fluctuations move out, the distances between minima in density become larger. This is due to the merging of shells, since these do not have the same velocity. Figure 2 shows the structure of the AGB wind 31000 years after we introduced the variable wind. One sees that the spacing of the rings increases with distance from the star. When this density pattern is projected on the sky the highest frequency structures disappear. As a consequence it is easier to distinguish the shells.

### 3.2 The dynamical evolution of the rings

Our standard model is for a star with a Zero Age Main Sequence (ZAMS) mass of 1.0 $`M_{}`$ and post-AGB mass of 0.565 $`M_{}`$. We call $`t=0`$ the start of this simulation, which corresponds to the start of the post-AGB evolution. The simulation follows the evolution of the circumstellar material over a time span of 15000 years and is performed on a grid of 6000 points. The post-AGB wind first collides with the slow AGB wind at a distance of $`5\times 10^{15}`$ cm, the inner edge of the grid lies at $`10^{14}`$ cm (6.6 AU).

In Fig. 3 the result after 400 years is shown. The effective temperature of the star is about 6000 K. There are enough energetic photons to ionize the fast wind. Because the velocity of the fast wind is 50 km/s and the mass loss rate of the star $`1.610^{12}M_{}`$ yr-1, the temperature in the shock, moving into the fast wind (the inner shock), is not that high, $``$ 3000 K. The pressure in the inner shock is actually lower than in the slow wind. As a consequence, the slow wind dissipates slightly inwards. The ionization front lies at the interaction between the fast and the slow wind.

Sixhunderd years later (Fig. 4), the effective temperature is 8200 K, and the mass loss rate went up by three orders of magnitude to $`10^9M_{}`$ yr-1, with a velocity to somewhat over 100 km/s. As a consequence, the temperature in the hot bubble is well over $`10^5`$ K. The pressure is now the same as in the slow wind. One can also see the ionization front moving into the slow wind. The temperature here is about 2500 K. The gas in front of the ionization front is not shocked yet and has a temperature around 100 K. One also sees that the shells merge while they move outward. At the place where the shells interact the temperature rises to a few hunderd degrees.

The ionization front has moved further out to 0.15 pc, 2100 years later (see Fig. 5), when the star’s effective temperature has risen to 20000 K. The part of the wind that is ionized has almost no density variations. Since the pressure of the gas in the ionized region is higher, the gas in front of the ionization front is pushed outward. This gas piles up in a thin shocked shell and with a temperature somewhat lower than the ionized region. This is typical of a D-type ionization front, where the speed of the front is slow enough for the gas to react to the increased pressure of the ionized region. This configuration also produces the typical positive velocity gradient in the ionized material, found in many PNe.

Almost 3700 years after the end of the AGB (see Fig. 6) the star’s effective temperature ($``$ 24000 K) reaches the point that the ionization front breaks through. In a few hundred years the whole nebula gets ionized and goes from ionization bounded to density bounded. The temperature in the ionized part of the slow wind rises to about 7000 K, and as a consequence the pressure also rises. The pressure being proportional to $`\rho T`$, follows the density profile and large pressure differences of an order of a magnitude occur. The gas tries to equalize this out, the effect of which can be seen in the velocity field. In the ionized region the velocity variations are over a shorter distance and have a larger amplitude than in the neutral region.

At 4500 years (see Fig. 7) the velocity variations have increased, and the velocity in the rings now range from 1 to 20 km s-1. These variations are much larger than in the period before ionization. To check the dependence on the initial velocity fluctuations, we did a simulation with the same initial density variations but with a constant velocity field. We found essentially the same velocity variations developing after ionization. This shows that the initial velocity variations are unimportant in this phase.

Considering the situation 3500 years later (see Fig. 8), we find that the density and velocity variations are still visible. These variations eventually equalize out, but even after 15000 years some traces of it can still be distinguished.

## 4 Synthesized observations

### 4.1 Emission lines and temperature

A radial density profile does not necessarily correspond to what would be observed. In order to compare to the observations, we calculated emissivities of certain lines, in \[O III\] 5007Å, the line in which the rings in NGC 6543 are most easily observed, as well as the \[O III\]4363Å line, whose ratio with the 5007 line gives us an ‘observed temperature’. We calculated the line emissivity from the density, temperature, and O2+ concentration at every point. The emission is calculated conform Sect. 3.5 and Tables 3.4 and 3.8 of Osterbrock (1989).

The radial \[O III\]5007Å and 4363Å emissivities at $`t`$=4500 years are shown in Fig. 9, together with the density and temperature. One can see that the density profile does not correspond to the \[O III\] profiles. This is due to the temperature dependence of the \[O III\] line which gives more weight to points of higher temperature, especially for the 4363Å line.

To compare to observations, we projected the emissivities onto the sky, and convolved them with a gaussian of FWHM corresponding to 0.1 arcsec at the distance of NGC 6543. We then took the ratio between the 5007 and 4363Å lines to determine the observed temperature. These projected emission and temperature profiles for the rings area are shown in Fig. 10.

The \[O III\] emission profile shows clear maxima and minima. The spacing between the maxima is about 0.02 pc. The dashed line shows the same in H$`\alpha `$. Clearly the two have different radial dependencies. Figure 2 of Balick et al. (2001) shows how the rings are much more pronounced in the \[O III\] profile than in the H$`\alpha `$ profile and how the H$`\alpha `$ emission falls off more rapidly, something which our model at least qualitatively reproduces.

The observed temperature varies between 7500 K and 10000 K, but most of this is due to a small outward gradient. Since the radial temperature profile does not show a temperature gradient, we must conclude that this is a projection effect. This is an interesting side-result, showing another way in which deriving electron temperatures from line ratios can be somewhat misleading. The true local variations are mostly of order 500 K. The biggest variations in temperature are seen in the outer part of the nebula, where most shocks occur at this time.

One thousand years later (Fig. 11), the temperature variations have decreased substantially. The largest variations are still seen in the outer region, where the shells are merging. Due to this merging the spacings between maxima have become larger.

At $`t`$=8000 years almost no temperature variations are seen, the amplitude has dropped below a few hunderd Kelvin (Fig. 12). Also in the density profile the rings have become less noticeable, although they can still be distinguished.

### 4.2 Lineprofiles

Combining the emissivities with the velocity information from the simulation, allows us to derive line shapes. We present line shapes at different positions along the rings in Fig. 13. Since the temperature of the gas is around 7500 K, the lines are also thermally broadened which we simulated by convolving with a gaussian of 6 km s-1 FWHM. The effects of seeing or slit width were not taken into account here.

One sees how the profiles are quite broad, up to 30 km s-1 FWHM. This is due to the fact that at a given position one observes several overlapping shells producing different velocity components. The profiles are clearly not gaussian, so the effect that the emission comes from different shells is visible.

These line profiles match the reported observations quite well. Bryce et al. (1992) and Balick et al. (2001) describe (but do not show) wide ($`30`$ km s -1) line profiles in the region of the rings. The widths match our model line profiles well, the low S/N of the spectroscopy of the faint rings makes further comparison difficult.

## 5 Discussion

The comparison with the observations of NGC 6543 in the previous section shows that our model matches a number of key observations:

1. The rings are clearly visible in the \[O III\] images.
2. The H$`\alpha `$ and \[O III\] emissivities have slightly different radial profiles.
3. The line profiles are broadened to $`30`$ km s -1.

This shows that our model gives a good description of the rings.

Comparing the spacing from our models with what is found in NGC 6543, we find a small descrepancy. Our spacing is 0.02 pc, about twice as large as in the observations (0.01–0.015 pc), assuming the distance to NGC 6543 is reliable. However, the exact spacing of the rings is influenced by many factors. First, the original spacing in the variable AGB mass loss depends on the density scale height at the base of the circumstellar envelope. Second, over time shells merge, increasing the spacing, and hence making it a function of time. A faster evolution will thus produce shorter distances between the shells.

The derived observable temperature (from the \[O III\]4363/5007 ratio) is found to vary somewhat, but not as much as claimed by Hyung et al. (2001), supporting the results of Balick et al. (2001) for an electron temperature comparable to those in the core nebula. We predict temperature variations of not much more than 1000 K, and perhaps a positive outward temperature gradient, depending somewhat on the phase NGC 6543 is in.

Our results show that the rings are a transient phenomenon, persisting not much longer than $``$3000 years after photo-ionization. This seems to be in line with the observed fact that most rings are found in their neutral state around proto-PNe or young PNe (e.g. NGC 7027). In this view the ionization in NGC 6543 is rather recent. In Hb 5 it is not known whether the rings are seen in scattered or intrinsic emission.

Another prediction from our models is that one will not find rings in so-called ‘attached haloes’ or shells around PNe, since these are thought to have been produced by the action of a D-type ionization front (Marten & Schönberner 1991; Mellema 1994), which as we saw, erases the density variations quite efficiently.

### 5.1 Origin of the rings

To produce our models we started with density and velocity variations produced through the mechanism suggested by Simis et al. (2001), which we slightly modified by taking a lower velocity, and scaling the density to obtain a certain average mass loss rate. Our simulations showed that for the structure during the PN phase, the initial velocity fluctuations are irrelevant, and so the important property is the nature of the density fluctuations. Studying the results of a number of simulations, we found that the key property is the amplitude of the density variations. At the start of the post-AGB phase, these have to be roughly one order of magnitude. Only then ionization leads to large enough pressure fluctuations to result in what could be loosely described as ‘supersonic turbulence’, which helps the longevity of the rings. Without the resulting shock waves the rings diffuse away in a few hundred years after ionization.

We would therefore expect that using density variations with similar amplitudes would give more or less the same results. In the model of García-Segura et al. (2001), for reasonable values of the magnetic field, the density initially varies less than in our simulations. However, since these shells are ‘frozen in’, they do not merge, and conserve their density contrast. At the time when the ionization front breaks through, the density contrast is likely to be similar to the one in our simulation. We therefore expect the rings from the MHD model to evolve similarly during the PN phase, and survive for about a few thousand years after ionization, producing broad line profiles as the flow is slowly erasing the rings.

Note that the magnetic pressure is able to support the rings only during the AGB and proto-PN phases. Upon ionization, the magnetic pressure becomes negligible in comparison to the gas pressure, and there is no real difference between the two models. This rather limits the usefulness of the magnetic field, since it is not clear whether it is really needed to maintain the fluctuations up to the time of ionization.

Binary interaction as proposed by Mastrodemos & Morris (1999) is probably not an explanation for the rings seen in planetary nebulae. The density contrasts they find in their models are much too low to survive for a long time after being ionized, and in fact already seem to be diffusing within a couple of orbital periods.

## 6 Conclusions

We numerically followed the evolution of AGB mass loss fluctuations during the post-AGB phase, paying special attention to the effects ionization has. We find that using mass loss fluctuations as produced by the mechanism described by Simis et al. (2001) easily survive up to the time ionization sets in for stars of ZAMS mass of 1 $`M_{}`$or larger. Upon ionization, the pressure fluctuations lead to shocks, which raises the average absolute velocity of the rings, and locally raises the temperature. After a couple of thousand years the rings will have mostly disappeared.

The ionized rings have the following observational characteristics:

* The rings are more pronounced and can be seen further out in \[O III\] then in H$`\alpha `$.
* Only slightly higher electron temperature ($`\pm 1000`$ K) than the rest of the PN, as measured by the \[O III\] line ratios.
* Broad line profiles due to the turbulent velocity field and the many overlapping shells.

We therefore predict that the electron temperature from accurate line ratios of ionized rings will show only a small difference with the electron temperature of the core nebula, showing the result from Hyung et al. (2001) to be spurious. We also predict that high signal-to-noise spectroscopy may reveal non-gaussian line shapes in the ring region. Lastly, according to our description, rings will only be found in either neutral or recently ionized haloes.

We found this type of evolution to be quite general, provided that the density fluctuations at the time of ionization are larger than a factor of 10. For lower values the rings diffuse away through gentle sound waves.

Using this requirement the wide binary model as published by Mastrodemos & Morris (1999) disqualifies as an explanation for the rings in ionized PNe. The MHD model of García-Segura et al. (2001) and the dust-driven wind instability from Simis et al. (2001) both fulfill the requirement. However, there is some dispute on the choice of period for the magnetic cycle of AGB stars, making the dust-driven wind instability model (with no imposed period), the more likely explanation.

###### Acknowledgements.

The research of GM has been made possible by a fellowship of the Royal Netherlands Academy of Arts and Sciences.


text:
# High temperature gate control of quantum well spin memory

## Abstract

Time-resolved optical measurements in (110)-oriented GaAs/AlGaAs quantum wells show a ten-fold increase of the spin-relaxation rate as a function of applied electric field from 20 to 80 kV cm-1 at 170 K and indicate a similar variation at 300 K, in agreement with calculations based on the Rashba effect. Spin relaxation is almost field-independent below 20 kV cm-1 reflecting quantum well interface asymmetry. The results indicate the achievability of voltage-gateable spin-memory time longer than 3 ns simultaneously with high electron mobility.

The longest possible spin memory and an ability to control the orientation or relaxation of non-equilibrium spin populations in semiconductor quantum wells (QWs) via an applied gate voltage will be the key to many spintronic applications. Recent experiments on (110)-oriented III-V QWs have demonstrated a predicted Oh99 (1, 2) dramatic increase of spin memory at room temperature, up to 20 ns in an $`n`$-type GaAs/AlGaAs QW Ada01 (3), but although various possible approaches exist, gate control has hitherto been more elusive. The dependence of carrier spin relaxation rates on gate-injected electron or hole concentration in QWs has been investigated by optical techniques Sn90 (4, 5). A second approach exploits the electric field dependence of the electron $`g`$-factor in parabolic potential wells to control spin precession in applied magnetic field Sal02 (6). Similarly gated ferromagnetism in transition-metal-doped III-V’s Oh00 (7) may allow control of spin orientation without an applied magnetic field. A third proposal is to manipulate the conduction band spin-splitting by an applied electric field through the Rashba effect Ra60 (8). This permits gate control of the exponential spin relaxation at high temperatures where electrons experience strong scattering, the so called collision-dominated regime, Lau02 (9) or even coherent spin reorientation under collision-free conditions at very low temperatures DD90 (10). For the collision-dominated case this last approach has the strong advantage of applicability at room temperature.

In this paper we report a time-resolved optical investigation of electron spin relaxation at high-temperatures in undoped (110)-oriented GaAs/AlGaAs QWs. We confirm the existence of a 100-fold increase of spin memory compared to (001)-oriented QWs at 300 K Oh99 (1) but we now demonstrate a ten-fold variation of the spin relaxation rate with the application of a modest electric field. This is in accord with theoretical expectations Lau02 (9) for the Rashba effect and strongly indicates the mechanism which currently limits the spin memory in (110)-oriented samples in zero field.

The spin relaxation in our experiments can be interpreted on the basis of the D’yakonov, Perel’ and Kachorovskii (DPK) DK86 (2, 11) mechanism of spin relaxation, as refined in a non-perturbative approach Lau02 (9, 12, 13), which dominates zero-field spin dynamics in GaAs/AlGaAs QW systems Fl02 (13, 14, 15, 16). In this model, spin reorientation is driven by precession of the individual electron spin vectors induced by the effective magnetic field represented by the conduction band spin splitting, which results from spin-orbit coupling and lack of inversion symmetry. The corresponding Larmor precession vector $`𝛀(𝐤)`$ varies in magnitude and direction according to the electron’s wavevector $`𝐤`$ so that scattering results in a randomly fluctuating precession vector DK86 (2, 11, 17). In the collision-dominated regime, $`|𝛀|\tau _p^{}1`$ with $`|𝛀|`$ the average precession frequency and $`\tau _p^{}`$ the momentum scattering time of an electron, spin relaxation is inhibited by scattering and has rate DK86 (2, 11, 13)

$$\tau _s^1=\mathrm{\Omega }^2\tau _p^{}.$$
(1)
The first factor in this expression is determined by the vector $`𝛀`$, which, in a QW, has three contributions, $`𝛀^{SIA}`$, $`𝛀^{BIA}`$ and $`𝛀^{NIA}`$ Fl02 (13). The natural interface asymmetry component $`𝛀^{NIA}`$ does not occur in GaAs/AlGaAs QWs and is not considered further here. The most interesting for gating spin relaxation is $`𝛀^{SIA}`$, the Rashba or structural inversion asymmetry (SIA) term. It may arise from some built-in asymmetry of the structure or be induced by an externally applied odd parity perturbation such as an electric field. For field $`𝐄_z`$ along the growth axis ($`z`$) it has the form ($`𝐄_z`$$`\times 𝐤`$) Ra60 (8, 13) and so is oriented in the QW plane for all in-plane $`𝐤`$. This induces precession of the electron spin away from the $`z`$-axis and therefore contributes a term to DPK spin relaxation along the axis, which can be varied with applied electric field Lau02 (9). However except in (110)-oriented QWs the effect will be small because $`𝛀^{BIA}`$, which is due to bulk inversion asymmetry (BIA) of the zinc blende structure and is only weakly field-dependant, is either wholly or partially in the QW plane and so swamps $`𝛀^{SIA}`$ Lau02 (9, 10, 17). For the special case of (110)-oriented QWs, illustrated in Fig. 1,

$`𝛀^{BIA}`$ lies approximately along the growth axis for all $`𝐤`$, causing only weak spin relaxation along the growth axis and thereby leading to greatly enhanced spin memory, as predicted by D’yakonov and Kachorovskii DK86 (2) and confirmed experimentally by Ohno et al. Oh99 (1). Spin relaxation along the growth axis therefore becomes more sensitive to $`𝛀^{SIA}`$ and should be significantly changed by the application of an external electric field Lau02 (9), an important possibility, which we investigate here.

Equation (1) also shows that the spin relaxation rate scales with $`\tau _p^{}`$ in a way analogous to motional narrowing. For undoped QWs, $`\tau _p^{}`$ is closely related to the momentum relaxation time $`\tau _p`$ which determines the electron mobility DK86 (2, 11, 13). This means that long spin memory and high mobility are, in general, mutually exclusive in low-doped $`n`$-type QWs but may be simultaneously achievable in the special case of (110)-oriented QWs.

Samples from three different wafers have been studied each grown on a semi-insulating (SI) GaAs substrate and each containing twenty 7.5 nm QWs separated by 12 nm barriers with aluminium fraction 0.4. Two wafers consisted only of undoped layers on (001)- and (110)-oriented substrates respectively. The third wafer was a $`pin`$ structure on a (110)-oriented substrate with the undoped QWs grown between two 0.1 $`\mu `$m undoped buffer layers of AlGaAs and with layers of 2$`\times 10^{18}`$ cm-3 $`n^+`$ doped GaAs and of 1.2$`\times 10^{18}`$ cm-3 $`p^+`$ doped AlGaAs respectively below and above. The orientation of the (110)-oriented substrates was accurate to $`\pm 0.5^{}`$. Measurements were carried out at 300 K on all three wafers as grown and the $`pin`$ wafer was also processed into 400 $`\mu `$m mesa devices (Fig. 2a) and used for measurements in variable electric field at 170 K. Although the sequence of layers in the samples does not have inversion symmetry, for flat-band conditions the QWs and barriers possess local inversion symmetry where the confined electron wavefunctions have significant density. The doped layers of the $`pin`$ structure result in a built-in electric field $`𝐄_z`$ of about 25 kV cm-1. Figure 2b shows the photoluminescence (PL) spectrum of one of the mesa structures for zero applied bias at 170 K. The width, at 10 K, of the excitonic recombination line was $``$6.5 meV in each of our (110)-oriented wafers and 1.8 meV in our (001)-oriented wafer. Since we expect similar fluctuations in QW width for the two sets of multi-quantum wells (MQWs), this suggests that the interfaces of the (110)-oriented QWs are less perfect than for (001) growth and could be further improved by variation of growth conditions.

Electron spin evolution was obtained from the change of intensity $`(\mathrm{\Delta }R)`$ and polarisation rotation $`(\mathrm{\Delta }\theta )`$ of optical probe pulses reflected at near normal incidence from the sample surface along the growth axis, at a variable delay following excitation by 10 times more intense, nearly collinear circularly polarised pump pulses DD90 (10, 14, 16). The pulses were of 2 ps duration from a mode-locked Ti-sapphire laser tuned to the peak of the first electron-heavy-hole PL transition (Fig. 2b) giving excitation density $``$109 cm-2. Figure 2c shows the $`I`$-$`V`$ characteristics of the device in the dark and illuminated with experimental laser intensity, $``$350 $`\mu `$W. Such low powers are necessary to avoid screening of the electric field in the MQW region. Since this incident intensity could generate a photocurrent up to 30 $`\mu `$A, the magnitude of the negative bias current shows that in fact a very small fraction, $``$2$`\%`$, of the photo-carriers is swept out by the bias but even this may be exaggerated because the rapid increase of current below $``$1.5 volts indicates avalanche multiplication. The absorbed pump photons thus generate ‘cold’ excitons which dissociate into free carriers on a sub-picosecond timescale with spins polarised along the growth axis and which remain confined in the QWs. The hole spins rapidly relax whereas the electron spin-relaxation is much slower Br98 (14, 16). On a timescale longer than $``$1 ps, $`\mathrm{\Delta }R`$ gives a measure of the population of photoexcited carriers and $`\mathrm{\Delta }\theta `$ a measure of the $`z`$-component of electron spin so that their ratio gives the pure longitudinal spin dynamics of the electrons Br98 (14, 16).

Figure 3 shows the observed $`\mathrm{\Delta }\theta `$ signals for the three as-grown wafers at 300 K. There is a dramatic difference between the decay for the (001)-oriented wafer (A) and the undoped (110)-oriented wafer (B), which essentially reproduces the findings of Ohno et al. Oh99 (1). However the decay for the $`pin`$ (110)-oriented wafer (C), where there is a built-in electric field of about 25 kV cm-1 (at 300 K), is significantly faster than for the undoped (110)-oriented wafer (B). Measurements were made at various pump intensities to allow extrapolation to zero power and when the decay of the $`\mathrm{\Delta }R`$ signal is also included in each case the spin relaxation times are found to be $`32\pm 1`$ ps, $`3.5\pm 0.2`$ ns and $`0.85\pm 0.02`$ ns for A, B and C wafers respectively. These figures are consistent with our theoretical predictions Lau02 (9) if we assume electron mobilities $``$0.4 m2 V-1s-1, $``$0.3 m2 V-1s-1 and $``$0.43 m2 V-1s-1 respectively. Although we do not have direct mobility measurements, these are reasonable room temperature values for such samples where optical phonon scattering and interface roughness are likely to be dominant Br98 (14). Figure 3 therefore gives a strong indication that there is a significant variation of spin relaxation rate with electric field at 300 K as predicted theoretically Lau02 (9).

Measurements on the mesa devices at different bias voltages at 170 K (Fig. 4a) support this conclusion (measurements at 300 K for the same range of voltages were prevented by excessive current in reverse bias leading to destruction of the device). The $`\mathrm{\Delta }\theta `$ signals show single-exponential decay and there is a strong variation of decay rate with voltage. The $`\mathrm{\Delta }R`$ signals also decay exponentially and Fig. 4b shows the extracted decay times ($`\tau _R`$) as a function of bias. For large negative bias $`\tau _R`$ increases significantly reflecting the reduced electron-hole overlap at high electric fields and showing that recombination is predominantly radiative Pol85 (18).

Figure 4c shows the spin relaxation rate vs bias voltage and vs the corresponding electric field at the QWs obtained using the layer thicknesses in the $`pin`$ structure and the band gap of GaAs at 170 K. Up to 20 kV cm-1 the spin relaxation rate is almost constant but then increases by about a factor of 10 by 80 kV cm-1. The solid dots are results of our non-perturbative calculations Lau02 (9, 12, 13) based on the DPK mechanism for a symmetrical (110)-oriented QW. The calculations assumed electron mobility of 0.6 m2 V-1s-1, which is consistent with the assumed 300 K value in the $`pin`$ wafer, $``$0.43 m2 V-1s-1, for mobility limited by optical phonon scattering in this temperature range. The variation for electric fields greater than 20 kV cm-1 (Fig. 4c) is very well fitted by the theory but at lower fields there is clearly a contribution to the spin relaxation not included in the calculations which predict a rate almost 100 times lower than we have observed at zero field.

The measured decay rates in this region were found to increase linearly with pump power, consistent with a small influence of the Bir, Aronov and Pikus (BAP) relaxation mechanism BAP76 (19) due to exchange interaction with the photo-excited holes as also found by Adachi et al. Ada01 (3). This does not, however, explain the observed discrepancy with theory since extrapolation to zero power would only make a reduction of about 10$`\%`$. Nor can it be explained in terms of the BAP mechanism with electrically injected holes; even for flat band conditions in forward bias, where the injected hole concentration would be greatest, the calculated concentration is at least two orders of magnitude too low to explain the additional spin relaxation BAP76 (19). Another possible cause of the additional spin relaxation is random built-in asymmetry in the QWs, which would give a field-independent SIA contribution. Such asymmetry could result from alloy fluctuations or from differences of top and bottom interface morphology for the GaAs/AlGaAs QWs and for a given $`𝐤`$ would generate an $`𝛀^{SIA}`$ component varying randomly in magnitude and orientation from point to point on a QW with associated non-zero mean square entering Eq. (1). To gauge this possibility, we have repeated our calculation in zero field for a QW with one perfect interface and the other containing a two-monolayer compositional gradient. The relaxation rate is $``$30 times greater than for a perfectly symmetrical QW (solid square in Fig. 4c). Even this small perturbation of the symmetry produces a change approaching our discrepancy between theory and experiment supporting the idea that the spin memory is actually limited by the sample perfection which, as indicated by the width of luminescence lines may be somewhat inferior to that in (001)-oriented samples. At sufficiently high reverse bias the increase of the relaxation rate due to the applied field becomes dominant over this random built-in asymmetry.

In conclusion, we have demonstrated that the enhanced high-temperature electron spin memory in (110)-oriented QWs may be varied by at least a factor 10 by application of modest gate bias voltages. The variation is consistent with the Rashba effect and the DPK spin-relaxation mechanism. On the basis of this mechanism, the measurements for low electric fields indicate that a further significant enhancement of spin memory may be achieved by modification of the growth techniques to optimise the interface morphology. This may be expected to give at 300 K spin memory longer than 10 ns Lau02 (9, 3) and straightforward voltage bias control simultaneously with high electron mobility. This combination of properties will facilitate a variety of spintronic devices.

###### Acknowledgements.

We wish to thank Professor E. L. Ivchenko, Mr M. M. Glazov and Dr. K. V. Kavokin for many illuminating discussions. This work was supported by EPSRC and DARPA/ARO.


text:
# 1 Introduction

## 1 Introduction

Ubiquitous of signals and processes with $`1/f`$ power spectral density at low frequencies has led to speculations that there might exist some generic mechanism underlying production of $`1/f`$ noise. The generic origins of two popular noises: white noise (no correlation in time, $`S(f)\mathrm{}1/f^0`$) and Brownian noise (no correlation between increments, $`S(f)\mathrm{}1/f^2`$) are very well known. It should be noted, that Brownian motion is the integral of white noise and that operation of integration of the signal increases the exponent by 2 while the inverse operation of differentiation decreases it by 2. Therefore, $`1/f`$ noise can not be obtained by simple procedure of integration or differentiation of the convenient signals. There are no simple, even stochastic, equations generating signals with $`1/f`$ noise. Note in this context also to the concept of the fractional Brownian motion and to the half-integration of a white noise signals used for generation of processes with $`1/f`$ noise . These and similar mathematical algorithms, procedures and models for generation of the processes with $`1/f`$ noise  are, however, sufficiently specific, formal or unphysical. They can not, as a rule, be solved analytically and they do not reveal the origin as well as the necessary and sufficient conditions for appearance of $`1/f`$ type fluctuations. Physical models of $`1/f`$ noise in some physical systems are usually very specialized, complicated and they do not explain the internal origin of the omnipresent processes with $`1/f^\delta `$ spectrum.

A lot of contributions are available in the literature concerning the origin of $`1/f`$ noise. On the web (http://linkage.rockefeller.edu/wli/1fnoise), Wentain Li is collective bibliography of flicker noise. Sufficiently comprehensive bibliography of the contributions concerning the modeling of $`1/f`$ noise may be find in references \[3–7\].

This work is a continuation of series of papers devoted to the modeling $`1/f`$ noise in simple systems  and search of necessary conditions for appearance of the signals with power spectrum at low frequencies like $`S(f)\mathrm{}1/f^\delta `$ ($`\delta 1`$) \[6–9\]. In papers  an analysis of the necessary conditions for appearance of $`1/f`$ type fluctuations in the simple systems consisting of few or even one particle and affected by random perturbations is presented. Later, a simple analytically solvable model of $`1/f`$ noise has been proposed , analyzed  and generalized . The model reveals main features, and parameter dependencies of the power spectrum of $`1/f`$ noise.

Here considering signals and currents as consisting of pulses generalizations and development of the model  are presented. The paper includes derivation of the expression for the correlation function, analysis of the examples of different signals and exhibition of the necessary conditions for appearance of $`1/f`$ type power spectrum in the signals consisting of pulses. It is shown that intrinsic origin of $`1/f`$ noise is a Brownian motion of the pulse interevent time.

## 2 The model

Let us consider currents or signals represented as sequences of random (but correlated) pulses $`A_k(tt_k)`$. Function $`A_k(tt_k)`$ represents the shape of the $`k`$-pulse of the signal in the region of the pulse occurrence time $`t_k`$. The signal or intensity of the current of particles in some space cross section may, therefore, be expressed as

$$I\left(t\right)=\underset{k}{}A_k\left(tt_k\right).$$
$`(1)`$
It is easy to show that the shapes of the pulses mainly influence the high frequency, $`f\mathrm{\Delta }t_p`$ with $`\mathrm{\Delta }t_p`$ being the characteristic pulse length, power spectrum while fluctuations of the pulse amplitudes result, as a rule, in the white or Lorentzian but not $`1/f`$ noise . Therefore, we restrict our analysis to the noise due to the correlations between the pulse occurrence times $`t_k`$. In such an approach we can replace the function $`A_k\left(tt_k\right)`$ by $`a\delta \left(tt_k\right)`$. Here $`\delta \left(tt_k\right)`$ is the Dirac delta function, $`a=_{\mathrm{}}^+\mathrm{}A_k\left(tt_k\right)dt`$ is the average area of the pulse and the brackets $`\mathrm{}`$ denote the averaging over realizations of the process. In such an approach the signal (1) my be expressed as

$$I\left(t\right)=a\underset{k}{}\delta \left(tt_k\right).$$
$`(2)`$
This model also corresponds to the flow of identical point objects: electrons, photons, cars and so on. On the other hand, fluctuations of the amplitudes $`A_k`$ may result in the additional noise but can not reduce $`1/f`$ noise we are looking for.

### 2.1 Power spectrum

Power spectral density of the signal $`I\left(t\right)`$ is

$$S\left(f\right)=\underset{T\mathrm{}}{lim}\frac{2}{T}\left|\underset{t_i}{\overset{t_f}{}}I\left(t\right)e^{i2\pi ft}dt\right|^2$$
$`(3)`$
where $`T=t_ft_i`$ is the observation time.

We can also introduce the autocorrelation function

$$\mathrm{\Phi }\left(s\right)=\frac{1}{T}\underset{t_i}{\overset{t_fs}{}}I\left(t\right)I\left(t+s\right)dt$$
$`(4)`$
and use the Wiener-Khintchine relations:

$$S\left(f\right)=4\underset{T\mathrm{}}{lim}\underset{0}{\overset{T}{}}\mathrm{\Phi }\left(s\right)\mathrm{cos}\left(2\pi fs\right)ds,$$
$`(5)`$
$$\mathrm{\Phi }\left(s\right)=\underset{0}{\overset{\mathrm{}}{}}S\left(f\right)\mathrm{cos}\left(2\pi fs\right)df.$$

Substitution of Eq. (2) into Eq. (3) results the power spectral density of the signal expressed as sequence of pulses

$$S\left(f\right)=\underset{T\mathrm{}}{lim}\frac{2a^2}{T}\left|\underset{k=k_{\mathrm{min}}}{\overset{k_{\mathrm{max}}}{}}e^{i2\pi ft_k}\right|^2$$

$$=\underset{T\mathrm{}}{lim}\frac{2a^2}{T}\underset{k=k_{\mathrm{min}}}{\overset{k_{\mathrm{max}}}{}}\underset{q=k_{\mathrm{min}}k}{\overset{k_{\mathrm{max}}k}{}}e^{i2\pi f\mathrm{\Delta }(k;q)}$$
$`(6)`$
where $`\mathrm{\Delta }(k;q)t_{k+q}t_k`$ is the difference of pulse occurrence times $`t_{k+q}`$ and $`t_k`$ while $`k_{\mathrm{min}}`$ and $`k_{\mathrm{max}}`$ are minimal and maximal values of index $`k`$ in the interval of observation $`T=t_ft_i`$.

For a stationary process Eq. (6) yields

$$S\left(f\right)=\frac{2a^2}{T}\underset{q=N}{\overset{N}{}}\left(N+1\left|q\right|\right)e^{i2\pi f\mathrm{\Delta }\left(q\right)}.$$

Here $`N=k_{\mathrm{max}}k_{\mathrm{min}}`$, the brackets $`\mathrm{}`$ denote the averaging over realizations of the process and over the time (index $`k`$) and a definition $`\mathrm{\Delta }\left(q\right)=\mathrm{\Delta }\left(q\right)=\mathrm{\Delta }(k;q)`$ is introduced. For abbreviation of the equations we have omitted the mark of the limit ”$`\underset{T\mathrm{}}{lim}`$” here and further on in expressions for the power spectrum $`S\left(f\right)`$.

The average over the distribution of $`\mathrm{\Delta }\left(q\right)`$ may be expressed as

$$e^{i2\pi f\mathrm{\Delta }\left(q\right)}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}e^{i2\pi f\mathrm{\Delta }\left(q\right)}\mathrm{\Psi }_\mathrm{\Delta }\left(\mathrm{\Delta }\left(q\right)\right)d\mathrm{\Delta }\left(q\right)$$

$$=\chi _{\mathrm{\Delta }\left(q\right)}\left(2\pi f\right).$$

Here $`\mathrm{\Psi }_\mathrm{\Delta }\left(\mathrm{\Delta }\left(q\right)\right)`$ is the distribution density of $`\mathrm{\Delta }\left(q\right)`$, and $`\chi _{\mathrm{\Delta }\left(q\right)}\left(2\pi f\right)`$ is the characteristic function of the distribution $`\mathrm{\Psi }_\mathrm{\Delta }\left(\mathrm{\Delta }\left(q\right)\right)`$. Therefore

$$S\left(f\right)=2a^2\underset{q=N}{\overset{N}{}}\left(\overline{\nu }\frac{\left|q\right|}{T}\right)\chi _{\mathrm{\Delta }\left(q\right)}\left(2\pi f\right)$$
$`(7)`$
where $`\overline{\nu }=`$$`\underset{T\mathrm{}}{lim}\frac{N+1}{T}`$ is the mean number of pulses per unit time.

When the sum $`\underset{q=N}{\overset{N}{}}\left|q\right|\chi _{\mathrm{\Delta }\left(q\right)}\left(2\pi f\right)`$ converges and $`T\mathrm{}`$ we have the power spectrum from Eq. (7) in the form

$$S\left(f\right)=2a\overline{I}\underset{q=N}{\overset{N}{}}\chi _{\mathrm{\Delta }\left(q\right)}\left(2\pi f\right).$$
$`(8)`$
Here $`\overline{I}\overline{I\left(t\right)}=\overline{\nu }a`$ is the average current.

### 2.2 Correlation function

Substitution of Eq. (2) into Eq. (4) yields the correlation function of the signal (2)

$$\mathrm{\Phi }\left(s\right)=\frac{a^2}{T}\underset{k,q}{}\delta \left(t_{k+q}t_ks\right).$$

After summation over index $`k`$ we have

$$\mathrm{\Phi }\left(s\right)=\overline{I}a\underset{q}{}\delta \left(\mathrm{\Delta }\left(q\right)s\right)$$

where the brackets $`\mathrm{}`$ denote again the averaging over realizations of the process and over the time (index $`k`$), as well. Such averaging coincides with the averaging over the distribution of the time difference $`\mathrm{\Delta }\left(q\right)`$

$$\mathrm{\Phi }\left(s\right)=\overline{I}a\underset{q}{}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\psi _\mathrm{\Delta }(\mathrm{\Delta })\delta \left(\mathrm{\Delta }s\right)d\mathrm{\Delta }$$

$$\mathrm{\Phi }\left(s\right)=\overline{I}a\delta \left(s\right)+\overline{I}a\underset{q0}{}\psi _\mathrm{\Delta }(s).$$
$`(9)`$
Here $`\psi _\mathrm{\Delta }(\mathrm{\Delta })`$ is the distribution density of $`\mathrm{\Delta }\left(q\right)`$. Substitution of Eq. (9) into Eq. (5) yields expressions (7) and (8).

## 3 Examples

Consider some examples of the signals represented by Eq. (2).

(i) Periodic signal expressed as $`I\left(t\right)=a_k\delta \left(tk\tau \right)`$ generates the power spectrum

$$\begin{array}{c}S\left(f\right)=2a^2\underset{T\mathrm{}}{lim}\frac{\mathrm{sin}^2\left(\pi \left(N+1\right)\tau f\right)}{T\mathrm{sin}^2\left(\pi \tau f\right)}2\overline{I}^2\delta \left(f\right),f\tau ^1.\end{array}$$

(ii) Perturbed periodic signal represented by Eq. (2) with the times series expressed as recurrence equations $`t_kt_{k1}\tau _k=\overline{\tau }+\sigma \epsilon _k`$ with $`\left\{\epsilon _k\right\}`$ being a sequence of uncorrelated normally distributed random variables with zero expectation and unit variance and $`\sigma `$ being the standard deviation of this white noise . For this model we have $`\mathrm{\Delta }\left(q\right)=q\overline{\tau }+\sigma \underset{l=k+1}{\overset{k+q}{}}\epsilon _l`$ and the characteristic function

$$\chi _{\mathrm{\Delta }\left(q\right)}\left(2\pi f\right)=e^{\mathrm{i2}\pi f\mathrm{\Delta }\left(q\right)\frac{1}{2}\left(2\pi f\right)^2\sigma _\mathrm{\Delta }^2}$$
$`(10)`$
where $`\mathrm{\Delta }\left(q\right)=q\overline{\tau }`$ and the variance $`\sigma _\mathrm{\Delta }^2`$ of the time difference $`\mathrm{\Delta }(q)`$ equals $`\sigma _\mathrm{\Delta }^2\mathrm{\Delta }\left(q\right)^2\mathrm{\Delta }\left(q\right)^2=\sigma ^2\left|q\right|`$. Substitution of Eq. (10) into Eq. (8) yields the Lorentzian spectrum

$$S\left(f\right)=\overline{I}^2\frac{4\tau _{\mathrm{rel}}}{1+\tau _{\mathrm{rel}}^2\omega ^2}$$
$`(11)`$
where $`\omega =2\pi f`$ and $`\tau _{\mathrm{rel}}=\sigma ^2/2\overline{\tau }`$.

(iii) Time difference $`\mathrm{\Delta }\left(q\right)=\underset{l=k+1}{\overset{k+q}{}}\tau _l`$ as a sum of uncorrelated interevent times $`\tau _l`$. According to Eqs. (6)–(8) we have in this case

$$S\left(f\right)=2a\overline{I}\left[1+2Re\underset{q=1}{\overset{N}{}}e^{i2\pi f\tau }^q\right]$$

$$S\left(f\right)=2a\overline{I}\left[1+2Re\frac{\chi _\tau \left(\omega \right)}{1\chi _\tau \left(\omega \right)}\right].$$
$`(12)`$
For instance, substitution at $`f\overline{\tau }^1`$ and $`f\sigma ^1`$ of Eq. (10) with $`q=1`$ into Eq. (12) results in Eq. (11).

(iv) For the Poisson process

$$\chi _\tau \left(2\pi f\right)=\frac{1}{1\mathrm{i2}\pi f\overline{\tau }},Re\frac{\chi _\tau \left(\omega \right)}{1\chi _\tau \left(\omega \right)}=0$$

and we have from Eq. (12) only the shot noise

$$S\left(f\right)=2a\overline{I}=S_{\mathrm{shot}}.$$
$`(13)`$
(v) Brownian motion of the interevent time $`\tau _k`$ with some restrictions, e.g., with the relaxation to the average value $`\overline{\tau }`$,

$$\tau _k=\tau _{k1}\gamma \left(\tau _{k1}\overline{\tau }\right)+\sigma \epsilon _k,$$
$`(14)`$
when the pulse occurrence times $`t_k`$ are expressed as

$$t_k=t_{k1}+\tau _k.$$
$`(15)`$
According to Eq. (6) the power spectrum of the signal (2) with the pulse occurrence times $`t_k`$ generated by Eqs. (14) and (15) for sufficiently small parameters $`\sigma `$ and $`\gamma `$ in any desirably wide range of frequencies, $`f_1=\gamma /\pi \sigma _\tau <f<f_2=1/\pi \sigma _\tau `$, is 1/f-like , i.e.,

$$S\left(f\right)=\overline{I}^2\sqrt{\frac{2}{\pi }}\frac{\overline{\tau }\mathrm{exp}\left(\overline{\tau }^2/2\sigma _\tau ^2\right)}{\sigma _\tau f}.$$
$`(16)`$
Here $`\sigma _\tau ^2=\sigma ^2/2\gamma `$ is the variance of the interevent time $`\tau _k`$.

## 4 Origin of $`1/f`$ noise

The origin for appearance of $`1/f`$ fluctuations in the model described in Eqs. (14) and (15) is related with the relatively slow, Brownian, fluctuations of the pulse interevent time. For this reason, the variance $`\sigma _\mathrm{\Delta }^2`$ of the time difference $`\mathrm{\Delta }(k;q)`$ for $`\left|q\right|\gamma ^1`$ is a quadratic function of the time difference and, consequently, of the difference $`q`$ of the pulse serial numbers $`k`$ , i.e.,

$$\sigma _\mathrm{\Delta }^2(k;q)=\sigma _\tau ^2\left(k\right)q^2.$$
$`(17)`$
Substitution of Eqs. (10) and (17) into Eq. (8) yields $`1/f`$ spectrum (16).

### 4.1 Generalization

For slowly fluctuating interevent time, the time difference $`\mathrm{\Delta }(k;q)`$ may be expressed as

$$\mathrm{\Delta }\left(q\right)=\underset{l=k+1}{\overset{k+q}{}}\tau _lq\tau $$
$`(18)`$
where $`\tau =\left(t_{k+q}t_k\right)/q`$ is the average interevent time in the time interval $`(t_k,t_{k+q})`$, a slowly fluctuating function of the arguments $`k`$ and $`q`$. In such an approach, the power spectrum according to Eq. (6) is

$$S\left(f\right)=2\overline{I}a\underset{q}{}e^{i2\pi fq\tau }$$
$`(19)`$
where

$$e^{i2\pi fq\tau }\underset{\mathrm{}}{\overset{+\mathrm{}}{}}e^{i2\pi fq\tau }\psi _\tau \left(\tau \right)d\tau =\chi _\tau \left(2\pi fq\right)$$

is the characteristic function of the distribution density $`\psi _\tau (\tau )`$ of the interevent time $`\tau `$. Therefore, the power spectrum according to Eq. (19) may be expressed as

$$S\left(f\right)2\overline{I}^2\overline{\tau }\mathrm{\Psi }_\tau \left(0\right)/f.$$
$`(20)`$
Here the property $`_{\mathrm{}}^+\mathrm{}\chi _\tau \left(x\right)dx=2\pi \mathrm{\Psi }_\tau \left(0\right)`$ of the characteristic function has been used.

### 4.2 Correlation function of $`1/f`$ noise

The correlation function of $`1/f`$ noise in the approximation (18) may be calculated according to Eq. (9), i.e.,

$$\begin{array}{c}\mathrm{\Phi }\left(s\right)=\overline{I}a\underset{q}{}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\psi _\tau (\tau )\delta \left(q\tau s\right)d\tau \\ =\overline{I}a\delta \left(s\right)+\overline{I}a\underset{q0}{}\psi _\tau (\frac{s}{q})\frac{1}{\left|q\right|}.\end{array}$$
$`(21)`$
For the Gaussian distribution of the interevent time $`\tau `$

$$\psi _\tau \left(\tau \right)=\frac{1}{\sqrt{2\pi }\sigma _\tau }\mathrm{exp}\left(\frac{\left(\tau \overline{\tau }\right)^2}{2\sigma _\tau ^2}\right)$$

the correlation function (21) reads as

$$\mathrm{\Phi }\left(s\right)=\frac{\overline{I}a}{\sqrt{2\pi }\sigma _\tau }\underset{q}{}\mathrm{e}^{\frac{\left(sq\overline{\tau }\right)^2}{2\sigma _\tau ^2q^2}}\frac{1}{\left|q\right|}.$$
$`(22)`$
It should be noted that the deviation of the variance $`\sigma _\mathrm{\Delta }^2`$ for large $`q`$ from the quadratic dependence (17) and approach to the linear function $`\sigma _\mathrm{\Delta }^2=2D_{t_k}\left|q\right|`$ ensures the convergence of sums (21) and (22) and, consequently, results in the Lorentzian power spectrum (11) at $`f0`$ . Here $`D_{t_k}`$ is the ”diffusion” coefficient of the pulse occurrence time $`t_k`$, related with the variance $`\sigma _{t_k}^2`$ of the pulse occurrence time as $`\sigma _{t_k}^2=2D_{t_k}k`$. For the model (14)–(15) $`D_{t_k}=\sigma ^2/2\gamma ^2`$.

The power spectra calculated according to Eq. (5) with the correlation functions (21) and (22) are expressed as Eq. (20) and Eq. (16), respectively.

## 5 Conclusions

From the above analysis we can conclude that the intrinsic origin of $`1/f`$ noise is a Brownian fluctuations of the interevent time of the signal pulses, similarly to the Brownian fluctuations of the signal amplitude resulting in $`1/f^2`$ noise. The random walk of the interevent time in the time axis is a property of the randomly perturbed or complex systems with the elements of selforganization.

## Acknowledgments

The author acknowledges useful discussions with Prof. A. Bastys, Prof. R. Katilius, Prof. A. Matulionis and Dr. T. Meškauskas.