text:
# Kibble-Zurek Mechanism in the Ginzburg Regime: Numerical Experiment in the Ising Model

## Abstract

Kibble-Zurek mechanism is a theory of defect formation in a non-equilibrium continuous phase transition. So far the theory has been successfully tested by numerical simulations and condensed matter experiments in a number of systems with small thermal fluctuations. This paper reports first numerical test of the mechanism in a system with large thermal fluctuations and strongly non-mean-field behavior: the two dimensional Ising model. The theory predicts correctly the initial density of defects that survive a quench from the disordered phase. However, before the system leaves the Ginzburg regime of large fluctuations most of these defects are annihilated and the final density is determined by the dynamics of the annihilation process only.

PACS 11.27.+d, 05.70.Fh, 98.80.Cq

Introduction.— In a system with a continuous phase transition an adiabatic change of a parameter of the system, like e.g. temperature, pressure or a coupling constant in a Hamiltonian, can drive the system from a disordered phase to an ordered one. A classic example is the paramagnet-ferromagnet transition in the two dimensional (2D) Ising model. Thermodynamics of continuous phase transitions has been intensively explored over many years. Two mayor achievements: the solution of the Ising model by Onsager and the renormalization group of Wilson were rewarded with a Nobel Prize in physics. The RG formalism revealed deep connections between statistical mechanics and quantum field theory.

A candidate theory of non-equilibrium phase transitions is the Kibble-Zurek mechanism (KZM) . Kibble pointed out  that in a non-equilibrium transition there is no time to fully develop the long range order characteristic for the ordered phase. As a result, the final state of the system is a mosaic of finite size ordered domains with different orientations of the order parameter in every domain. In a topologically non-trivial case this disorder takes the form of a finite density of topological defects. This qualitative idea was quantified more by Zurek . Zurek mechanism is a combination of two basic facts: (1) a divergent correlation length

$$\xi \xi _0|ϵ|^\nu ,$$
(1)
where $`ϵ`$ is a dimensionless distance from the critical point, $`\nu `$ is a critical exponent, and $`\xi _0`$ is a microscopic length scale, and (2) the critical slowing down or divergent relaxation time

$$\tau \tau _0|ϵ|^y.$$
(2)
Here $`\tau _0`$ is a microscopic time scale. Because of the divergent relaxation time any finite rate transition is a non-adiabatic phase transition: sufficiently close to the critical point (where $`ϵ=0`$) the system is too slow to react to the changing external parameter $`ϵ(t)`$. Close to $`ϵ=0`$ we can linearize

$$ϵ(t)=\frac{t}{\tau _Q}.$$
(3)
The relaxation time (2) equals the transition rate $`|ϵ/\frac{dϵ}{dt}|`$ at $`ϵ_Z(\tau _0/\tau _Q)^{\frac{1}{y+1}}`$ when the correlation length (1) is

$$\xi _Z\xi _0\left(\frac{\tau _Q}{\tau _0}\right)^{\frac{\nu }{y+1}}.$$
(4)
This Zurek length is the average size of the ordered domains after the phase transition and it determines the initial density of topological defects frozen into the ordered phase after a non-adiabatic continuous phase transition.

The original motivation for Kibble and Zurek were symmetry breaking phase transitions in cosmology. The random topological defects arising in such transitions might provide initial seeds for structure formation in the early Universe . However, the universality of phase transitions makes these ideas also relevant for a wide variety of condensed matter systems where they can be verified by experiment.

The KZM prediction (4) is supported by a number of numerical simulations . However, as a result of finite numerical resources these numerical data are limited to fast quenches (small $`\tau _Q`$) with a large $`ϵ_Z`$ in the regime of small fluctuations where one can use the mean field (MF) value of the critical exponent $`\nu _{\mathrm{MF}}=\frac{1}{2}`$. KZM is also supported by experiments in systems with small fluctuations like superfluid helium 3 , low $`T_c`$ superconductors , and fast quenches in high $`T_c`$ superconductors . In contrast, experiments in systems with large fluctuations like helium 4  or slow quenches in high $`T_c`$ superconductors  are inconclusive. Rivers suggested  that vortices created in the helium 4 experiment  disappear in a faster than expected annihilation. Due to technical difficulties the analytic calculations in Ref. eventually resort to a linearization equivalent to the mean-field theory. It is suggested there that beyond this linearized theory close to the critical point the annihilation rate is divergent. However, simulations in Ref. show that this effect may be not as dramatic as anticipated in Ref.. These authors suggest that because of the critical slowing down the annihilation rate close to $`ϵ=0`$ may in fact vanish. Due to limited numerical resources the numerical evidence in Ref. is rather indirect. To summarize, the problem of KZM in the Ginzburg regime of large fluctuations has been recognized  but is far from being settled.

At the moment we do not have any condensed matter or numerical experiment supporting KZM for large fluctuations and at the same time this is the regime where KZM in principle can be questioned on general grounds. The argument leading to Eq.(4) implicitly assumes that close to the critical point the divergent correlation length $`\xi `$ in Eq.(1) is the only relevant length scale. However, as is well known  but not quite generally appreciated, if $`\xi `$ were the only length scale, then, on dimensional grounds, all critical exponents would take their mean field values. As they do not (for example, in the 2D Ising model $`\nu =1`$ instead of the mean field $`\nu _{\mathrm{MF}}=\frac{1}{2}`$), then both $`\xi `$ and the microscopic $`\xi _0`$ must be relevant. With two relevant length scales the dimensional argument alone is not sufficient to determine the initial density of defects.

In this paper I report first numerical test of KZM for large fluctuations. As the critical regime is numerically demanding (large $`\xi `$ means large lattice and large $`\tau `$ means long time) I chose the simplest possible model - the celebrated 2D Ising model. This simple model has $`\nu =1`$ clearly different from the mean field $`\nu _{\mathrm{MF}}=1/2`$, and it has no regime where the MF theory might be at least remotely accurate. It is a perfect testing ground for KZM.

Ising model with Glauber dynamics.— Hamiltonian of the ferromagnetic Ising model in 2D is

$$H=\underset{i,j}{}S_iS_j.$$
(5)
Spins $`S_i\{1,+1\}`$ sit on a 2D $`N\times N`$ lattice with periodic boundary conditions, $`i,j`$ means a pair of nearest neighbor sites. The microscopic lengthscale $`\xi _0=1`$ is the lattice spacing. In all the following numerical simulations a $`1024\times 1024`$ lattice was used.

To study non-equilibrium phase transitions the Ising model has to be supplemented with dynamics. The standard choice is Glauber dynamics also known as Monte-Carlo with a heat bath . In the Glauber algorithm every time step consists of the following sub-steps:

* choose a random spin $`S_i`$ from the lattice,
* calculate its local field $`h_i=_{j\mathrm{n}.\mathrm{n}.i}S_j`$,
* calculate a probability $`P=\mathrm{exp}(\beta h_i)`$,
* choose a random number $`r[0,1]`$,
* if $`r>P`$ then set $`S_i=+1`$, else set $`S_i=1`$.

Here $`\beta `$ is an inverse temperature of the heat bath. This algorithm relaxes the state of the Ising model towards thermal equilibrium at a given $`\beta `$ . On average it takes $`N^2`$ steps to update the state of $`N^2`$ spins on the lattice. These $`N^2`$ steps define the microscopic time scale $`\tau _0`$ which I set equal to $`1`$.

The Ising model with Glauber dynamics belongs to the same universality class as the $`\varphi ^4`$ model with noise $`\eta `$

$$\tau _0\frac{}{t}\varphi =\xi _0^2^2\varphi \lambda (\varphi ^21)\varphi +\eta $$
(6)
so often employed in the numerical simulations of KZM . Here the continuum real field $`\varphi `$ is a coarse grained lattice spin $`S_i`$. The Ising model is an efficient “molecular dynamics” version of the $`\varphi ^4`$ field theory (6).

Relaxation time.— In order to estimate the exponent $`y`$ in Eq.(2) the relaxation time $`\tau `$ was measured for several values of $`\beta <\beta _c`$. For each $`\beta `$ the Ising model was initially prepared in a fully polarized state with all $`S_i=1`$, and then its average magnetization $`M=_iS_i/N^2`$ relaxed towards the equilibrium at $`M=0`$, see the insert in Fig.1. Each magnetization decay was fitted with an exponent $`M=\mathrm{exp}(t/\tau )`$. The best fits of $`\tau `$ are shown in the double logarithmic Fig.1 as a function of $`\beta _c\beta `$. The slope of the linear fit in Fig.1 gives an estimate of $`y=2.09\pm 0.02`$.

Quenches.— Phase transitions were simulated with a linear ramp of the inverse temperature

$$\beta (t)=1.5\frac{t}{\tau _Q}.$$
(7)
The initial state at $`t=0`$ was a state with random mutually uncorrelated spins - the state of equilibrium at $`\beta =0`$. Fig.2 shows density of domain walls separating positive $`S_i`$ from negative $`S_i`$ as a function of $`\beta `$ for a number of different transition times $`\tau _Q`$. The critical point is $`\beta _c=0.4407`$. For large $`\tau _Q`$ the density plots approach the equilibrium density $`n_{eq}(\beta )`$. Note that the equilibrium density $`n_{eq}(\beta )`$ of domain walls remains nonzero even for $`\beta >\beta _c`$. This is the critical Ginzburg regime of large fluctuations. A non-equilibrium transition with a finite $`\tau _Q`$ results in an additional non-equilibrium density $`dn(\beta )=n(\beta )n_{eq}(\beta )>0`$. KZM predicts that

$$dn_{\mathrm{KZM}}(\beta )\xi _Z^1=\tau _Q^{\frac{\nu }{y+1}}\tau _Q^{0.324\pm 0.003}.$$
(8)
Before this prediction is compared with the numerical data in Fig.2, let me digress on annihilation of domain walls.

Defect annihilation.— First example is annihilation of defects from an initially totally random spin configuration. The initial $`dn(t=0)`$ decays in time. Fig.3 shows the equilibrating $`n(t)`$ for several values of $`\beta >\beta _c`$. Each decay is fitted with a solid line that follows the power law $`dn(t)=(\tau _a/t)^{1/2}`$ with an exponent of $`1/2`$ known from the theory of phase ordering kinetics . The best fits are $`\tau _a=0.86\pm 0.05,0.93\pm 0.05,0.64\pm 0.05`$ for $`\beta =0.47,0.60,1.0`$ respectively. They are more or less constant in the considered range of temperatures: as the critical point is approached the time scale for annihilation $`\tau _a`$ neither diverges (as suggested in Ref.) nor vanishes (as suggested in Ref.), but remains finite and close to the microscopic $`\tau _0=1`$,

$$\tau _a\tau _0.$$
(9)
The quench time $`\tau _Q`$ determines the time available for defect annihilation. At late times after the transition, when most of the initial KZM domain walls are already annihilated, we expect the scaling

$$dn(\beta )\left(\frac{\tau _0}{\tau _Q}\right)^{\frac{1}{2}}.$$
(10)
It also follows from a phenomenological equation: $`\tau _0\frac{d}{dt}dn(t)=\frac{1}{2}dn^3(t)`$. Its solution is

$$dn(t)=\frac{dn(0)}{\sqrt{1+\frac{t}{\tau _0}dn^2(0)}}.$$
(11)
Note that at late times $`dn(t)`$ is forgetting its initial value $`dn(0)=dn_{\mathrm{KZM}}`$. This solution is an ilustration of the exact result (10) from Ref..

Second example is annihilation of domain walls from an initial state of equilibrium at $`\beta >\beta _c`$. The initial state was prepared by starting from fully polarized spins, all $`S_i=1`$, and then heating them up at $`\beta =0.45`$ for a time of $`10^5`$ sufficient to reach thermal equilibrium with $`n_{eq}(0.45)=0.20`$. Then at $`t=0`$ $`\beta `$ was suddenly increased (the heat bath was cooled) to $`\beta =0.55`$. Fig.3 shows $`n(t)`$ decaying towards the new equilibrium at $`n_{eq}(0.55)=0.075`$. This decay is much faster than for random initial spins because the equilibrium domain walls in the Ginzburg regime at $`\beta >\beta _c`$ are just boundaries of bubbles of the minority spin-down phase in the spin-up polarized ferromagnet. The bubbles together with their walls decay soon after the temperature is turned down.

KZM versus annihilation.— Figure 4 is a double logarithmic plot of the non-equilibrium density $`dn(\beta )`$ in Fig.2 as a function of $`\tau _Q`$ for a number of $`\beta `$s. The slope at the critical $`\beta _c=0.4407`$ is $`0.315\pm 0.007`$. This slope is consistent with the KZM slope (8) of $`0.324`$ and very different from a mean-field KZM slope of $`0.65`$ for $`\nu _{\mathrm{MF}}=1/2`$. The initial non-equilibrium density of domain walls is determined by KZM.

In contrast, similar slopes for $`\beta =1.0`$ and $`1.5`$ are $`0.45\pm 0.01`$ and $`0.48\pm 0.01`$ respectively, and they are consistent with the phase ordering kinetics exponent of $`1/2`$ in Eq.(10). Apparently at later times the system forgets the initial density $`dn_{\mathrm{KZM}}`$ and $`dn(\beta )`$ is determined solely by the dynamics of defect annihilation.

Indeed, the circles in Fig.4 show $`dn(\beta =1.5)`$ for a simulation where $`\beta (t)`$ is ramped up like in Eq.(7), but starting from the initial $`\beta _0=0.6>\beta _c`$. The spins were random at the initial $`\beta _0`$. The circles sit on the solid line which is a fit to $`dn(\beta =1.5)`$ obtained from a full quench like in Eq.(7). The annihilation dominated $`dn(\beta )`$ at later times is not sensitive to the details of the KZM of defect formation, compare Eqs.(10,11).

However, the defects that survive annihilation at later times are KZM defects quenched in from the disordered phase. As we have already seen, compare Fig.3, that annihilation of the Ginzburg domain walls is much faster than annihilation of defects from the initially random spin state. The latter state contains large domain walls while in the former domain walls are boundaries of bubbles of a minority spin phase. The points in Fig.(4) connected by a dashed line show $`dn(\beta =1.5)`$ after a quench starting from the equilibrium state at $`\beta _0=0.47`$ in the Ginzburg regime. These densities are orders of magnitude lower than densities from the full quenches starting at $`\beta =0`$: Ginzburg defects do not survive annihilation.

Conclusion.— I presented first numerical test of the Kibble-Zurek mechanism (KZM) in the Ginzburg regime of large thermal fluctuations. In this regime both the Zurek length $`\xi _Z`$ and the microscopic length $`\xi _0`$ are relevant length scales that determine the density of defects. However, the density of non-equilibrium defects frozen into the ordered phase by a quench from the disordered phase is determined by $`\xi _Z`$ only. This initial density of defects is gradually annihilated and when the system leaves the Ginzburg regime the density of defects is no longer sensitive to the details of the KZM, but it is determined by the dynamics of the annihilation process only. In particular, the dependence of the density on the transition rate is determined by an exponent that comes from the theory of phase ordering kinetics and not from the KZM. The only way to see the KZM scaling (8) directly is to measure the amount of disorder close to the critical point where the non-equilibrium KZM density is largely obscured by the prevailing equilibrium thermal fluctuations. However, the defects that survive the annihilation are the KZM defects quenched in from the high temperature phase, the defects quenched in from the Ginzburg regime decay much faster. The surviving defects are a clear, though indirect, signature of the KZM.

Acknowledgements.— I appreciate discussions with Henryk Arodź, Józef Sznajd, and Wojciech Żurek. This research was supported in part by the KBN grant 2 P03B 092 23 and in part by COSLAB.


format: text

text:
# Coulomb gap and variable range hopping in a pinned Wigner crystal

## Abstract

It is shown that pinning of the electron Wigner crystal by a small concentration of charged impurities creates the finite density of charged localized states near the Fermi level. In the case of residual impurities in the spacer this density of states is related to nonlinear screening of a close acceptor by a Wigner crystal vacancy. On the other hand, intentional doping by a remote layer of donors is a source of a long range potential, which generates dislocations in Wigner crystal. Dislocations in turn create charged localized states near the Fermi level. In both cases Coulomb interaction of localized states leads to the soft Coulomb gap and ES variable range hopping at low enough temperatures.

The growth of GaAs-heterostructures mobilities has made possible to study collective properties of the hole gas with a very small concentration. The average distance between holes $`r_s`$, calculated in units of the hole Bohr radius, has reached 80 and exceeded the theoretically predicted freezing point of the hole liquid, $`r_s=38`$, where it should become the Wigner crystal. The transport of such low density systems with and without magnetic field was studied recently. If a strong magnetic field applied along the hole plane the low temperature conductivity of the Wigner-crystal-like phase seems to be related to the variable range hopping (VRH) Noh . Thus, it seems timely to discuss ohmic transport of the Wigner crystal. Of course, any Wigner crystal is always pinned by impurities and can not slide as a whole in a small electric field. Interstitials and vacancies of the Wigner crystal provide conductivity but they have large activation energies and freeze out at low temperatures. On the other hand, the VRH conductivity requires nonzero density of states of charge excitations near the Fermi level, both in the case of Mott law or Efros-Shklovskii (ES) law. If pinning centers produce such ”bare” density of states then the long range Coulomb interaction creates the Coulomb gap around the Fermi level and leads to the ES variable range hopping. In this paper, I study the origin of the bare density of states of the pinned Wigner crystal, the Coulomb gap on its background and the resulting VRH conductivity at low temperatures.

Let me start from a clean electron Wigner crystal on uniform positive background. Consider an additional electron added to the perfect Wigner crystal and let the lattice relax locally. This addition costs the interstitial energy $`0.14e^2/\kappa a`$, where $`a`$ is the lattice constant or the optimal triangular lattice, $`\kappa `$ is the dielectric constant of GaAs Fisher . On the other hand, the energy cost of the extraction of an electron and the corresponding relaxation around the vacancy is equal to the vacancy energy $`0.26e^2/\kappa a`$. Thus, the density of states of relaxed excitations (electronic polarons ES84 ), consists of two delta-like peaks at energies $`0.26e^2/\kappa a`$ and $`0.14e^2/\kappa a`$. The lower peak corresponds to occupied states, the upper one contains empty states. Between them there is the hard (completely empty) gap of the width $`0.40e^2/\kappa a`$. The Fermi level of electronic polarons (the zero of energy) is situated in this gap, therefore the density of states at the Fermi level is zero.

In typical samples the Wigner crystal is pinned by charged impurities: donor and acceptors. Even if there is no intentional doping and the two-dimensional electrons gas (2DEG) is created in a heterostructure by a distant metallic gate, there is always a small concentration of residual donors and acceptors. Many high mobility heterostructures are intentionally doped by a donor layer, which is separated from 2DEG by the spacer with a large width $`sa`$. This is an additional source of disorder. In the beginning of this paper we neglect this source and deal with residual impurities only. We return to the role of doped layer in the second part of this paper.

Let me argue that close residual donors and acceptors create localized states at the Fermi level. I use theory of pinning of the Wigner crystal by charge impurities developed in Ref. Rouzin, . As in that paper I talk about Wigner crystal of electrons and concentrate on close residual acceptors which according to Ref. Rouzin,  produce stronger pinning than donors. (For the hole crystal the role of acceptors is of course played by donors.) Because the spacer width, $`s`$, and the lattice constant of the Wigner crystal, $`a`$, are much larger than the lattice constant of AlGaAs the authors of Ref. Rouzin,  assumed that acceptors are randomly distributed in space with three-dimensional concentration $`N`$. Each acceptor is negative because it has captured one electron from 2DEG.

It was found Rouzin  that effect of an acceptor dramatically depends on its distance, $`d`$, from the plane of the Wigner crystal. If $`d`$ is larger than $`0.68a`$ the electron crystal adjusts to acceptor in such a way that in the ground state an interstitial position of the crystal is exactly above the acceptor. In this case, interaction of the acceptor with the Wigner crystal can be calculated in the elastic approximation. On the other hand, when $`d`$ is smaller than $`0.68a`$ acceptor creates a vacancy in the crystal, which positions itself right above the acceptor. This means that the close acceptor effectively builds in the crystal cite pushing away to infinity its electron. We can consider the latter acceptor as an empty localized state for electron, while the former acceptor can be considered as a localized state occupied by an electron. Their energies are equal at $`d=0.68a`$, so that all acceptors with $`d>0.68a`$ are occupied and all acceptor with $`d<0.68a`$ are empty. Thus, acceptors with $`d=0.68a`$ are at the Fermi level and the bare density of states near the Fermi level

$$g_B\frac{Na}{e^2/\kappa a},$$
(1)
where $`Na`$ is an estimate of the two-dimensional concentration of acceptors located within distance $`a`$ (only they provide substantial pinning) and $`e^2/\kappa a`$ is the characteristic energy of their pinning. The nonzero density of states $`g_B`$ makes possible VRH conductivity of the pinned Wigner crystal. For example, an electron from an acceptor with $`d=0.68a+0`$ can hop to a distant one with $`d=0.68a0`$.

The long range Coulomb interaction between localized states leads to the soft Coulomb gap at the Fermi level ES . Shape of the Coulomb gap depends on the interaction, $`U(r)`$, of the hopping electron in the final state of the hop and the empty place it has left behind (excitonic term) ES ; ES84 . If the interaction can be described by the standard Coulomb law $`U_C(r)=e^2/\kappa r`$, where $`r`$ is the two-dimensional distance between initial and final state, then the Coulomb gap has the standard for two-dimensional case form

$$g(E)=\frac{2}{\pi }\frac{|E|\kappa ^2}{e^4}.$$
(2)
The attraction potential $`U(r)`$ generally speaking is different from $`U_C(r)`$ due to screening of the Coulomb interaction by elastic deformations of the Wigner crystal or in other words by its polarization. In order to calculate Fourier image $`U(q)`$ of $`U(r)`$ we have to introduce the Larkin length $`L`$ at which acceptors destroy the long range order of the crystal. In the two-dimensional Wigner crystal pinning by close acceptors is so strong that this length is close to the average distance between them, $`L(Na)^{1/2}`$. Using $`L`$, we can write

$$U(q)=2\pi /q\kappa ϵ(q),$$
(3)
where

$$ϵ(q)=1+\frac{2}{qr_D}\frac{q^2}{q^2L^2}.$$
(4)
is given by interpolation between large and small $`q`$ cases. Here $`r_D0.32a`$ is the linear Debye screening radius of the Wigner crystal. (The asymptotics of Eq. (4) at large and small $`q`$ are similar to those of the dielectric constant of a two-dimensional electron gas in a strong magnetic field Kukushkin ; Aleiner ; Koulakov , where the Larmour radius plays the role of $`L`$.)

In the real space the Fourier image of Eq. (4) results in the Coulomb potential $`U(r)=U_C(r)`$ only at $`rL^2/a=1/Na^2`$. One can say that by these distances all electric lines of a probe charge located in the reference point in the plane of Wigner crystal leave the plane. This leads to the standard Coulomb gap Eq. (2), but only at energies $`|E|<e^2Na^2/\kappa `$. At smaller distances $`r1/Na^2`$ potential $`U(r)`$ is weaker than $`U_C(r)`$ and grows only logarithmically with decreasing $`r`$. But this does not lead to faster (exponential) growth of the density of states at $`|E|e^2Na^2/\kappa `$ as it would if $`g_B`$ were very large. The reason is that at energy $`|E|e^2Na^2/\kappa `$ the Coulomb gap density of states, Eq. (2), already reaches the bare density of states, $`g_B`$, given by Eq. (1). The Coulomb gap is just a depletion of the density of states on the background of the bare density. Therefore, the density of states in the Coulomb gap can not be larger than $`g_B`$ ES84 . This means that at $`|E|e^2Na^2/\kappa `$ the density of states saturates at the level of Eq. (1).

The role of distance $`L^2/a=1/Na^2`$ can be interpreted in another way. When we add an electron to the clean Wigner crystal it makes an interstitial but its charge spreads to infinity. In the a pinned Wigner crystal an added charge is smeared in the disc of the finite radius $`R`$. We can find $`R`$ minimizing the sum of the Coulomb energy $`e^2/\kappa r`$ of the disc of radius $`r`$ and the shear energy $`\mu (u/L)^2r^2`$ necessary for the disc dilatation by the area $`a^2`$. Here $`\mu e^2/a^3`$ is the shear modulus and $`ua^2/r`$ is the necessary displacement. It is important that the characteristic distance of the variation of the shear displacement (which percolates between strongly pinning acceptors) is the average distance between close acceptors $`L(Na)^{1/2}`$. Minimizing this sum we find that $`R=L^2/a=1/Na^2`$. This again tells us that at distances larger than $`L^2/a=1/Na^2`$ we deal with the unscreened Coulomb potential.

At low enough temperatures the Coulomb gap leads to Efros-Shklovskii(ES) law of the temperature dependence of the variable range hopping ES ; ES84

$$\sigma =\sigma _0\mathrm{exp}[(T_{ES}/T)^{1/2}].$$
(5)
Here $`\sigma _0`$ is a prefactor, which has only an algebraic $`T`$-dependence and

$$T_{ES}=Ce^2/\kappa \xi ,$$
(6)
where $`C`$ is a numerical coefficient close to 6 and $`\xi `$ is the localization length for electron tunnelling with energy close to the Fermi level in the Wigner crystal. We will discuss the value of $`\xi `$ in the end of the paper.

When with increasing temperature the width of the band of energies used for ES hopping $`(TT_{ES})^{1/2}`$ reaches $`e^2Na^2/\kappa `$ ES law is replaced by the Mott conductivity

$$\sigma =\sigma _0\mathrm{exp}[(T_M/T)^{1/3}],$$
(7)
where $`T_M=C_M/(g_B\xi ^2)`$. Transition from the ES law to the Mott law happens at $`TT_{ES}^3/T_M^2T_{ES}(Na^2\xi )^2`$. Thus, while ES law does not depend on the acceptor concentration $`N`$ the range of ES law shrinks when concentration of acceptors decreases.

This example emphasizes universality of the Coulomb gap and ES law. The Coulomb gap was derived for lightly doped semiconductors, where disorder is as strong as interactions ES . Later Efros Efros  suggested a model, where electrons are located on sites of a square lattice with +1/2 charges, the number of sites being twice larger than number of electrons. Random energies of sites are uniformly distributed in the band $`A`$ (they are measured in units of Coulomb interaction of two electrons on nearest sites). In this model the Coulomb gap does not survive in the small disorder case $`A<<1`$ because in this case in spite of small disorder positive charges of empty sites and negative charges of occupied sites alternate in the perfect NaCl-like order. One could, therefore, say that the Coulomb gap is the property of strong disorder only. The above example of the Wigner crystal on the uniform background pinned by rare strong impurities however shows that such impression may be misleading. Even in the case when concentration of acceptors $`N`$ is small (one can call it the weak disorder case) the Coulomb gap and ES law survive.

One can interpret what happens with the Wigner crystal in the terms of the Efros model. At strong disorder when $`Na^31`$ our model is close to the Efros model with $`A1`$. When parameter $`Na^3`$ becomes much smaller than unity and the disorder becomes weak, the energy scatter of states created by impurities stays at the level of $`e^2/\kappa a`$, while the long Coulomb interaction between these states becomes much smaller, of the order of $`e^2(Na)^{1/2}kappa`$. Introducing a two-dimensional random lattice with cites occupied by acceptors we can come to a renormalized Efros model with $`A1/(Na^3)^{1/2}1`$. This brings us again to the Coulomb gap at small energies.

Until now we talked about the role of residual acceptors. Let us qualitatively discuss intentionally doped heterojunctions where donors are situated in the narrow layer parallel to the plane of the 2DEG (delta-doping) at a large distance $`sa`$ from 2DEG. Random distribution of donors in this layer creates fluctuations of their potential of all wavelengths, but only harmonics of the random potential with wavelengths larger than $`s`$ reach 2DEG. Let us assume that the two-dimensional concentration of charged donors in the doped layer is $`N_2`$. In principle $`N_2`$ can be different from the two-dimensional concentration of electrons in the Wigner crystal $`n`$ (2DEG can be compensated by acceptors or created by illumination). If $`N_2n`$ the Wigner crystal screens external potential by small, purely elastic deformations. Due to these deformations energies of an interstitial and a vacancy depend on a coordinates. Therefore, both delta-like peaks of the density of states are somewhat smeared. Their width, however, is much smaller than the hard gap between them. Thus, the hard gap is preserved and no states at the Fermi level appear. When the concentration of donors, $`N_2`$, grows, the amplitude $`u`$ of displacements of electrons in the Wigner crystal with wavelength of the order $`s`$ reaches the lattice constant of the crystal $`a`$. Such large deformations resolve themselves by creation of dislocations FisherF ; Fertig , because the energy price of the dislocation core becomes smaller than an elastic energy, which is eliminated by the dislocation. This happens when the donor concentration, $`N_2`$, is of the order of electron concentration, $`n`$. Indeed, if we cover the layer of donors by squares with the side $`s`$ the typical fluctuation of number of charges in a square is equal $`(N_2s^2)^{1/2}`$. Potential of the charge of this fluctuation reaches 2DEG practically without compensation by oppositely charged fluctuations in neighboring squares. As a result at $`n=N_2`$ 2DEG has to provide $`(Ns^2)^{1/2}=s/a`$ new electrons to screen random potential in a square. This may be done by an additional raw of electrons in the square or, in other words, by two dislocations. Thus approximately one dislocation appears in a square with a side $`s`$. This picture is actually a simple visualization of the Larkin domain. Below we concentrate on the case $`Nn`$. In this case, the length $`s`$ plays the role of the Larkin length $`L`$.

An isolated dislocation brings the electronic polaron state right to the Fermi level. Indeed, if we add an electron from the Fermi level to the end of an additional raw terminated by a dislocation and let the rest of electrons relax, the dislocation just moves along this additional raw by one lattice constant $`a`$ and the energy remains unchanged. Similarly, if we extract an electron from the end of the raw to the Fermi level, the dislocation moves in the opposite direction and the energy does not change. This means that if we neglect the interaction of dislocation with long range ptential of donors and interaction between dislocations they create a delta peak of the density of states right at the Fermi level. This peak is normalized on the concentration of dislocations. The Fermi level is pinned in the middle of this peak.

In the long range fluctuating potential of donors the peak of the density of dislocation states is smeared, because dislocations strongly interact with the gradient of potential. The peak is smeared also due to the interaction between the dislocations. When an additional electron is absorbed by a dislocation and, as a result, the dislocation moves by one lattice site, the logarithmic interaction with other dislocation changes.

To understand the role of Coulomb interaction in the density of states of the pinned Wigner crystal one has to concentrate on the fate of the charge of an added electron. In the clean Wigner crystal, if we have a single dislocation and move it by one lattice constant adding a new electron, the charge of this electron spreads to infinity, leaving the dislocation neutral.

If we have a gas of dislocations in the positions fixed by disorder and their interaction, a single charge can only spread to the finite distance, $`R`$. It can be estimated in the way we did for residual randomly distributed acceptors. If we assume that dislocation are fixed in space by fluctuating donor potential the extra electron charge spreads optimizing the sum of its Coulomb energy and the energy of the shear deformation. Optimization leads to $`Rs`$. Thus, the Coulomb potential of charges inside the pinned Wigner crystal is valid at distances in the plane $`rs`$.

This system is clearly similar to the Coulomb glass and therefore has the Coulomb gap. Indeed, as we mentioned above the derivation of the Coulomb gap is based upon the observation that when a localized electron is transferred to another localized state one should take into account its 1/r Coulomb attraction with the hole it has left (the excitonic effect). We claimed above that that the Coulomb interaction is valid if $`rs`$. Therefore, the linear in energy Coulomb gap appears at the Fermi level. The width of this gap is $`e^2/\kappa R=e^2/\kappa s`$. At $`sa`$ the Coulomb gap occupies only small fraction of the energy range between interstitial and vacancy peaks. Away from the Coulomb gap the density of states is almost constant. Thus, in the case of intentional $`\delta `$-doping by remote donors we again arrive to the Coulomb gap of density of states and correspondingly to ES law at low temperatures.

Let us discuss the value of localization length $`\xi `$ in Eq. (6). In classical Wigner crystal $`\xi a/r_s^{1/2}`$. This is a quite small value which leads to too large $`T_0`$ and very large resistances in the range of ES law. However, close to the melting point of the Wigner crystal $`\xi `$ can be much larger, making an observation of ES law in pinned Wigner crystal more realistic. Magnetoresistance, observed in Ref. Noh,  may be related to ES variable range hopping conductivity.

In conclusion, I emphasize again that I am dealing with the Wigner crystal which in the absence of impurities slides on the positive background. Impurities pin this Wigner crystal and lead to ES VRH conductivity due to hops between its pieces. This situation is similar to what happens in a system of many quantum dots situated in random electrostatic potential of stray charges or in a system of many electron puddles with random positive charge of their background. Similar physics was recently theoretically studied in quasi-one-dimensional systems Fogler .

I am grateful to M. M. Fogler, A. L. Efros and S. Teber for many important discussions. This paper is supported by NSF DMR-9985785.


text:
# Complementarity between Gauge-Boson Compositeness and Asymptotic Freedom — with scalar matter

## Abstract

We derive and solve the compositeness condition for the SU($`N_c`$) gauge boson coupling with $`N_s`$ scalar fields at the next-to-leading order in $`1/N_s`$ and the leading order in $`\mathrm{ln}\mathrm{\Lambda }^2`$ ($`\mathrm{\Lambda }`$ is the compositeness scale). It turns out that the argument of gauge-boson compositeness (with a large $`\mathrm{\Lambda }`$) is successful only when $`N_s/N_c>22`$, in which the asymptotic freedom fails, as is in the previously investigated case with fermionic matter.

The gauge theories have been widely applied not only in high energy physics, but also in various branches of physics, and the gauge bosons often appear as composites of matter fields cg geo . For example, some people considered the composite models of weak bosons and other gauge bosons cgw , some developed theory of induced gravity with composite graviton indgra , some interpreted the vector mesons (which is known to be composite) as gauge bosons with hidden local symmetry hidden , and others extensively studied dynamically induced ‘connections’ (i. e. gauge fields) in models of spacetime emb ; qm , and in molecular, nuclear, and solid-state systems geo . A gauge boson interacting with matter can be interpreted as composite under the compositeness condition $`Z_3=0`$ cc , where $`Z_3`$ is the wave-function renormalization constant of the gauge boson. The compositeness condition gives relations between the effective coupling constants and the compositeness scale $`\mathrm{\Lambda }`$. In our previous papers ccA cgb2 , we derived and solved the compositeness condition at the next-to-leading order in $`1/N`$ in various models, and investigated its implications, where $`N`$ is number of the matter field species. In particular in cgb2 , considering compositeness of the SU($`N_c`$) gauge boson coupling with $`N_f`$ fermion fields, we found a “complementarity” between the gauge-boson compositeness and asymptotic freedom af  of the gauge theory. Namely, the gauge-boson can be considered as a composite due to $`Z_3=0`$ only when $`N_f/N_c>11/2`$, in which the asymptotic freedom fails. Thus it is urgent for us to investigate if such a complementarity is accidental or really universal, by, for example, examining other cases of theories. In this paper, we perform a similar investigation for the SU($`N_c`$) gauge boson coupling with $`N_s`$ scalar fields instead of the fermion fields. We find that the complementarity again holds, i.e. that the gauge-boson is considered as a composite only when $`N_s/N_c>22`$, in which the asymptotic freedom fails.

We consider the SU($`N_c`$) gauge theory for the gauge boson $`G_\mu ^a(a=1,\mathrm{},N_c^21)`$ coupling with $`N_s`$ complex scalar fields $`\varphi _j(j=1,\mathrm{},N_s)`$, each of which belongs to the fundamental representation of SU($`N_c`$):

$``$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(G_{\mu \nu }^a\right)^2+{\displaystyle \underset{j}{}}\left[|𝒟_\mu \varphi _j|^2m^2|\varphi _j|^2\lambda (|\varphi _j|^2)^2\right]`$ (2)
$`{\displaystyle \frac{1}{2\alpha }}\left(^\mu G_\mu ^a\right)^2+^\mu \eta ^a\left(_\mu \eta ^agf^{abc}\eta ^bG_\mu ^c\right)`$
with $`G_{\mu \nu }^a=_\mu G_\nu ^a_\nu G_\mu ^a+gf^{abc}G_\mu ^bG_\nu ^c`$ and $`𝒟_\mu \varphi _j=_\mu \varphi _jigT^aG_\mu ^a\varphi _j`$, where $`g`$ and $`\lambda `$ are the coupling constants, $`m`$ is the mass of $`\varphi _j`$, $`f^{abc}`$ is the structure constant of SU($`N_c`$), $`T^a`$ is the representation matrix for the basis of the associated Lie algebra su($`N_c`$), $`\alpha `$ is the gauge fixing parameter, and $`\eta ^a(a=1,\mathrm{},N_c^21)`$ is the Fadeev-Popov ghost. In order to absorb the ultraviolet divergences arising from quantum fluctuations, we renormalize the fields, mass, and coupling constants as

$``$ $`=`$ $`{\displaystyle \frac{1}{4}}Z_3\left(G_{\mathrm{r}\mu \nu }^a\right)^2+{\displaystyle \underset{j}{}}\left[Z_2|𝒟_\mu \varphi _{\mathrm{r}j}|^2Z_mm_\mathrm{r}^2|\varphi _{\mathrm{r}j}|^2Z_\lambda \lambda _\mathrm{r}(|\varphi _{\mathrm{r}j}|^2)^2\right]`$ (4)
$`{\displaystyle \frac{Z_3}{2Z_\alpha \alpha _\mathrm{r}}}\left(^\mu G_{\mathrm{r}\mu }^a\right)^2+Z_\eta ^\mu \eta _\mathrm{r}^a\left(_\mu \eta _\mathrm{r}^a{\displaystyle \frac{Z_1}{Z_2}}g_\mathrm{r}f^{abc}\eta _\mathrm{r}^bG_{\mathrm{r}\mu }^c\right)`$
with $`G_{\mathrm{r}\mu \nu }^a=_\mu G_{\mathrm{r}\nu }^a_\nu G_{\mathrm{r}\mu }^a+(Z_1/Z_2)g_\mathrm{r}f^{abc}G_{\mathrm{r}\mu }^bG_{\mathrm{r}\nu }^c`$ and $`𝒟_\mu \varphi _{\mathrm{r}j}=_\mu \varphi _{\mathrm{r}j}i(Z_1/Z_2)g_\mathrm{r}T^aG_{\mathrm{r}\mu }^a\varphi _{\mathrm{r}j}`$, where the quantities with the index “r” are the renormalized ones, and $`Z_1`$, $`Z_2`$, $`Z_3`$, $`Z_\eta `$, $`Z_m`$, $`Z_\lambda `$, and $`Z_\alpha `$ are the renormalization constants. Since the degree of freedom of the $`Z_\alpha `$ can be absorbed by that of the non-physical parameter $`\alpha _\mathrm{r}`$, we simply fix it as $`Z_\alpha =Z_3`$.

If we impose the compositeness condition $`Z_3=0`$, the kinetic term of $`G_{\mathrm{r}\mu }^a`$ disappear, and becomes a non-dynamical auxiliary field, whose Euler equation reduces to a constraint. The constraint, however, causes a self-interaction of $`\varphi _\mathrm{r}`$, which gives rise to a composite gauge field. The kinetic term is reproduced through the quantum effects of the matter fields. This can be treated in the most rigorous way by following the usual renormalization procedure with (4), and by imposing the compositeness condition $`Z_3=0`$ after that. As was argued in the case of the model with fermionic matters cgb2 , the compositeness condition itself spoils the ordinary perturbation expansion in the coupling constant $`g_\mathrm{r}`$, and the appropriate expansion parameter is the inverse of the number of the matter-field species, $`1/N_s`$. Here we calculate $`Z_3`$ based on the Lagrangian (4) in the leading and the next-to-leading orders in $`1/N_s`$, solve the condition $`Z_3=0`$ for the coupling constant $`g_\mathrm{r}`$, and consider its implications with a particular attention on the complementarity between gauge-boson compositeness and asymptotic freedom.

In $`1/N_s`$ expansion, $`\lambda `$ should be $`O(1/N_s)`$, because otherwise diagrams with the more loops become the larger. We assign $`\lambda O(1/N_s)`$. We first choose the renormalization constant $`Z_m`$ so as to cancel out the constant (momentum-independent) contributions from the scalar self-mass parts order by order. Then the renormalization constant $`Z_3`$ should be chosen so as to cancel out all the divergences in the diagram A in Fig. 1 at the leading order in $`1/N_s`$ and in the diagrams B – N at the next-to-leading order. There is no mixing between the gauge boson line with one-scalar-loops inserted (Fig. 2a) and the scalar-loop chain due to $`\lambda _\mathrm{r}|\varphi _{\mathrm{r}j}|^4`$ interaction (Fig. 2b), because the sub-diagram in Fig. 3a vanishes. In calculation, we adopt the minimal subtraction prescription in the dimensional regularization. To avoid the absurdity of vanishing coupling constants, we should keep $`ϵ=(4d)/2`$ at a very small but non-vanishing value, where $`d`$ is the analytically continued number of the spacetime dimensions. We retain only the part leading in $`ϵ`$, which amounts to retaining the terms of $`O(ϵ^0)`$ because the compositeness condition implies $`g_\mathrm{r}^2=O(ϵ)`$ and each diagram has at most the order of magnitude of $`g^{2n}/ϵ^n=O(ϵ^0)`$.

Cancellations between the sub-diagrams in Fig. 3b extremely facilitate the calculations. Because the one-scalar-loop (denoted by $`\mathrm{\Pi }_0^{\mu \nu }`$) inserted into a gauge-boson propagator (with momentum $`q`$) (Fig. 2a) is divergenceless (i.e. $`q_\mu \mathrm{\Pi }_0^{\mu \nu }=0`$), it suppresses the $`\alpha _\mathrm{r}`$-dependent part of the diagram. After lengthy calculations similar to those indicated in Ref. cgb2 , we obtain the contributions from the A – H to $`Z_3`$ as are shown in Table 1, where we also cited the corresponding results for the model with fermionic matters in the previous paper cgb2 . The diagrams I and J in Fig. 1 contribute no leading divergences, because the scalar-loop chain behaves like $`(q^2)^{lϵ}`$, ($`q`$ is the total momentum transfer through the chain, and $`l`$ is the number of loops in the chain.) and consequently the integration over $`q`$ give rise to an extra factor of $`ϵ`$ in comparison with the leading contributions in $`ϵ`$. The diagram K is not leading in $`ϵ`$ because the sub-diagram Fig. 3c converges, and the diagrams L – N cancel out because the sub-diagrams in Fig. 3d cancel,

Next, we renormalize the sub-diagram divergences by subtracting the divergent counter parts of (i) each scalar loop inserted into the gauge-boson lines in D – H, (ii) the scalar self-energy part in D, (iii) the scalar-scalar-gauge-boson vertex parts in D, E, G, and H, and (iv) the three-gauge-boson vertex part in G and H. The contributions from the counter parts cancel out for the diagrams with odd loops, and amount to minus twice the original terms for those with even loops. Collecting all the contributions in Table 1 and the counter terms, we finally obtain the compact expression

$`Z_3`$ $`=`$ $`1{\displaystyle \frac{1}{6}}N_sg_\mathrm{r}^2I+{\displaystyle \frac{11}{3}}N_cg_\mathrm{r}^2I{\displaystyle \frac{\alpha _\mathrm{r}}{2}}N_cg_\mathrm{r}^2I(1{\displaystyle \frac{1}{6}}N_sg_\mathrm{r}^2I)`$ (6)
$`+{\displaystyle \frac{3}{2}}N_c({\displaystyle \frac{6}{N_s}}g_\mathrm{r}^2I)\mathrm{ln}(1{\displaystyle \frac{1}{6}}N_sg_\mathrm{r}^2I)+O({\displaystyle \frac{1}{N_s^3}}),`$
where $`I=1/16\pi ^2ϵ`$. It is remarkable that this is much similar in form to the corresponding result for the previous model with fermionic matters (cgb2 ), in spite that the contributions from the corresponding diagram differ. In fact (6) is what is obtained by replacing $`N_f`$ by $`N_s/4`$ in the result Eq. (12) in cgb2 . We note that $`Z_3`$ does not depend on the scalar self-coupling constant $`\lambda _\mathrm{r}`$.

The compositeness condition $`Z_3=0`$ with expression (6) can be solved for $`g_\mathrm{r}^2`$ by iterating the leading-order solution into itself. The solution is rather simple:

$`g_\mathrm{r}^2={\displaystyle \frac{6}{N_sI}}\left[1+{\displaystyle \frac{22N_c}{N_s}}+O({\displaystyle \frac{1}{N_s^2}})\right].`$ (7)
The logarithmic term in (6) is suppressed in the solution (7) by the factor which vanishes in iteration of the leading solution. It is interesting that the solution does not depend on the gauge parameter $`\alpha _\mathrm{r}`$, in spite of the fact that $`Z_3`$ does. The $`\alpha _\mathrm{r}`$-dependent term in (6) is also suppressed in the solution in the same way that the logarithmic term is. The solution of the compositeness condition should be gauge-independent because it is a relation among physically observable quantities.

Because the above argument relies on $`1/N_s`$ expansion including iteration, the absolute value of the next-to-leading contribution should not exceed that of the leading one. If we apply it to (7), we obtain

$`N_s>22N_c.`$ (8)
The allowed region of $`N_s/N_c`$ by (8) is complementary to that for asymptotic freedom in the gauge theory af , just like in the case of the previous model with fermionic matters cgb2 . When the gauge theory is asymptotically free, the next-to-leading contributions to the compositeness condition are so large that the gauge bosons cannot be composites of the above type. On the contrary, when the theory is asymptotically non-free, the next-to-leading order contributions are suppressed, and the gauge bosons can be interpreted as composite.

Though the complementarity in the scalar matter system is very parallel to that in the fermionic system, they are much different in several points. Because the scalar mass is not protected by any symmetry unlike fermions’, quadratic divergences arise in the scalar self-mass. Though they are not explicit in the dimensional regularization, they should be interpreted to exist from the physical points of view. Then they should be renormalized by unnatural fine tuning. Another characteristic of the scalar matter system is the existence of the self coupling constant $`\lambda `$. At low energies in comparison with the cutoff scale, $`\lambda `$ does not appear in the compositeness condition, as is mentioned above. It could, however, affect our result in various way. For example, it runs with the energy scale, and may blow up at some scale, invalidating the approximation at the order in $`1/N_s`$. On the other hand, the gauge coupling constant $`g`$ remains fixed even for large scale $`\mu `$ (but below the cutoff $`\mathrm{\Lambda }m_\mathrm{r}e^{1/ϵ}`$), because the renormalization group beta function

$`\beta (g_\mathrm{r})\mu {\displaystyle \frac{g_\mathrm{r}}{\mu }}=ϵg_\mathrm{r}+{\displaystyle \frac{1}{6}}N_sg_\mathrm{r}^3I{\displaystyle \frac{11}{3}}N_cg_\mathrm{r}^3I.`$ (9)
vanishes due to solution (7) of the compositeness condition within the present approximation. In general, however, $`g`$ would also blow up together with $`\lambda `$ at the order where $`Z_3`$ depends on $`\lambda `$.

It is straightforward to extend this result to the case where both the fermionic and scalar matters exist. We have only to replace $`N_s`$ in (6) and (7) by $`N_s+4N_f`$, and consequently the complementality still exists because the same replacement rule holds also for the asymptotic freedom. It is tempting to extend it to supersymmetric models. For this case, there are a few problems to be overcome. First we should take into account the effects of the gaugeno, the superpartner of the gauge boson. It gives rise to further quantum corrections through intermediate virtual states, and at the same time it should be a composite since it is a partner of the composite gauge boson. Furthermore it has an additional interaction through gaugeno-scalar-fermion vertex, which would give rise to complexities of diagrams contributing to $`Z_3`$, while asymptotic freedom seems to be free of the effects of these additional interactions. These problems are now under investigation.

We would like to thank Professor H. Terazawa for useful discussions. This work is supported by Grant-in-Aid for Scientific Research, Japanese Ministry of Culture and Science.

Table 1. Contributions to $`Z_3`$ from the diagrams A – H
The diagrams A – H are defined in Fig. 1 in Ref. cgb2  for the model with spinorial matters and in Fig. 1 in this paper for the model with scalar matters. $`I=1/16\pi ^2ϵ`$.

| diagram | with spinorial matter | with scalar matter |
| --- | --- | --- |
| A | $`{\displaystyle \frac{2}{3}}N_fg_\mathrm{r}^2I`$ | $`{\displaystyle \frac{1}{6}}N_sg_\mathrm{r}^2I`$ |
| B, C | $`\left({\displaystyle \frac{13}{6}}{\displaystyle \frac{\alpha _\mathrm{r}}{2}}\right)N_cg_\mathrm{r}^2I`$ | $`\left({\displaystyle \frac{13}{6}}{\displaystyle \frac{\alpha _\mathrm{r}}{2}}\right)N_cg_\mathrm{r}^2I`$ |
| D, D, E | $`{\displaystyle \frac{1}{3}}\alpha _\mathrm{r}N_cN_fg_\mathrm{r}^4I^2`$ | $`{\displaystyle \frac{3}{2}}N_c{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{6}}\right)^l{\displaystyle \frac{N_s^l(g_\mathrm{r}^2I)^{l+1}}{l(l+1)}}`$ |
|  |  | $`+{\displaystyle \frac{1}{12}}\alpha _\mathrm{r}N_cN_sg_\mathrm{r}^4I^2`$ |
| F, G, H | $`{\displaystyle \frac{3}{2}}N_c{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{2}{3}}\right)^l{\displaystyle \frac{N_f^l(g_\mathrm{r}^2I)^{l+1}}{l(l+1)}}`$ | $`{\displaystyle \frac{1}{6}}\alpha _\mathrm{r}N_cN_sg_\mathrm{r}^4I^2`$ |
|  | $`{\displaystyle \frac{2}{3}}\alpha _\mathrm{r}N_cN_fg_\mathrm{r}^4I^2`$ |  |


text:
# Cooperative mixing induced surface roughening in bilayer metals: a possible novel surface damage mechanism

## I Introduction

Low-energy ion bombardment (around some keV) is a widely used technique in surface analysis, sputter cleaning, depth-profiling, micromachining, TEM specimen preparation etc. A tremendous experience on various effects of ion irradiation on solid surfaces has been amassed . There are also several theoretical descriptions of the process especially if we consider homogeneous material. The most visible effect of the ion bombardment, which can be studied easily experimentally (e.g. by means of TEM, STM, AFM) is the change in the topography of the surface. It is generally accepted that the surface morphology produced by subjecting the surface to a given ion fluence is the result of competition between various processes resulting in smoothing and roughening. The roughening can produce somewhat periodic structures and/or rough structures which can be characterized by various scaling lows (a detailed reference of the experimental results are given in ref. ).

Molecular-dynamics (MD) simulations have also demonstrated the existence of the above mechanisms and proved to be capable to describe the development of various surface features during low energy ion sputtering and bombardment .

Three categories of atomic mass transport can be classified when surfaces are subjected to ion beam irradiation : (i) Atoms are ejected from the top layers, usually from the surface and srictly speaking this is what we call ”sputtering” . (ii) The displacement of atoms within their own layers creating vacancies, eroded surface, etc. (iii) Intermixing between layers .

These transport phenomena usually lead to the change of the surface morphology when the thermal spike  is located close to the free surface. It has been shown recently that in this special situation the surface damage mechanisms could be chategorized to four separate mechanisms : ballistic collisions, liquid flow , microexplosions  and coherent displacements . Recent extensive MD and experimental studies showed that the ion beam irradiation leads to surface modifications that are more extensive than just simple sputtering . Ion irradiation of interfacial metal bilayers are much less studied by MD  although the literature of interfaces in crystalline solids is rich .

In this paper we will report on MD simulation of ion irradiation induced alteration in the case of 4 Ti atomic layer/Pt substrate system. This system is specific in that sense that we place the interface close to the free surface ($`10\text{Å}`$) hence the thermal spike effects occur not only close to the free surface but also in the vicinity of the interface. The irradiation conditions were: ion energy 1 keV at grazing angle of incidence ($`83^{}`$). Though the applied conditions can be considered gentle still strong mixing was observed. We will conclude that because of the presence of two unlike atoms enhanced material transport resulting in roughening and mixing occurs. We would like to show that in the case of a sample made of thin layers of different atoms with large difference in the atomic mass, not ion sputtering causes the surface erosion but it is rather induced by interfacial mixing. Morover the mixing of atoms with different atomic mass and the cavitation at the surface is strongly coupled via the ion induced atomic transport normal to the surface leading to a cooperative phenomenon. This process might be a novel surface damage mechanism never reported before.

We irradiate a Ti/Pt system using several ion impacts consecutively. The reason for choosing such an irradiation mode is evident, since practically in all applications a series of ions hit the surface, even for low ion doses . We show through MD simulations of 1 keV Ar ions impacting the Ti surface of the Ti/Pt system that near surface interfacial mixing has a profound effect on the mass transport.

## II Sample preparation and simulation technique

Classical molecular dynamics simulations (MD) were used to simulate the ion-solid interaction using the PARCAS code developed by Nordlund et al. . Here we only shortly summarize the most important aspects. Its detailed description is given in  and details specific to the current system is in a recent communication .

The sample consists of 37289 atoms for the interface (IF) system with 4 Ti top layers and a bulk which is Pt. The lattice constants for Pt is $`a3.92\text{Å}`$ and for Ti $`a2.95`$ and $`c4.68`$ Å. At the interface (111) of the fcc crystal is parallel to (0001) of the hcp crystal and close packed directions are parallel. The interfacial system as a heterophase bicrystal and a composite object of two different crystals with different symmetry are created as follows: the hcp Ti is put by hand on the (111) Pt bulk and various structures are probed and are put together randomly. Finally that one is selected which has the smallest misfit strain prior to the relaxation run. The remaining misfit is properly minimized below $`6\%`$ during the relaxation process so that the Ti and Pt layers keep their original crystal structure. The corresponding Ti-Ti and Pt-Pt interatomic distances are $`2.89`$ and $`2.73`$ Å at the interface. The Ti (hcp) and Pt (fcc) layers at the interface initially are separated by $`2.8`$ Å and allowed freely to relax during the simulations. The variation of the equilibrium Ti-Pt distance within a reasonable $`2.63.0`$ Å interval in the interatomic potential does not affect the final results significantly. We find the average value of $`d2.65`$ Å Ti-Pt distance in the various irradiation steps and also in the nonirradiated system after a careful relaxation process. We believe that the system is properly relaxed and equilibrated before and between the irradiation steps. We found no vacancies at the interface or elsewhere in the system after the relaxation procedure. We visually checked and found no apparent screw dislocations, misorientations or any kind of undesired distorsions (in general stress-generator dislocations) at the interface or elsewhere in the system which could result in stressed induced rearrangements in the crystal during the simulations. The temperature scaled softly down towards zero at the outermost three atomic layers during the cascade events . The simulations are repeated in a much larger system consisting of $`72000`$ atoms in order to rule out any kind of size effects. The simulation cell is of particular interest because it has been shown  that at the applied sputtering conditions (1 keV Ar+ with grazing angle of incidence) the thermal spike is typically formed at the interface and thus the emerging processes can be easily studied.

The initial velocity direction of the impacting atom was $`7`$ degrees with respect to the surface of the crystal. Further details are given in Ref. . A series of up to $`10`$ ion impacts are employed in order to study the effects of multiple ion irradiations, e.g. to simulate high dose effect. Between impacts the simulation time was $`20`$ ps. This time was found to be sufficient to reach a relaxed system after each ion impact (e.g. the temperature approaches zero). Thus direct interaction between the cascade processes of various irradiation steps are excluded. It should also be remarked that the damage rate is extremely high due to the high dose effect and therefore some recovery process might occur between the irradiation steps. Although we did not find significant recrystallization effects up to several hundred ps. Thermally activated process might induce, however, further recovery in the system. We have done relaxation simulations on the irradiated sample at 300 K and found no recrystallization of the amorphised top layers in the irradiated zone up to several hundred ps. In any case, we do not think that the alloy phase would phase separate either at room temperature or on a longer time scale. We will demonstrate in the Discussion section that the the alloy phase is energetically competitive with the phase separated system (based on the cohesive energies).

To obtain a representative statistics, the impact position of the incoming ion is varied randomly within a $`10\times 10`$ Å2 area. We find considerable variation of mixing as a function of the impact position of the recoil. Using a $`20\times 20`$ Å irradiation area (with 80000 atoms or more) we checked that a further increase of the irradiated area does not affect the final results heavily, except that the larger the area, the larger the number of incoming ions needed to get the same extent of interfacial mixing (IFM) and surface roughening. This is because individual ions with impact positions far from each other do not result in the high dose effect ($`10^{15}`$ Ar+ ions/cm2) we are interested in in this article. Impact positions with sufficient proximity to each other (that is within the characteristic area $`5\times 5`$ Å2, $`10^{15}`$ ion/cm2) lead to the phenomenon called mixing induced enhancement of surface roughening (MIESR). The overlap of these individually eroded surface areas might lead to the macroscopic phenomenon of sputter removal when a sufficiently high dose of ions are bombarded to the surface in that case when sputtering is the only active mechanism for surface damage. However, the final morphology of the sputtered surface might also be affected by MIESR if sufficiently large interfacial mixing occurs close to the free surface.

We saw no further recovery process after 20 ps simulations although we have extended several runs until 500 ps. Thermally activated processes do not also play a significant role up to moderate temperatures (room temperature or so). Also we set the quenching rate to a value (5 K/Fs) which seems to be a reasonable choice. Reasonable variation of the quenching rate does not affect the cooling of the system and does not lead to the recovery of the original crystal order in the upper layers.

During the simulations we observe extensive depression of the bombarded surface with respect to the nonsputtered surface. This depression did not result in a flat bottom, rather a rough surface had formed. We estimate the extent of cavitation by counting the number of missing atoms in the cavity using the same procedure introduced for vacancies .

Vacancies were recognized in the simulations using a very simple analysis. A lattice site with an empty cell was considered to be a vacancy, where the diameter of the cell sphere is set to $`2\text{Å}`$ (the radius of a sphere around the relocated atom around its original position) which is about $`60\%`$ of the average atomic distance in this system. We find this criterion is suitable for counting the number of vacancies.

We used a many body potential, the type of an embedded-atom-method given by Cleri and Rosato , to describe interatomic interactions. This type of a potential gives a very good description of lattice vacancies, including migration properties and a reasonable description of solid surfaces and melting . Since the present work is mostly associated with the elastic properties, melting behaviors, surface, interface and migration energies, we believe the model used should be suitable for this study. The Ti-Pt interatomic potential of the Cleri-Rosato  type is fitted to the experimental heat of mixing of the corresponding alloy system. Further details are also given elsewhere .

An atom is labeled mixed if it moved beyond the interface by more then $`d/21.4\text{Å}`$ along the $`z`$-axis, where $`d`$ is the interfacial Ti-Pt distance. The broadening at the interface (or interfacial mixing, IFM) is calculated as follows:

$$b=N_{mix}^1\underset{i=1}{\overset{N_{mix}}{}}\mathrm{\Delta }z_i,$$
(1)
where $`N_{mix}`$ is the number of mixed atoms and $`\mathrm{\Delta }z_i`$ is the distance of the mixed atom $`i`$ from the interface ($`z`$-axis depth position). $`b`$ is a useful quantity for comparison with measurements only in that case when $`\mathrm{\Delta }z_i`$ is larger then the interlayer distance. The reason is that the measured broadening is dominated by those intermixed atoms which move beyond the interface by more than the inter-layer distance (2nd layer mixed atoms). Using this condition we get a reasonable agreement with Auger depth profiling .

## III Results and discussion

The detailed analysis of the mixing process in the Ti/Pt interfacial system after a single ion irradiation process has been given recently . Here we restrict ourselfs to the study of the effect of multiple ion impacts on mixing as well as on surface erosion.

In FIGs (1)- (2) the $`xz`$ cross-sectional and $`xy`$ top plots of the Ti/Pt system are shown after the $`5`$th irradiation. It is evident that a considerable amount of material has been removed from the irradiated area. Since the height of the non-irradiated surface did not change apart from some adatoms the material removal produces a feature like a crater. It should be added, however, that the feature formed has a rough surface. The volume of the crater is much larger than which is expected considering only sputtered atoms. The number of sputtered atoms is less by an order of magnitude then the number of intermixed Ti atoms which are pushed into the Pt bulk. These atoms cause the considerable interface broadening visible in FIGs (1) and (4). Indeed, the majority of removed atoms are the intermixed Ti atoms (injected to the bulk). In FIGs (3) we show the corresponding pure components after the $`5`$th bombardment and see no apparent surface morphology development. This also shows that the degree of the material removal by sputtering is much weaker than the amount of the missing atoms in the case of the multilayered material.

Since the difference between the pure material and the Ti/Pt system is that in the case of the latter a two component system forms due to ion sputtering, we should conclude that the presence of the second component seriously affects the mass transport and thus most likely leads to morphology development as well. Let us check more closely the correspondence between IFM and mass transport (which is measured by the volume of the apparent crater). Our primary intention is then to point out that the enhancement of surface cavitation in Ti/Pt is due to interfacial mixing. MIESR is strongly connected to a two-way material transport (with counterflow motion) which might be understood as a radiation-induced flow of Pt atoms to the surface and the injection of Ti to the bulk.

We see the strong enhancement of IFM due to the repeated irradiation of the surface in FIG (4) where the broadening at the IF is plotted against the number of irradiation steps (upper figure). Anomalously large intermixing in metal bilayers have recently been seen as well by Buchanan et al. . Surprisingly, we see the saturation of mixing (broadening) at the 5th step. The magnitude of mixing at saturation shows no dependence on heat of mixing as it was obtained in the first irradiation step . Hence MIESR is not governed by thermal spikes.

In the lower part of FIG (4) the number of ”missing” atoms in the top layer (vacancies) is plotted against the number of irradiation steps. The number of missing atoms is proportional roughly to the volume of the surface cavity which are mostly injected to the bulk. Interestingly these features coincide with each other: the number of intermixed Ti atoms (the ”missing” atoms) saturates also at the $`5`$th step.

The volume of the cavity rapidly increases in the first three irradiation steps and a crater similar to that of reported in earlier studies by Bringa et al. at higher impact energies with single ion bombardment, forms . In these cases cratering is attributed to the high deposited energy close to the free surface at certain nonchanneling directions. In our case, however, after repeated irradiation of the surface, there seems to be a strong correlation between interfacial intermixing and material removal from the surface (material injection to the bulk). Pt atoms intermixed to the top layers, however, do not occupy the empty sites in the cavity. Rather they form a dense phase at the surface of the cavity. The concentration of atoms in the irradiated region ($`10\times 10\times 10`$ Å3) increases from the initial $`6.1\times 10^{22}`$ atom/cm3 to $`6.8\times 10^{22}`$ atom/cm3 in the $`5`$th irradiation step. The atomic concentration $`n`$ in TiPt alloy is $`7.2\times 10^{22}`$ atom/cm3, therefore the increase of $`n`$ in the irradiated region indicates the formation of a broad alloy phase at the interface. However, the saturation of the curves in FIGs (4) indicates that the simultaneouss process of intermixing and surface cavitation is limited. A possible control factor can be the increase in the number of the undercoordinated surface atoms which is always unfavored. Another factor could be the saturation of the broadened interface. No more intermixing of Ti atoms is observed possibly due to the denseness of the alloy phase formed in the original Pt bulk which does not permit the further transport of Ti atoms to the Pt bulk. Anyhow the maximal broadening of the IF seems to be limited to $`b10\text{Å}`$. This restriction on broadening seems not to be affected by chemical forces since it was found that the heat of mixing did not affect the IFM .

In order to understand more fundamentally the mechanism of MIESR we vary the atomic mass of the top layer constituent ($`m_{top}`$), keeping the same potential which can be done easily in MD simulations . According to our expectations we find the suppression of IFM and surface cavitation when $`m_{top}/m_{Pt}1`$.

We plot in FIG (5) the number of Ti atoms which was removed from the top layers by MIESR (cavitation) as a function of the atomic mass ratio obtained from simulations in the third irradiation step. In this step we obtained extensive mass transport especially for lighter atoms ($`m_{top}211`$ g/mol, $`\delta 0.010.05`$, see also FIG (1)). Notable feature is in lower FIG (1) the enhanced cavitation with small mass ratio in the third irradiation step.

It can obviously be seen that the lighter the atoms of the top layers, the larger the extent of the mass transport. At the same time we see no obvious effect of heat of mixing on IFM . It is therefore quite probable that IFM and surface roughening is governed primarily by the atomic mass difference between the constituent atoms of various layers. Mass transport from the surface layer is strongly enhanced for multilayers with light atoms in the thin top layer. This finding supports our previous conclusion  that the ballistic mechanism governs the IFM.

This is interesting because we are not aware of similar results where the mass transport is strongly coupled to interfacial mixing. Moreover the bulk process (mixing) and surface roughening are coupled via the enhanced material transport normal to the surface leading to a cooperative phenomenon. The importance of this finding is to point out that the surface of materials which include more components close to the surface might behave in a different way than the pure components. The surface erosion is governed not primarily by sputtering but rather by mixing. A specific feature of this system is that we found no effect of chemical forces (heat of mixing) on IFM . The process is governed by the large difference of atomic mass between Ti and Pt.

In order to understand the energy saldo of mixing we plot in FIG (6) the cohesive energy (binding energy) of various atoms as a function of their depth position. The data points are collected from the $`10`$th irradiation step. The atomic cohesive energy $`E_{coh}`$ is averaged with that of every atom’s nearest neighborsrs to reduce the atom-level scatter in the values. It can obviously be seen that the cohesive energy of Pt atoms decreases when the atoms move towards the surface while Ti atoms gain energy when mixed to the Pt bulk. The energy gain per Ti atom is $`1`$ eV on average while approximately the same energy loss occurs for the intermixed Pt atoms. Above the depth position of $`30`$ Å the widening of the cohesive energy ”window” can also be seen for Pt. This is due to the formation of an alloy phase at the broadened interface. Close to the free surface no features of alloy phase formation can be seen, e.g. the $`E_{coh}`$ values are rather close to the atomic value of pure Ti. Also the $`E_{coh}`$ of the intermixed Pt atoms decreases at the surface which is due to their decreased coordination number ($`<12`$) with respect to the bulk value of $`12`$.

Finally it should also be mentioned that the ion induced surface cavity growth might have many common features with the growth of cavities found on boundaries in stressed polycrystals or with cracking and in general with various fracture processes at heterophase interfaces . This is because we do not rule out the possibility of a cavity growth mechanism induced by a stress field generated by the thermal spike. This is, however, beyond the scope of the present article and needs further studies.

## IV Conclusions

We have shown that the mass transport for the studied double layer system is different from that of the corresponding pure elements. We have also shown that the observed apparent cavitation is in correlation with the interfacial mixing and it increases with increasing mass difference of the constituents. The most important consequence of this finding is that a reliable theory describing the surface morphology development of a system containing two components should also account for the mixing near the surface. Presently available theories are missing this point. This novel surface damage mechanism might be applicable for those multicomponent systems (e.g. thin films) where the material is subjected to low energy ion-irradiation and the interface is close to the free surface. The mechanism is regressively reliable for systems with mass ratio approaching $`1`$.

## V acknowledgment

This work is supported by the OTKA grants F037710 and T30430 from the Hungarian Academy of Sciences and from EU contract no. ICAI-CT-2000-70029.


text:
# Regular and Random Magnetic Resonance Force Microscopy Signal with a Cantilever Oscillating Parallel to a Sample Surface

## I Introduction

Inspired by the recent progress in magnetic resonance force microscopy (MRFM) based on oscillating cantilever-driven adiabatic reversals (OSCAR) technique , we consider here the theory of the OSCAR MRFM. In the OSCAR MRFM the oscillating cantilever with the ferromagnetic particle on its tip produces oscillating dipole field along the direction of the permanent external magnetic field. The rf magnetic field is applied in the plane perpendicular to the permanent external field. In the rotating system of coordinates the effective magnetic field changes its direction with the frequency of the oscillating cantilever. If the conditions of adiabatic reversal are satisfied for a spin of a sample then the spin will follow the direction of the effective magnetic field. The region of the sample where the effective magnetic field periodically reverses its direction is called the resonant slice. The spins of the resonant slice experience cyclic adiabatic reversals and produce the back effect on the cantilever: they cause the shift of the cantilever vibration period. In this paper we call this shift ”the MRFM signal”. In previous works  we have studied the “perpendicular oscillations”: the cantilever tip oscillates perpendicular to the surface of a sample. In its advance to a single-spin detection the OSCAR MRFM technique relays on the “parallel oscillations”, which allow one to reduce a spacing between the cantilever tip and the sample. Thus, the OSCAR MRFM theory should be extended to the parallel oscillations setup. This paper represents such an extension. We show that the thermal high frequency cantilever vibrations cause approximately exponential decay of the initial regular MRFM signal. However the MRFM signal does not disappear. It transforms into the random MRFM signal. The Fourier spectrum of the random signal has a characteristic peak, which can be used for the identification of the signal.

## II Equations of motion

We assume that a ferromagnetic spherical particle is attached to the cantilever tip and oscillates along the $`x`$-axis which is parallel to the surface of the sample. (See Fig. 1.) In the equilibrium position the center of the ferromagnetic particle is at the origin. The permanent external magnetic field $`\stackrel{}{B}_{ext}`$ points in the positive $`z`$-direction, and the rf magnetic field $`\stackrel{}{B}_1`$ rotates in the $`xy`$-plane. The motion of the center of the ferromagnetic particle can be described by the equation

$$\ddot{x}_c+x_c+\dot{x}_c/Q=f(\tau ),$$
$`(1)`$
where $`x_c`$ is the dimensionless coordinate of the center (in terms of the oscillation amplitude of the cantilever, $`x_0`$), $`Q`$ is the effective quality factor of the cantilever, and we use the dimensionless time $`\tau =\omega _ct`$, where $`\omega _c`$ is the cantilever fundamental frequency. The dimensionless force $`f(\tau )`$ is produced by all magnetic moments of the “resonant slice” of a sample:

$$f(\tau )=\underset{k=1}{\overset{\mathrm{N}}{}}\eta _k\mu _{kz},$$
$`(2)`$
$$\eta _k=\frac{\mu _0}{4\pi }\frac{3m\mu }{k_cx_0^5}\frac{(\stackrel{~}{x}_k/\stackrel{~}{r}_k)(5z_k^2/\stackrel{~}{r}_k^21)}{\stackrel{~}{r}_k^4}.$$

Here N is the number of magnetic moments in the resonant slice, $`\stackrel{}{\mu }_k`$ is the $`k`$-th magnetic moment (in terms of its magnitude $`\mu `$ which is the same for all magnetic moments), $`k_c`$ is the effective spring constant of the cantilever, $`m`$ is the magnetic moment of the ferromagnetic particle, $`\stackrel{~}{x}_k=x_kx_c`$, $`\stackrel{~}{r}_k=[(x_kx_c)^2+y_k^2+z_k^2]^{1/2}`$, $`x_k`$, $`y_k`$, $`z_k`$ are the coordinates of the $`k`$-th magnetic moment in terms of the oscillation amplitude $`x_0`$. We assume in (2) that only $`z`$-component of $`\stackrel{}{\mu }_k`$ influences the motion of the cantilever.

The motion of the $`k`$-th magnetic moment of a sample in the rotating frame can be described by the equations:

$$\dot{\mu }_{kx}=\mathrm{\Delta }_k\mu _{ky},$$
$`(3)`$
$$\dot{\mu }_{ky}=\mathrm{\Delta }_k\mu _{kx}\epsilon \mu _{kz},$$

$$\dot{\mu }_{kz}=\epsilon \mu _{ky},$$

$$\mathrm{\Delta }_k=(\gamma B_{ext}\omega )/\omega _c+\frac{\mu _0}{4\pi }\frac{\gamma m}{\omega _cx_0^3}\frac{3z_k^2/\stackrel{~}{r}_k^21}{\stackrel{~}{r}_k^3}.$$

Here $`\gamma `$ is the gyromagnetic ratio of the spins in the sample, $`\omega `$ is the frequency of the rf field, $`\epsilon =\gamma B_1/\omega _c`$, and we neglect the transverse components of the dipole field produced by the ferromagnetic particle.

The thermal vibrations of the cantilever tip with the frequencies close to the Rabi frequency $`\gamma B_1`$ cause the deviation of the magnetic moment from the direction of the effective magnetic field and the decay of the regular MRFM signal . To take into account this effect we make a substitution in the expression for $`\mathrm{\Delta }_k`$ in (3):

$$x_cx_c+\delta x_c,$$
$`(4)`$
where $`\delta x_c`$ describes the thermal vibrations of the cantilever tip

$$\delta x_c=\underset{n}{}\frac{2a_n}{x_0}\mathrm{cos}(\omega _n\tau +\mathrm{\Psi }_n).$$
$`(5)`$
Here $`a_n`$ is the amplitude of the thermal vibrations for a cantilever mode $`n`$, $`\omega _n`$ is the eigenfrequency of the mode $`n`$ (in the units of $`\omega _c`$), $`\mathrm{\Psi }_n`$ is its phase. In our simple model we assume that the amplitude $`a_n`$ can be found from the equipartition theorem

$$m_c\omega _n^2a_n^2=2k_BT,$$
$`(6)`$
where $`m_c`$ is the cantilever mass, and the phase $`\mathrm{\Psi }_n`$ is the random function of time. The factor 2 in (5) appears as the amplitude of the tip vibrations for a uniform cantilever is twice the amplitude of the mode . (For a non-uniform cantilever this factor depends on $`n`$ .) We ignore the influence of the thermal vibrations on the value $`\eta _k`$ in (2).

## III Numerical simulations

In our numerical simulations we used the following parameters: the fundamental cantilever frequency $`\omega _c/2\pi =7`$kHz, the effective spring constant $`k_c=10^4`$N/m, the effective quality factor $`Q=5\times 10^4`$, the amplitude of the cantilever vibrations $`x_0=10`$nm (to model the action of the feedback technique in the OSCAR MRFM, our computer algorithm increased the value $`x_c`$ to “1” every time when the cantilever tip passed the maximum value of $`x_c`$), the radius of the ferromagnetic spherical particle $`R=200`$nm, the magnetic moment of the sphere $`m=2.5\times 10^{14}`$J/T, the distances from the bottom of the sphere to the sample surface and to the center of the resonant slice $`d_1=220`$nm and $`d_2=300`$nm (see Fig. 1), the amplitude and the frequency of the rf field $`B_1=3\times 10^4`$T and $`\omega /2\pi =3`$GHz, the magnetization of the sample is $`0.9`$A/m. The resonant slice boundaries were found from the condition $`\mathrm{\Delta }_k(x_c=\pm 1)=0`$.

Note that for two spins with coordinates ($`x,y,z`$) and ($`x,y,z`$) change of the $`z`$-component of the dipole field caused by the cantilever displacement has an opposite sign. Let, for example, the $`z`$-component of the effective field for these two spins is zero when the ferromagnetic particle is at the origin ($`x_c=0`$). If $`x_c0`$ the effective magnetic field on the first spin is opposite to that on the second spin.

If the initial direction of the two spins relative to the external magnetic field is the same, these spins will have an opposite direction relative to the local effective field. Thus, the two spins induce the MRFM signal of the opposite sign. If the spins are uniformly distributed in the resonant slice, the net MRFM signal disappears. That is why in our simulations we assume that spins occupy only the resonant semi-slice $`x>0`$. The cross-sectional area of the semi-resonant slice in the $`xz`$-plane is shown in Fig. 2.

The magnetic moments in our model are randomly distributed inside the semi-resonant slice. We assume that initially all magnetic moments point in the positive $`z`$-direction, and the cantilever tip is in its right end point, $`x_c=1`$.

In our computer simulations the phases $`\mathrm{\Psi }_n`$ in (5) were changed randomly between $`0`$ and $`2\pi `$, with random time intervals between two successive “jumps”. The time interval between the phase jumps was taken randomly between $`8.3T_R`$ and $`14T_R`$, where $`T_R=2\pi /\gamma B_1`$ is the Rabi period. In Eq. (5), we took into consideration 25 cantilever modes in the vicinity of the Rabi frequency.

Next, we describe the results of our simulations. The main contribution to the MRFM signal is associated with the magnetic moments near the surface of the sample. The regular MRFM signal decays approximately exponentially. The regular MRFM signal is followed by the non-dissipating random signal. The random signal has its characteristic signature: the major peak in its Fourier spectrum.

Figs 3-5 illustrate the results of our computer simulations. Fig. 3 shows the decay of the regular signal $`\mathrm{\Delta }T/\mathrm{\Delta }T_0`$, where $`\mathrm{\Delta }T`$ is the shift of the cantilever period caused by the magnetic force, $`\mathrm{\Delta }T_0`$ is the initial value of $`\mathrm{\Delta }T`$. Fig. 4 demonstrates the random signal which follows the regular signal shown in Fig. 3. Finally, Fig. 5 shows the Fourier spectrum of the random signal with its characteristic peak.

## Conclusion

We have studied theoretically the OSCAR MRFM signal for the case when the cantilever tip oscillates along the axis parallel to the surface of the sample. In this case, the main contribution to the signal is associated with spins located in the part of the resonant slice near the surface of the sample. The main features of the OSCAR MRFM signal are the following: the regular (approximately exponential) decay of the MRFM signal is followed by the non-dissipating random signal. The Fourier spectrum of the random signal has a characteristic peak which can be used for identification of the signal.

## Acknowledgments

This work was supported by the Department of Energy under the contract W-7405-ENG-36 and DOE Office of Basic Energy Sciences, and by the DARPA Program MOSAIC.


text:
# HD 161796 (IRAS 17436+5003)

## 1. The Object

HD 161796 (IRAS 17436+5003) is a low metallicity oxygen-rich post-AGB star of F3Ib located at high Galactic latitude (e.g., Hrivnak et al. 1989). Due to the star’s relative proximity ($`1`$ kpc), its circumstellar shell has been relatively well studied at various wavelengths in order to understand its mass loss history.

## 2. The Outer Shell

The circumstellar shell around HD 161796 was marginally resolved by IRAS at 100 $`\mu `$m, and the IRAS map seems to show a spherically symmetric shell. Using ISO linear scan observations across the shell at 90 and 160 $`\mu `$m, Speck et al. (2001) have recently shown that the shell is indeed extended out to about 400<sup>′′</sup>. Even though the ISO scan was done only in one direction, it seems reasonable that this cold ($`20`$ K) dust shell possesses generally spherically symmetric structure. This parsec-sized, nearly spherical cold dust shell corresponds to the part of the circumstellar shell that was created by an ancient AGB mass loss roughly $`2\times 10^5`$ years ago (Fig. 1, left).

## 3. The Inner Shell

While the outer shell represent the old AGB mass loss history, the inner shell is a manifestation of the most recent mass loss history that may provide clues for the latest structure formation in the circumstellar shell. By directly probing the distribution of warm ($`100200`$ K) dust grains, Skinner et al. (1994) have revealed a two-peaked structure in the emission core in their deconvolved 12 $`\mu `$m image. Most recently, Gledhill et al. (2003) have taken N and Q band images of the object using Gemini-N and resolved the two-peaked structure in the mid-IR emission core (Fig. 1, right). The two-peaked mid-IR structure is typically interpreted as limb-brightened edges of a dust torus (or an equatorially enhanced dust distribution) embedded in the innermost region of the circumstellar shell.

The dust shell structure can also be investigated indirectly through dust-scattered star light in the optical and near-IR. As part of their HST imaging survey of PPNs, Ueta et al. (2000) have detected an elliptically elongated shell surrounding the mid-IR core (Fig. 1, middle). These non-spherically symmetric structure of the shell seem to have developed during the last $`500`$ years before the mass loss was terminated nearly 300 years ago (Meixner et al. 2002).

## 4. The Overall Shell Structure

These images of the shell at various wavelengths probe different parts of the shell, and hence, can be thought of as “time slices” that would help to reconstruct the entire AGB mass loss history of the star. Based on the shell structure we see at different epochs, the overall shell structure seems to be explained by a “layered shell” model, in which the shell consists of three layers that represent different epochs of mass loss (Fig. 2, left).

The outermost layer of cold dust corresponds to the part of the shell created by the earliest AGB mass loss, in which the wind was spherically symmetric. The innermost layer of warm dust, on the other hand, corresponds to the part of the shell generated by equatorially enhanced (or axisymmetric/toroidal) mass loss at the end of the AGB phase. In addition, the mid-layer seen in the dust-scattered star light corresponds to one particular epoch during which the shell geometry transforms from spherical to axisymmetric (toroidal) symmetry.

2.5-dimensional dust radiative transfer calculations have been performed to test this layered shell model (Meixner et al. 2002; Ueta et al. 2003). The best-fit model has successfully reproduced the spectral energy distribution from UV to far-IR and characteristic morphologies at various wavelengths from the optical through mid-IR (Fig. 2, right).

## 5. The “Twist” in the Shell

Although the overall shell structure can be understood by the layered shell model consisting of three distinct structures, the model is not enough to explain the detailed structure of the shell, especially the apparent deviation from axisymmetry. If we define the axis of the shell based on the toroidal structure in the mid-IR core, the axis would have the position angle of roughly 30. However, the elliptical structure of the dust-scattered reflection nebulosity would delineate the axis of symmetry oriented at roughly 0 position angle. By comparing these two “time slices” (i.e., mid-IR and optical images representing different epochs of mass loss), this apparent “shift” of the shell orientation seems to occur rather abruptly.

However, recent mid-IR observations have captured dust emission of the innermost core and of the low emission halo regions in a single image. So, the axis of symmetry now appears to shift continuously in the counter-clockwise direction over the course of the star’s mass loss history (Gledhill et al. 2003; Fig. 3, left). Interestingly, if one follows the density distribution further away from the central star probing even colder parts of the shell through CO gas in the mm-wavelengths, the twisting of the shell seems to be continuing for quite some time (roughly 1000 years; Fong et al. 2003; Fig. 3, right).

## 6. Summary

As we have seen, imaging of the circumstellar shells at various wavelengths is equivalent to taking “snapshots” of the mass distribution at various epochs of the mass loss history. This method has been quite effective to observationally establish the geometrical transition of the circumstellar shells. Recent imaging of HD 161796 at the mid-IR has probed the outer shell of cold dust and has suggested that the axis of symmetry rotates continuously while the shell morphology assumes more equatorially enhanced structure.

Therefore, it is very interesting to increase the “time resolution” of the early to intermediate phases of the mass loss history represented by rather cold (around $`20100`$ K) part of the dust shells (corresponding to the mass loss epochs on the order of $`10^3`$ years ago). With far-IR to sub-mm telescopes/missions such as SIRTF, SOFIA, ASTRO-F, VLTI, Herschel, and SMA, we will have a plenty of opportunity to observe these critical regions of the circumstellar shells at higher resolution and with higher sensitivities. We can then follow the early phases of the mass loss history and address how spherically symmetric circumstellar shells develop the equatorially-enhanced structure. This line of research may also provide clues to explain how elliptical and bipolar PPNs assume their respective structure by comparing their early mass loss history.

## References

Fong, D., Meixner, M., Sutton, E. C., Zalucha, A., & Welch, W. J. 2003, ApJ, submitted

Gledhill, T. M., & Yates, J. A. 2003, MNRAS, 343, 880

Hrivnak, B. J., Kwok, S., & Volk, K. 1989, ApJ, 346, 265

Meixner, M., Ueta, T., Bobrowsky, M., & Speck, A. K. 2002, ApJ, 571, 936

Skinner, C. J., Meixner, M., Hawkins, G. W., Keto, E., Jernigan, J. G., & Arens, J. F. 1994, ApJ, 423, L135

Speck, A. K., Meixner, M., & Ueta, T. 2001, in Post-AGB Objects as a Phase of Stellar Evolution, eds. R. Szczerba & S. K. Górny (Dordorecht: Kluwer), 333

Ueta, T., & Meixner, M. 2003, ApJ, 586, 1338

Ueta, T., Meixner, M., & Bobrowsky, M. 2000, ApJ, 528, 861