text:
# Looking back at superfluid helium

## The discovery of superfluidity

In january 1938, J.F. Allen and A.D. Misener Allen1  published the experimental evidence that the hydrodynamics of liquid helium was not classical below 2.2 K. In the same issue of Nature, Kapitza introduced the word “superfluid” to qualify this anomalous behavior Kapitza . Soon after this discovery, F. London suggested that superfluid helium forms a macroscopic liquid matter wave, as a consequence of “Bose-Einstein condensation” (BEC). 65 years later, it is generally accepted that London was right, but it happenned to be very difficult to prove. On the contrary, when BEC was discovered in cold alkali vapors (in 1995), the experimental evidence was immediately very clear. The superfluidity of these cold gases was soon demonstrated as well. It is thus interesting to look back at some properties of superfluid helium now. As we shall see, not everything is completely understood in the properties of superfluid helium. In this short review, I have selected only a few properties of superfluid helium. For a complete description, one can look at the historical book written by J. Wilks Wilks , or the one by P. Nozières and D. Pines on Bose liquids Nozieres . Here, I only wish to recall some aspects of the early history of this subject in order to show that fundamental questions came up very soon in this field. Then, I whish to briefly discuss three particular topics:

1- boiling and evaporation because this is what makes superfluid helium look different from normal helium in a dewar.

2- the mechanisms which determine the “critical velocity”

3- the existence of a condensate in superfluid helium and the pressure variation of the superfluid transition temperature.

In the years 1928-32 in Leiden, W.H. Keesom and his co-workers had already realized that liquid helium exhibited surprising properties below about 2.2 K Keesom1927 . They had found a peak in the temperature variation of the specific heat at this temperature. The shape of this peak resembling the greek letter lambda, Keesom called the superfluid transition temperature the lambda point or “lambda temperature $`T_\lambda `$”. He also called “helium I” the liquid above $`T_\lambda `$ and “helium II” the liquid below $`T_\lambda `$. In his laboratory in Leiden, Keesom had also discovered that helium II was able to flow through very tiny pores as soon as in 1930 Keesom1930 . In 1932 in Toronto, J.C. McLennan, H.D. Smith and J.O. Wilhelm then discovered that liquid helium ceased boiling when cooled down through $`T_\lambda `$ McLennan . This was soon attributed to an exceptionnally high thermal conductivity, as evidenced by successive experiments by W.H. Keesom again Keesom1936  and by J.F. Allen in Cambridge Allen1937 . J.F. Allen had been hired with R. Peierls to replace Kapitza there, because Kapitza had been forced by Stalin in 1934 to stay in Moscow and to quit his research position in CambridgeSFP . Since the high thermal conductivity was attributed to some kind of convection or turbulence, researchers looked at the flow of liquid helium. In 1935 in Toronto, J.O. Wilhelm, A.D. Misener and A.R. Clark had measured the viscosity of liquid helium with a torsion pendulum and found that it decreased sharply below $`T_\lambda `$ Wilhelm . When A.D Misener moved to Cambridge in order to prepare his doctorate with J.F. Allen, they started together a systematic study of the flow of liquid helium through capillaries. Classical fluids obey the Poiseuille law: the flow rate is proportional to the pressure difference accross the capillary and to the fourth power of the capillary radius. Instead of this, Allen and Misener found that, below $`T_\lambda `$, the flow rate was not only high, it was independent of the pressure and independent of the capillary radius which they had changed by a factor 50! Clearly, this liquid was not classical. In a sense, this was also what Kapitza claimed in his article Kapitza , since he introduced the word “superfluid” in connection with the already know “superconductivity”. However, Kapitza did not justify his remarkable intuition that superfluidity and superconductivity were related phenomena, and his 1938 article contained no quantitative measurements. He performed remarkable experimental measurements during the following years.

In february 1938, J.F. Allen and H. Jones published another estonishing discovery. They had found a remarkable thermomechanical effect AllenJones . When heating superfluid helium on one side of a porous medium or a thin capillary, the pressure increased sufficiently to produce a little fountain at the end of the tube which contained the liquid. The “fountain effect” was another spectacular phenomenon which was impossible to understand within classical thermodynamics. This is really what triggered Fritz London’s thinking London . London new that, in this liquid, there were large quantum fluctuations responsible for the large molar volume, and he proposed that the superfluid transition at $`T_\lambda `$ was a consequence of Bose-Einstein condensation (BEC). 65 years later, the connection between superfluidity and BEC is still a matter of debate and study. Before discussing this (in Section III), I wish to say a few words about boiling and evaporation (section I), and , in Section II, I briefly review the problem of the critical velocity, in connection with Landau’s “two fluid model” Landau , the existence of “rotons” and quantized vortices. One of my motivations in section II is to allow a comparison with observations of a critical velocity and vortices in gaseous superfluids.

## I Boiling, evaporation and cavitation

The easiest way to see superfluidity is to look at liquid helium in a dewar, while pumping on it so as to cool it below $`T_\lambda `$ = 2.2 K (see Fig. 1). When crossing $`T_\lambda `$, the liquid stops boiling. This is because the thermal conductivity of liquid helium has suddenly increased, so that the temperature is very homogeneous. The walls are no longer warmer than the surface. They no longer provide efficient nucleation sites for bubbles. As a result, superfluid helium evaporates from its free surface, instead of showing bubble nucleation on hot defects. Superfluid helium looks quiet, and it is tempting, although not rigorous, to take this as an illustration of quantum order getting established. At least this shows a striking reality: a liquid made of very simple atoms can exist in two different states! It happens that I worked on the evaporation of superfluid helium Balibar78  so that i would like to add a few comments on the evaporation of superfluid helium, an interesting property which is rarely mentionned. We have verified the prediction by P.W. Anderson that, at low enough temperature, this evaporation is a quantum process which is analogous to the photoelectric effect Anderson . A photon incident on a metal surface can eject one electron elastically: the electron kinetic energy is the difference between the photon energy and the binding energy of the electron to the metal. Similarly, we have shown that atoms can be ejected from superfluid helium with a kinetic energy which is close to 1.5 K, the difference between the energy of incident rotons (whose minimum is 8.65 K) and the binding energy of atoms to the liquid, which is nothing but the latent heat of evaporation per atom (7.15 K). This phenomenon was later called “quantum evaporation” by A.F.G. Wyatt who made a long quantitative study of it. A few open issues remain in this subject, such as the exact measurement and calculation of evaporation probabilities by phonons or rotons as a function of energy, momentum and incidence angle, or the nature of possible inelastic processes as well.

Furthermore, one should not believe that boiling, more generally cavitation or bubble nucleation, is impossible in superfluid helium. On the contrary, this is another interesting subject on which we also worked Balibar02 . A strong local heat input can generate bubbles in superfluid helium, of course. Moreover, one can apply a large stress, that is a large negative pressure to liquid helium at low temperature and also look for the nucleation of bubbles. This is a way to test the internal cohesion of this liquid. We have observed this phenomenon near -9.5 bar, the stability limit of liquid helium under stress, also called the “liquid-gas spinodal limit”. It allowed us also to obtain some information on superfluidity at reduced density as we shall see briefly in Section III.

## II The critical velocity: rotons and vortices

The vanishing of boiling is spectacular but the existence of flow without dissipation is of course the more fundamental property which signals superfluidity. In may 1938, one month only after London had published his suggestion of a Bose-Einstein condensation, L. Tisza introduced a “two fluid model” which qualitatively explained the observations. He considered that “the atoms belonging to the lowest energy state do not take part in the dissipation. Thus the viscosity of the system is entirely due to the atoms in the excited states”. He understood that this liquid could be considered as a mixture of two components: a superfluid component which carries no entropy and has zero viscosity, plus a normal component, made of atoms in excited states, which is viscous and conducts heat. The ratio of one to the other component was determined by the temperature: at zero temperature there is only the superfluid and at the lambda transition there is only the normal component. The important point was to understand that these two components could move independently. This very pure liquid had two independent velocity fields, something London could not easily admit. From these considerations Tisza predicted that when superfluid helium flowed through a capillary, the normal component, being viscous, was blocked while the superfluid component could flow without dissipation. As a consequence, viscosity measurements could not give the same result when done in an open geometry, as in the torsion pendulum used in Toronto, or in a restricted geometry like a thin capillary or a porous medium. This resolved the apparent discrepancy between different experiments at that time. Moreover, he made an even more surprising prediction, namely that “a temperature gradient should arise during the flow of helium II through a thin capillary”. Tisza was right! Helium cools down when flowing so. Furthermore, he explained the fountain effect as the reverse of the previous one: the superfluid component moves in the direction of the hot side where the pressure increases. His 1938 article was only one column long. He developped his “two fluid model” in two subsequent articles Tisza2  where he made another fundamental prediction: thermal waves should exist where the superfluid and normal components move out of phase.

Nearly all the predictions by Tisza were proved correct by subsequent experiments. However, his model appeared slightly incorrect because the superfluid and normal components could not be defined in relation to a Bose gas. The two fluid model was subsequently developped by Landau in his 1941 article, with no reference, this time, to BEC. He considered energy states of the liquid instead of individual atoms. He attributed a mass density $`\rho _s`$ to the part of the fluid which occupies the ground state. The rest of the density $`\rho _n=\rho \rho _s`$ corresponded to the excited states. Furthermore, he had a remarkable intuition about the nature of these excited states. He assumed that single particle states were replaced by collective modes of two different kinds: phonons, i.e. sound quanta, and quantized vortices which he called “rotons”. Phonons had basically a linear dispersion relation $`\omega =ck`$ and rotons had a quadratic one with a gap ($`ϵ=\mathrm{\Delta }+p^2/2\mu `$), so that no excitation had a phase velocity $`ϵ/p`$ less than a certain minimum value $`v_c`$. From this assumption, Landau calculated the thermodynamics of helium II, the magnitude of thermal waves called “second sound”, etc. The major breakthrough in this article is the introduction of elementary excitations as quantized collective modes which led him to the crucial prediction of a critical velocity $`v_c`$. He explained that, due to the necessary conservation of energy and momentum, the superfluid component in motion at a velocity $`v`$ could not slow down by exciting collective modes if $`v`$ was less than $`v_c`$. In the absence of rotons, $`v_c`$ would have been the sound velocity c. Landau modified the dispersion relation of rotons in his 1947 article Landau47  in order to obtain a better agreement with thermodynamic properties. The phonon branch evolved continuously into a roton branch $`ϵ=\mathrm{\Delta }+(pp_0)^2/2\mu `$ around the finite momentum $`p_0`$. Given this phonon-roton spectrum, Landau’s critical velocity was predicted to be about 60 m/s at low pressure.

M. Cohen and R.P. Feynman Cohen  then proposed to verify Landau’s dispersion relation in neutron scattering experiments. This was done with great accuracy by several groups, in particular by Henshaw and Woods Henshaw . Of course, it was a strong support to Landau’s theory. In the mean time, superfluidity was proved not to exist in liquid helium 3 at similar temperatures, and this was a strong support to London and Tisza. As for the existence of Landau’s critical velocity, things appeared difficult. In fact, Landau had proposed this mechanism but he had not ruled out the possibility that other mechanisms could exist at lower velocity. The existence of quantized vortex lines was proposed in 1949 by Onsager Onsager  who predicted that the circulation of their velocity should be quantized as multiples of $`h/m`$. Independently in 1955, Feynman Feynman  explained that such vortices should lead to a much lower critical velocity of order $`v_c=(\mathrm{}/md)\mathrm{ln}(d/a)`$ where d is the size of the container (a vessel or a capillary diameter) and $`a`$ the atomic size of the vortex core. Feynman was right. The most beautiful evidence for the vortex mechanism was obtained by Avenel and Varoquaux Varoquaux  who could see quantized dissipation events in a flow through a small aperture (see Fig. 2). Moreover, their analysis of vortex nucleation shows a remarkable crossover from thermally activated nucleation above a crossover temperature of order 0.2 K to nucleation by quantum tunneling in the low temperature limit. A similar crossover was found in studies of superconducting Josephson junctions, also in our studies of cavitation in superfluid helium 4.

As for Landau’s critical velocity, it was also observed, but in cases where the system size was very small only. P. McClintock et al. measured the mobility of electrons moving in superfluid helium under special conditions. These electrons repell the electronic clouds of neighbouring atoms so as to form little empty cavities around them, whose typical radius is 10 Å. At low pressure and small electric field, these electron bubbles trap vortex rings as was checked from by their energy vs momentum curve. At high pressure (25 bar, near the melting pressure) and high electric field (from 20 to 2000 V/cm), the electron bubbles could in fact be accelerated to velocities slightly larger then Landau’s prediction. This was understood by Bowley and Sheard as a consequence of rotons being emitted in pairs McC . This whole subject would need a much longer review, but the comparison with gaseous superfluids is interesting. Superfluid gases have no rotons in their excitations. Landau’s critical velocity is thus associated with phonons. It is observed if the moving object is small (an atom), but if it is larger (a laser beam for exemple), then dissipation occurs at lower velocity. Just as in superfluid helium the nucleation of vortices is invoked to understand the low value of the critical velocity in this case. The comparison between quantized vortices in superfluid helium and in superfluid gases has been beautifully extended to the observation of similar vortex arrays in both cases (see J. Dalibard and C. Salomon, this conference).

## III The critical temperature and the condensate fraction in liquid helium

London had noticed that $`T_\lambda `$ was close to the transition temperature of an ideal Bose gas with the same density.

$$T_{BEC}=\left(\frac{2\pi \mathrm{}^2}{1.897mk_B}\right)n^{2/3}$$
(1)
Inserting a number density $`n=\mathrm{\hspace{0.25em}2.18}\times 10^{22}\mathrm{cm}^3`$ in Eq. 1 leads to $`T_{BEC}`$ = 3.1 K while $`T_\lambda `$ = 2.2 K. London had also noticed that the molar volume of liquid helium was large, due to the large kinetic energy corresponding to the quantum fluctuations of atoms in the cage formed by their neighbours. He had further noticed that the anomaly of the specific heat of liquid helium near $`T_\lambda `$ was not very different from the cusp one should expect in the case of the BEC of an ideal gas. He thus claimed that , for the understanding of the superfluid transition in helium, “it seems difficult not to imagine a connection with the condensation of the Bose-Einstein statistics”. However, in the next sentence, he immediately added that “a model which is so far away from reality that it simplifies liquid helium to an ideal gas” cannot lead to the right value of its specific heat.

In fact, the difficulty is much deeper than a question of specific heat only. As claimed by A.J. Leggett Leggett , it is not rigorously proven that a condensate, i.e. a macroscopic occupation of a ground state, should necessarily exist in a Bose liquid. It is not obvious either that there is a continuous path from the ideal Bose gas to the strongly interacting Bose liquid via the weakly interacting Bose gas where the existence of a condensation has been demonstrated both theoretically and , now, experimentally. The total absence of reference to BEC in Landau’s work is obviously intentional. Landau deliberately ignored London prize . He must have considered that there was no continuity from the gas to the liquid. It is now generally accepted that London was right, but, as we shall see, this issue is quite subtle.

Let us keep a few important steps only in the long history of this controversy. In 1947, Bogoliubov Bogoliubov  showed that, in a Bose gas with weak repulsive interactions, below the BEC transition, the low energy excitations are collective modes with a non-zero velocity. This is what allows such a weakly interacting Bose gas to be superfluid. Bogoliubov justified Landau’s assumption that there are no single particle excitations at low energy, so that a superfluid current cannot dissipate its kinetic energy if its velocity is less than a critical value, here the velocity of this collective mode, which is identical to the ordinary sound velocity.

In my opinion, the next crucial step is the historical paper by O. Penrose and L. Onsager Penrose  who generalized BEC to an interacting system. The single particle states are no longer eigen states. They considered the one-body density matrix $`\rho _1(r)`$ and its eigenvalues in the limit where $`r`$ tends to infinity. They correspond to the occupation of the various eigenstates of the whole liquid. If the occupation of the ground state is macroscopic, meaning that the “condensate fraction” $`n_0`$ is of order unity (when $`r`$ tends to infinity), then BEC has taken place. Since the momentum distribution is the Fourier transform of this one-body density matrix, a finite $`n_0`$ means a delta function at $`p`$ = 0 in the momentum distribution: there is a ground state with zero momentum which is occupied by a macroscopic fraction of the total number of particles. As explained by P. Sokol Sokol , the physical interpretation of $`\rho (r)`$ is the overlap of a wavefunction of the system when a particle is removed at the origin and replaced a distance $`r`$ away. The finite value of $`n_0`$ indicates a long range coherence which is often called “ODLRO” for “Off diagonal long range order”. The corresponding eigen function is the ground state wave function which is the order parameter for the transition. It is remarkable to see that, in the same article, Penrose and Onsager used Feynmann’s wave function for the condensate, approximated liquid helium as a gas of hard spheres, and ended with $`n_0`$ 8 % at low pressure, in agreement with more recent results. The ground state wave function has an amplitude and a phase $`\varphi `$. The superfluid velocity is

$$v_s=\frac{\mathrm{}}{m}\varphi $$
(2)
More theoretical efforts have been done in order to calculate $`n_0`$. In his review Sokol , Sokol cites several numerical methods from “Path integral Monte Carlo” (PIMC) to “Green’s function Monte Carlo” (GFMC) which are more or less consistent with the rough estimate first made by Penrose and Onsager. His GFMC predicts 9.2 % at low temperature and pressure. As for measuring the condensate fraction in liquid helium it has been another challenge. Contrary to what happens with trapped Bose gases, where the condensate is directly observed and easily measured, there is no experimental evidence for the existence of a condensate in liquid helium. The quantitative analysis of “Deep Inelastic Neutron scattering” experiments is very delicate. Reliable experimental values for $`n_0`$ were obtained in the last decade only. They depend not only on the assumption that a condensate exists, but also on the particular shape of the distribution function for excited states with non-zero momentum. Given this “caveat”, the presently accepted experimental value for $`n_0`$ in liquid helium at zero pressure and at low temperature is 10 $`\pm `$ 1.5 % , in agreement with theory.

It is now interesting to consider the pressure (or density) dependence of the condensate fraction. As shown by Sokol Sokol  again, theory and experiment agree about this also: $`n_0`$ decreases from about 10 % at $`P`$ = 0 bar where the density of liquid helium is $`\rho `$ = 0.145 $`\mathrm{g}.\mathrm{cm}^3`$ to about 5 % only at the melting pressure $`P_m`$ = 25.3 bar where $`\rho `$ = 0.172 $`\mathrm{g}.\mathrm{cm}^3`$ (Remember that liquid helium does not solidify at low temperature except if a pressure greater than 25.3 bar is applied). One understands the decrease of $`n_0`$ as a result of the exchange between atoms becoming more difficult as the density increases in this range. A probably equivalent explanation is given by Bauer et al. Bauer  who write that “exchange is inhibited by the increase in effective mass of the particles”. This variation is consistent with the transition temperature $`T_\lambda `$ being also a decreasing function of pressure. One now realizes that this behaviour is opposite to the prediction from Eq. 1. Does it mean that interactions always decrease the critical temperature for condensation? Certainly not, as shown by Gruter, Ceperley and Laloe Gruter . Gruter et al. calculated $`T_{BEC}`$ in a hard sphere fluid as a function of the parameter $`na^3`$ which describes the interaction strength ($`n`$ is the number density and $`a`$ is either a scattering length in the dilute case or a hard core in the dense limit of liquid helium). They found that the variation of $`T_{BEC}`$ is non-monotonic, a quite remarkable result. At low density, some interactions help the system to be homogeneous and to establish a long range coherence. On the contrary, large interactions make exchange difficult at large density. On Fig. 3, one sees that they were able to put the transition temperature of liquid helium on the same curve as their prediction for the critical temperature of BEC in dilute Bose gases. This curve would have surprised Landau !

I would like to make two further remarks on the pressure variation of $`T_\lambda `$ in liquid helium. Suppose now that, by some clever mean, one could study liquid helium at a density lower than its equilibrium value. Would $`T_\lambda `$ come closer to the ideal gas value $`T_{BEC}`$? We have shown that this is possible by using acoustic pulses of high intensity and short duration Balibar02 . In the absence of walls and impurities, it is possible to drive liquid helium to a metastable state at a negative pressure which approaches the “spinodal limit” at - 9.5 bar. Under such conditions, the density of liquid helium can be lowered to 0.10 $`\mathrm{g}.\mathrm{cm}^3`$, about 30 % less than in equilibrium at the saturated vapor pressure. Apenko Apenko  and Bauer Bauer  have calculated $`T_\lambda `$ in this metastable region of the phase diagram (see Fig. 4). They both predicted that $`T_\lambda `$ should reach a maximum value of about 2.2 K for a density of 0.12 $`\mathrm{g}.\mathrm{cm}^3`$ corresponding to -8 bar. In summary, in liquid helium at negative pressure, $`T_\lambda `$ ceases to be a decreasing function of pressure and comes closer to the ideal gas behaviour. This is consistent with our observation of a cusp in the cavitation pressure at 2.2 K Balibar02 . We naturally attributed this cusp to the crossing of the lambda transition at negative pressure, and we found it where it is predicted by recent theories, not on a linear extrapolation of the behavior at positive pressure. It is the sign that rotons are no longer dominating the thermodynamics of liquid helium near its spinodal limit.

We are also using the same acoustic technique to study overpressurized liquid helium. Our main motivation is to look for a possible limit of instability for the metastable liquid at high pressure with respect to the formation of the crystalline phase which is the stable one. In the course of this search, we have already reached about 120 bar without seing nucleation of helium crystals. This shows that it is possible to study liquid helium at much higher densities than previously thought. In his review , Sokol Sokol  considers this region as “unaccessible” but he was not aware of our experiments. It is thus not absurd to consider what happens to superfluidity in liquid helium at high density. A crude extrapolation of Sokol’s calculation of $`n_0`$ would predict that it becomes vanishingly small near a density of about 0.19 $`\mathrm{g}.\mathrm{cm}^3`$. Given the known equation of state Balibar02 , it would correspond to a pressure of order 55 bar, i.e. 30 bar above the equilibrium melting pressure. We have driven liquid helium up to a much higher density. It is hard to imagine that, if $`n_0`$ is very small $`T_\lambda `$ is not very small as well. How exactly does the lambda line extrapolate in this metastable region? Is it possible that a helium glass exists ? It is certainly not easy to study such speculative properties but we might obtain some information on it from our experiments. It would also be interesting to calculate the properties of a quantum hard sphere model in the very high density limit where the system jams, not only in the very dilute case where the system approaches the ideal gas.

With these few comments and ideas, I hope that I have provided some information for a comparison between superfluid liquids and superfluid gases. My purpose was also to show that a few questions deserve further study.


format: text

text:
# Spin Relaxation Caused by Thermal Excitations of High Frequency Modes of Cantilever Vibrations

## I Introduction

The method of oscillating cantilever-driven adiabatic reversals (OSCAR) is a promising method for spin detection in magnetic resonance force microscopy (MRFM) . Just like the ordinary MRFM technique , the magnetic moment of a sample driven by a resonant microwave magnetic field interacts with a ferromagnetic particle. If the ferromagnetic particle is attached to the cantilever tip, the magnetic moment of the sample changes the cantilever vibrations. This change can be detected through measurement of the cantilever vibration parameters. In the OSCAR technique the cantilever is driven by an external force with a feedback loop designed to keep the amplitude constant. The cantilever vibrations together with the resonant microwave magnetic field cause a cyclic adiabatic inversion of the magnetic moment of the sample. In turn, this cyclic inversion causes a shift of the cantilever vibration frequency, which is measured.

The main challenge in the OSCAR technique is the sharp decrease of the MRFM relaxation time when the distance between the ferromagnetic particle and the sample is below 1 micrometer . Thus, an important problem of the MRFM technique is to understand the nature of the spin relaxation in the submicron region and to find a way to reduce the relaxation rate. One of the possible mechanisms involved in MRFM relaxation in the submicron region is thermal currents in the metallic ferromagnetic particle . Another possible mechanism is thermal excitation of the high frequency modes of the cantilever vibrations .

In this paper, we explore this second mechanism. We show that the thermal excitations of the high frequency modes may provide a noticeable contribution to the experimentally observed spin relaxation of the MRFM signal. We also show that the use of a nonuniform cantilever can reduce the spin relaxation rate caused by this mechanism.

## II Model and equations of motion

The geometry of the problem is shown in Fig. 1. This geometry is essentially the same as that considered in . However, now we consider not just a single magnetic moment below the cantilever tip, but a uniform distribution of the magnetic moments across the sample which corresponds to the experimental conditions in . The ferromagnetic particle attached to the cantilever tip effectively interacts with all magnetic moments inside the resonant slice where the Larmor frequency matches the frequency of the microwave magnetic field, $`B_1`$. (The effective thickness of the resonant slice equals twice the amplitude of the cantilever vibrations.)

The tip vibrations associated with the fundamental mode of the cantilever can be described as a motion of an oscillator with a quality factor $`Q`$:

$$\ddot{z}_c+z_c+\dot{z}_c/Q=f(\tau ),$$
$`(1)`$
$$f(\tau )=\underset{k}{}\eta _k\mu _{kz},$$

$$\eta _k=\frac{\mu _0}{4\pi }\frac{3m\mu }{k_cz_0^5}\frac{\stackrel{~}{z}_k(5\stackrel{~}{z}_k^23\stackrel{~}{r}_k^2)}{\stackrel{~}{r}_k^7}.$$

Here we used the dimensionless quantities: $`\tau =\omega _ct`$, where $`\omega _c`$ is the “unperturbed” cantilever frequency of the fundamental mode; $`z_c`$ is the $`z`$-coordinate of the center of ferromagnetic particle in the units of the cantilever amplitude $`z_0`$; $`x_k`$, $`y_k`$, $`z_k`$ are the coordinates of the $`k`$-th magnetic moment of the sample in the same units; $`\stackrel{~}{z}_k=z_kz_c`$; $`\stackrel{~}{r}_k^2=x_k^2+y_k^2+\stackrel{~}{z}_k^2`$. $`\mu _0`$ is the permeability of free space; $`\stackrel{}{\mu }_k`$ is the $`k`$-th magnetic moment in units of its magnitude $`\mu `$ (which is the same for all magnetic moments); $`k_c`$ is the cantilever spring constant, $`m`$ is the magnetic moment of the ferromagnetic particle. The function $`f(\tau )`$ describes the action of the magnetic moments in the resonant slice on the cantilever.

The motion of the $`k`$-th macroscopic magnetic moment in the rotating reference frame is given by

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

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

In Eqs. (2) the following notation was used:

$$\mathrm{\Delta }_k=(\gamma B_{ext}\omega )/\omega _c+\frac{\mu _0}{4\pi }\frac{\gamma m}{\omega _cz_0^3}\frac{3\stackrel{~}{z}_k^2\stackrel{~}{r}_k^2}{\stackrel{~}{r}_k^5},$$
$`(3)`$
$$\epsilon =\gamma B_1/\omega _c.$$

Here $`\gamma `$ is the electron gyromagnetic ratio in the sample; $`B_{ext}`$ is the external static magnetic field; $`B_1`$ and $`\omega `$ are the amplitude and frequency of the rotating rf magnetic field. Note, that $`\mathrm{\Delta }_k`$ is the $`z`$-component of the rotating-frame effective field $`\stackrel{}{B}_{eff}^k`$ (in units $`\omega _c/\gamma `$), and $`\epsilon `$ is the $`x`$-component of the effective field in the same units, $`\stackrel{}{B}_{eff}^k=(\epsilon ,0,\mathrm{\Delta }_k)`$.

The upper and the lower boundaries of the resonant slice are determined by the equation

$$\mathrm{\Delta }_k=0,atz_c=\pm 1.$$

In our numerical experiments, the magnetic moments have been distributed uniformly inside the resonant slice. We used the following initial conditions

$$z_c=1,\mu _{kz}=1,$$
$`(4)`$
i.e. the magnetic moments are oriented approximately along the effective field in the rotating frame. To model the feedback technique in the OSCAR MRFM, our computer algorithm increased the value of $`z_c`$ to 1 every time the cantilever passed the upper point. The period of the cantilever oscillations was determined as the time interval between the instants of the $`z_c`$ maximum values.

To model the thermal excitations of the high frequency cantilever modes we use the following approach. The cantilever energy $`E`$ can be written as a sum of the energies of all vibrational modes. The high frequency modes can be described approximately by the formula

$$\omega _n=\frac{c_nt_c}{l^2}\left(\frac{Y}{12\rho }\right)^{1/2},$$
$`(5)`$
$$c_n=[\pi (n0.5)]^2,n>1,$$

where $`t_c`$ is the thickness of the cantilever, $`l`$ is its length, $`Y`$ is Young’s modulus, and $`\rho `$ is the density of the cantilever. Combining the expressions

$$\frac{m_c\omega _n^2a_n^2}{2}=k_BT,\omega _c^2=4k_c/m_c,$$
$`(6)`$
where $`m_c`$ and $`k_c`$ are the mass and the spring constant of the cantilever, we estimate the thermal amplitude of the $`n`$-th mode to be

$$a_n=\frac{1}{\mathrm{\Omega }_n}\left(\frac{k_BT}{2k_c}\right)^{1/2},\mathrm{\Omega }_n=\omega _n/\omega _c.$$
$`(7)`$
The amplitude of the cantilever tip oscillations, $`b_n`$, for any high frequency mode, $`n`$, is twice the amplitude of the mode :

$$b_n=2a_n=(2k_BT/k_c\mathrm{\Omega }_n^2)^{1/2}.$$

To describe the influence of the noise on the spin dynamics we replace the coordinate $`z_c`$ with $`z_c+\delta z_c`$ in the expression for $`\stackrel{~}{z}_k`$ and $`\stackrel{~}{r}_k`$ in (3), where

$$\delta z_c=\underset{n}{}\frac{b_n}{z_0}\mathrm{cos}(\mathrm{\Omega }_n\tau +\mathrm{\Psi }_n),$$
$`(8)`$
and $`\mathrm{\Psi }_n`$ being a random phase.

## III Results of numerical simulations

We solved the system of equations (1), (2) changing $`z_cz_c+\delta z_c`$ in Eqs. (2) according to (8). (We did not take into account the influence of the thermal noise on $`\eta _k`$ in Eqs. (1). As our simulations demonstrate, the influence of the thermal noise does not cause a significant direct contribution to the cantilever vibrations through the parameter $`\eta _k`$, but causes a dephasing of magnetic moments in the resonant slice.) The number of magnetic moments in the resonant slice was 200 ($`1k200`$). Our simulations show that the results do not change significantly when this number is increased to 400.

The parameters in our numerical simulations were taken from experiment :

$$B_{ext}=140\mathrm{mT},k_c=0.014\mathrm{N}/\mathrm{m},\omega _c/2\pi =21.4\mathrm{kHz},z_0=28\mathrm{nm},$$
$`(9)`$
$$Q=2\times 10^4,\mathrm{m}=1.5\times 10^{12}\mathrm{J}/\mathrm{T}.$$

The distance between the bottom of the ferromagnetic particle and the surface of the sample, $`d_1`$, (see Fig. 1) was $`d_1=700`$ nm, and the radius of the ferromagnetic particle was equal to $`d_1`$. The average magnetization of the sample was 0.89 A/m. The rotating rf field amplitude, $`B_1`$, in experiments was no less than 0.15 mT. The cantilever temperature, T, did not exceed 20 K. (For these values of parameters the value of $`d_2`$ in Fig. 1 was found to be $`d_2=875`$ nm, and the value of the magnetic field gradient at the center of the resonant slice was $`1.4\times 10^5`$ T/m.)

We have studied the decay of the OSCAR signal, $`\mathrm{\Delta }T/\mathrm{\Delta }T_0`$, where $`\mathrm{\Delta }T`$ is the shift of the period of the cantilever vibrations due to the influence of the magnetic moments of the sample, and $`\mathrm{\Delta }T_0`$ is the initial value of $`\mathrm{\Delta }T`$. In our model, the decay of the OSCAR signal is associated with cantilever noise near the Rabi frequency $`\omega _R=\gamma B_1`$ which is the smallest resonant frequency in the rotating reference frame.

High frequency cantilever modes may generate significant noise near the Rabi frequency causing a deviation of the magnetic moments from the direction of the effective field and a nonuniform dephasing between magnetic moments.

We write the amplitude of the cantilever tip thermal vibrations associated with the $`n`$-th mode as $`b_n=b_R\epsilon /\mathrm{\Omega }_n`$, where

$$b_R=\epsilon ^1(2k_BT/k_c)^{1/2}.$$

The phases $`\mathrm{\Psi }_n`$ in (8) were changed randomly between $`0`$ and $`2\pi `$ with random time intervals, $`\tau _\mathrm{\Psi }`$, between two successive changes of the phase. In particular, we put $`\tau _\mathrm{\Psi }=N\tau _0`$, where $`\tau _0`$ takes random values between $`2\pi /1.2\epsilon `$ and $`2\pi /0.8\epsilon `$, and $`N`$ is a free parameter of the model.

Our simulations show that the decay of the OSCAR signal is almost independent of $`N`$ for $`N<1000`$. (See Fig. 2.) We have found that for $`N<1000`$ the signal decay can be approximately described by an exponential function with relaxation time, $`\tau _m`$ (Fig. 2).

The relaxation time quickly decreases with an increase of the characteristic noise amplitude at the Rabi frequency, $`b_R`$ which, in turn, is proportional to $`\sqrt{T}/\epsilon `$. Taking the minimum value of $`\epsilon `$ ($`\epsilon =195`$, i.e. $`B_1=0.15`$mT) and the maximum value of $`T`$ ($`T=20`$K) from the range of parameters , we obtain the maximum value of $`b_R`$, $`b_R=1`$pm. The corresponding minimum value of the relaxation time is approximately $`\tau _m570`$ms. The experimental value found in , $`\tau _m=68`$ms, is smaller but of the same order of magnitude as our minimum theoretical value. Thus, we conclude that the thermal vibrations of the high frequency modes provide a noticeable contribution to the spin relaxation.

To study the spin relaxation we also considered values of parameters which provide relatively small relaxation time. This allowed us to reduce the simulation time and to determine the characteristic scaling properties of the relaxation process. As an example, Fig. 3 shows the decay of the OSCAR signal for various values of the cantilever amplitudes. One can see that the relaxation time $`\tau _m`$ decreases when the cantilever amplitude decreases. (See curves $`a`$ and $`c`$ in Fig. 3.) Note that if we take into consideration spins near the resonant slice, the value of $`\tau _m`$ slightly increases. (Compare curves $`b`$ and $`c`$ in Fig. 3.)

Fig. 4 shows the decay of the OSCAR signal for various values of the temperature and rotating magnetic field. One can observe the expected decrease of the relaxation time $`\tau _m`$ with an increase of the temperature or decrease of the rotating magnetic field.

Fig. 5 shows the decay of the OSCAR signal for various numbers of high frequency modes taken into consideration. One can see that the relaxation time $`\tau _m`$ slightly decreases when the number of high frequency cantilever modes increases from 2 to 22. A further increase of the number of modes does not essentially change the decay rate.

Our simulations show that the OSCAR signal generally consists of two parts: (a) a regular part which we studied in our paper and (b) a random part which could be observed after the decay of the regular signal. Numerical analysis of the random part of the OSCAR signal requires large computational times, and will be studied elsewhere.

## IV Reduction of the relaxation rate

In this section we discuss how to reduce the spin relaxation rate caused by thermal excitations of the high frequency cantilever modes. For this purpose we consider a cantilever with nonuniform thickness.

The mechanical energy of the cantilever can be written as

$$E=\frac{1}{2}_0^l\left[S\rho \left(\frac{z}{t}\right)^2+YI\left(\frac{^2z}{x^2}\right)^2\right]𝑑x.$$
$`(10)`$
Here $`S`$ is the cross-sectional area, $`I=wt_c^3/12`$, and $`w`$ is the width of the cantilever. The equation for the cantilever motion without an external force and damping is given by

$$S(x)\rho \frac{^2z}{t^2}=Y\frac{^2}{x^2}\left(I(x)\frac{^2z}{x^2}\right).$$
$`(11)`$
Using a new variable, $`\zeta =z\sqrt{S(x)}`$, we obtain the equation of motion in the form

$$\rho \frac{^2\zeta }{t^2}=\frac{Y}{\sqrt{S(x)}}\frac{^2}{x^2}\left(I(x)\frac{^2}{x^2}\frac{\zeta }{\sqrt{S(x)}}\right).$$
$`(12)`$
The general solution of the equation (12) can be written as

$$\zeta (x,t)=\underset{n}{}a_nf_n(x)\mathrm{exp}(i\omega _nt).$$
$`(13)`$
The eigenfunctions $`f_n(x)`$ are normalized in the following way:

$$_0^lS(x)f_n(x)f_m(x)𝑑x=V\delta _{nm},$$
$`(14)`$
where $`V=_0^lS(x)𝑑x`$ is the cantilever volume. The equation for the eigenfunctions is

$$\rho \omega _n^2f_n(x)=\frac{Y}{\sqrt{S(x)}}\frac{^2}{x^2}\left(I(x)\frac{^2}{x^2}\frac{f_n(x)}{\sqrt{S(x)}}\right).$$
$`(15)`$
We express the solution of Eq. (15) in the form

$$f_n(x)=\underset{k=1}{\overset{k_{max}}{}}c_{kn}\phi _k(x),$$
$`(16)`$
where $`\phi _k(x)`$ is the $`k`$-th eigenfunction for a cantilever with uniform cross-sectional area , satisfying the equation:

$$S\rho \omega _n^2\phi _n=YI\frac{^4\phi _n}{x^4},$$
$`(17)`$
$$_0^l\phi _n(x)\phi _m(x)𝑑x=l\delta _{nm}.$$

Taking $`k_{max}=70`$, we use a standard procedure to find the coefficients $`c_{kn}`$ and frequencies $`\omega _n`$. Finally, using the same approach as for the uniform cross-sectional area we can find thermal vibrations of the cantilever tip, $`\delta z_c`$, in (8) with new values of $`\mathrm{\Omega }_n`$ and $`b_n`$.

We found that it is possible to reduce the relaxation rate using a nonuniform thickness of the cantilever. (See also .) We considered a cantilever whose thickness increases near the location of the ferromagnetic particle:

$$\delta t_c(x)=\delta t_c(l)\mathrm{exp}[(xl)^2/a^2].$$
$`(18)`$
For the value $`a/l0.01`$, we obtained an almost tenfold increase in the relaxation time.

We also considered the case when the ferromagnetic particle was located away from the tip of the cantilever. We found that in this case the increase of the cantilever thickness near the location of the ferromagnetic particle causes a significant increase of the relaxation time. The details of these results will be published elsewhere.

## V Conclusion

We simulated the spin relaxation in OSCAR MRFM technique caused by the thermal excitation of the high frequency cantilever modes. Using a range of experimental parameters  we obtained the minimal relaxation time 570 ms which is greater than but of the same order of magnitude as the experimental value 68 ms. We demonstrated that an increase of the cantilever thickness at the location of the ferromagnetic particle can significantly increase the relaxation time.

Finally, we would like to mention two assumptions of our computational model: We assumed 1) that the mass of the ferromagnetic particle was small compared to the cantilever mass, and 2) that the ferromagnetic particle has a spherical shape. We are now simulating a modified cantilever-spin model which is not restricted by these two assumptions.

## Acknowledgments

We thank G.D. Doolen and D. Mozyrsky for discussions. 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:
# Non-Gaussianity of the derived maps from the first-year WMAP data

## 1 Introduction

With the first-year data release of the Wilkinson Microwave Anisotropy Probe (WMAP ), it has been proclaimed that we have entered the era of “precision cosmology”. The temperature fluctuations of the cosmic microwave background (CMB) radiation are believed to be the imprint of primordial density fluctuations in the early Universe which give rise to the large-scale structures we see today. Hence the data enable us to test the statistical character of the primordial fluctuations, making subsequent inferences on the topology and content of the Universe.

Although the WMAP team (Komatsu et al., 2003) claims that the signal is Gaussian with 95% confidence, the internal linear combination map released by the WMAP team is not up for CMB studies due to “complex noise properties” 111http://lambda.gsfc.nasa.gov/product/map/m\_products.html. Another group led by M. Tegmark has performed an independent foreground cleaning from the first-year WMAP data and made public their whole-sky CMB maps. Their foreground-cleaned map (hereafter FCM) and the Wiener-filtered map (WFM) are available online (Tegmark, 2003).

The FCM by the authors’ definition is such that the foreground contamination is removed as much as possible. As foregrounds are rather non-Gaussian, any residual after cleaning would manifest itself in the phase configuration. In this Letter we display the phases of the FCM and the WFM with color coding and implement our phase-mapping technique to test quantitatively the Gaussianity of both maps, based on the random-phase hypothesis of homogeneous and isotropic Gaussian random fields. Our phase mapping technique can play a crucial role as a qualitative criterion for component separation similar to the field of image reconstruction.

## 2 Gaussian Random Fields and the Random Phase Hypothesis

The statistical characterization of temperature fluctuation of CMB radiation on a sphere can be expressed as a sum over spherical harmonics:

$$\mathrm{\Delta }T(\theta ,\phi )=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_\mathrm{}mY_\mathrm{}m(\theta ,\phi ),$$
(1)
where $`a_\mathrm{}m=|a_\mathrm{}m|\mathrm{exp}(i\varphi _\mathrm{}m)`$. Homogeneous and isotropic Gaussian random fields (GRFs), as a result of the simplest inflation paradigm, possess Fourier modes whose real and imaginary parts are independently distributed. In other words, they have phases $`\varphi _\mathrm{}m`$ that are independently distributed and uniformly random on the interval $`[0,2\pi ]`$ (BBKS, 1986; Bond & Efstathiou, 1987). Thus the spatial variations should constitute a statistically homogeneous and isotropic GRF (BBKS, 1986) whose statistical properties are completely specified by its angular power spectrum $`C_{\mathrm{}}`$,

$$a_\mathrm{}ma_\mathrm{}^{^{}}m^{^{}}^{}=C_{\mathrm{}}\delta _{\mathrm{}\mathrm{}^{^{}}}\delta _{mm^{^{}}}.$$
(2)
The strict definition of a homogeneous and isotropic GRF requires that the amplitudes are Rayleigh distributed and the phases are random (Watts & Coles, 2003). At the same time, the Central Limit Theorem guarantees that a superposition of a large number of Fourier modes with random phases will be Gaussian. Therefore the random-phase hypothesis on its own serves as a definition of Gaussianity (BBKS, 1986).

## 3 Color-coded phase map of the derived WMAP maps

Tegmark, de Oliveira-Costa & Hamilton (2003) (TDH03) perform an independent foreground analysis from the WMAP data and provide a FCM and WFM. We first use a visual display of phases by colors to show phase associations (Coles & Chiang, 2000). In color image display devices, each pixel represents the intensity and color at that position in the image. Two color schemes are usually used for the quantitative specification of color, namely the Red-Green-Blue (RGB) and Hue-Saturation-Brightness (HSB) color schemes. Hue is the term used to distinguish between different basic colors (blue, yellow, red and so on). Saturation refers to the purity of the color, defined by how much white is mixed with it. Brightness indicates the overall intensity of the pixel on a grey scale. The HSB color model is particularly useful because of the properties of the ‘hue’ parameter, which is defined as a circular variable. Therefore we are mapping phases from 0 to $`2\pi `$ to the hue circle.

We have used the Healpix 222http://www.eso.org/science/healpix/ package to produce $`a_\mathrm{}m`$. In Fig. 3 we show the color-coded phase gradient $`D_{\mathrm{}}\varphi _{\mathrm{}+1,m}\varphi _{\mathrm{},m}`$ for the FCM and WFM. The vertical axis is the multipole $`\mathrm{}`$ up to $`\mathrm{}=600`$ and the horizontal the $`m`$ axis where $`m\mathrm{}`$. Due to the relation $`a_\mathrm{}m=a_{\mathrm{},m}^{}`$, only modes from non-negative $`m`$ are shown. Although the phase gradient (from neighboring modes) is the most primitive way of qualitatively checking phase correlations, the apparent presence of stripes shown in the FCM indicates strong coupling between modes of neighboring $`\mathrm{}`$ of the same $`m`$.

## 4 Phase mapping and the mean chi-square statistic of the derived maps

To test the Gaussianity of the FCM and the WFM based on the random-phase hypothesis, we apply a phase mapping technique (Chiang, Coles & Naselsky, 2002; Chiang, Naselsky & Coles, 2002) to quantify the degree of ‘randomness’ of the phases (i.e. Gaussian). The return map of phases is a bounded square in which all phase pairs of fixed separation $`(\mathrm{\Delta }m,\mathrm{\Delta }\mathrm{})`$ are mapped as points (see Fig. 4). For example, one single return map for phase pairs with separation $`(\mathrm{\Delta }m,\mathrm{\Delta }\mathrm{})=(0,1)`$ contains points with $`(x,y)`$ coordinate $`(\varphi _{\mathrm{},m},\varphi _{\mathrm{}+1,m})`$, i.e. all phase pairs from modes that are separated by $`\mathrm{\Delta }\mathrm{}=1`$. If the phases are random, we expect to have an ensemble of return maps of all possible separations, each of which should be a scatter plot. As such we are testing the ‘randomness’ on the most strict terms. After mapping phase pairs on to a return map, we can apply a mean $`\chi ^2`$ statistic on the return map, which is defined as

$$\overline{\chi ^2}=\frac{1}{M}\underset{i,j}{}\frac{\left[p(i,j)\overline{p}\right]^2}{\overline{p}},$$
(3)
where $`M`$ is the number of pixels on the return map, $`\overline{p}`$ is the mean value for each pixel on the discretized return map. Chiang, Naselsky & Coles (2002) have shown that for a homogeneous and isotropic GRF, return mapping of phases results in an ensemble of return maps, each with a Poisson distribution. The expectation value of the $`\overline{\chi ^2}`$ from such ensembles of Poisson-distributed maps is

$$\overline{\chi ^2}_\mathrm{P}=\frac{1}{4\pi R^2},$$
(4)
where $`R`$ is the scale of smoothing from a 2D Gaussian convolution in order to probe the spatial structure. The $`\overline{\chi ^2}_\mathrm{P}`$ will have a statistical spreading around $`\overline{\chi ^2}_\mathrm{P}`$ with a dispersion $`\mathrm{\Sigma }_\mathrm{P}`$ where

$$\mathrm{\Sigma }_\mathrm{P}^2=\frac{1}{\pi ^3R^2(M/2)}.$$
(5)
Figure 4 shows the histograms of the $`\overline{\chi ^2}`$ statistics from the ensemble of the return maps of the FCM and the WFM for 4 multipole bands. One of the advantages of phase mapping technique is that we are able to check Gaussianity in different multipole bands, in particular those corresponding to foreground contamination and noise. Here we present the $`\overline{\chi ^2}`$ statistic at 4 bands centered around $`\mathrm{}150`$, 290, 400 and 500: $`41<\mathrm{}<250`$ (roughly the first Doppler peak), $`224<\mathrm{}<350`$, $`350<\mathrm{}<450`$ and $`463<\mathrm{}<550`$. The solid dark and dotted gray curves are the WFM and FCM, respectively. In each panel the vertical line denotes the expectation value $`\overline{\chi ^2}_\mathrm{P}=(4\pi R^2)^1`$. The curves from the FCM are obviously skewed, and hence are manifestations of phase correlations (i.e. non-Gaussian).

In Fig. 4 we display the gross behavior of the distribution curves in terms of the arithmetic mean $`\overline{\chi ^2}`$ and the dispersion $`\mathrm{\Sigma }`$ from the mean chi-square statistic. The top panel is from the FCM and the bottom WFM. The contours mark 68% (solid curve) and 95% (dotted curve) CL regions from 2000 realizations of GRFs. The symbols corresponds to 4 multipole bands centered at $`\mathrm{}150`$, 290, 400 and 500. Note that the contour region in the bottom panel corresponds to a small section in the top panel. The phases of the 4 multipole bands from the FCM are all strongly correlated, so they are far away from the 95% CL region. The WFM, however, shows that phases below the first Doppler peak are random, with the other 3 multipole bands around the edge of 68% CL region.

We see evidence of non-Gaussianity, however, in the WFM of the following two bands centered $`\mathrm{}400`$ and 500. In the lower two panels of Fig. 4 there are points appearing at the tails above $`6\mathrm{\Sigma }_\mathrm{P}`$. On the other hand, among the 2000 realizations we simulate for GRFs, no mapping of phases reaches $`\overline{\chi ^2}`$ value over $`6\mathrm{\Sigma }_\mathrm{P}`$, setting the probability below 0.05% for a GRF to have such mapping. Phase mapping from the separation $`(\mathrm{\Delta }m,\mathrm{\Delta }\mathrm{})=(0,2)`$ produces $`\overline{\chi ^2}`$ value at $`7.3\mathrm{\Sigma }_\mathrm{P}`$ for the multipole band centered $`\mathrm{}400`$, also at $`6.5\mathrm{\Sigma }_\mathrm{P}`$ at $`(\mathrm{\Delta }m,\mathrm{\Delta }\mathrm{})=(1,2)`$. For the band $`\mathrm{}500`$, $`7.6\mathrm{\Sigma }_\mathrm{P}`$ appears at $`(\mathrm{\Delta }m,\mathrm{\Delta }\mathrm{})=(2,2)`$. These phase couplings are clear signs against the random-phase hypothesis, therefore a manifestation of non-Gaussianity.

We plot in Fig. 4 the CMB temperature map from only two multipoles $`\mathrm{}=350`$ and 352 (of all $`m`$’s) of the FCM and WFM. The choice of these specific multipoles of $`\mathrm{\Delta }\mathrm{}=2`$ from our previous calculation is to demonstrate non-Gaussian signals the correlated phases will produce in the map. The structures at $`\phi 0`$ and $`\pi `$ in the FCM, the residual signal after foreground cleaning, disappear after Wiener filtering.

## 5 Discussions

In this Letter we have tested non-Gaussianity of two maps: the Foreground-cleaned map and the Wiener-filtered map, which are processed by TDH03 from the WMAP data. Based on the random-phase hypothesis, we use a phase mapping technique to yield a statistic that has detected considerable non-Gaussian signals for both maps at most multipole bands. Our phase mapping technique is particularly useful in separating non-Gaussian contributions from different sources when various contaminations are present at different $`\mathrm{}`$ ranges. A multipole band which is considerably non-Gaussian could have an insignificant non-Gaussian contribution in the whole map and still produce an over-all Gaussian realization within a certain confidence level. As the uncertainties in foreground cleaning propagate through the data-processing pipelines to the accuracy of the angular power spectrum, it is therefore necessary to have effective methods in component separation. We believe that our phase mapping technique is a useful criterion to be incorporated into such methods. Our statistic based on phase mapping also holds great advantage when it comes to the issue of creating many whole-sky Gaussian realizations for Gaussian statistics. As our null hypothesis is that phases are random, we only need to put random phases (with Gaussian instrumental noise being automatically included) for each harmonic mode, which is easily done without any limit on the highest harmonic number $`\mathrm{}`$ from any pixelization scheme. It is worth mentioning that the upcoming Planck mission will have higher sensitivity and resolution, hence every step of data processing will be crucial in reaching such precision.

This paper was supported by Danmarks Grundforskningsfond through its support for the establishment of the Theoretical Astrophysics Center. We thank Max Tegmark et al. for providing their processed maps and making them public with openness. We thank Peter Coles and Max Tegmark for useful discussions. We also acknowledge the use of Healpix package (Górski, Hivon & Wandelt, 1999) to produce $`a_\mathrm{}m`$ and Fig. 4.


text:
# Guiding and reflecting light by boundary material

## Abstract

We study effects of finite height and surrounding material on photonic crystal slabs of one- and two-dimensional photonic crystals with a pseudo-spectral method and finite difference time domain simulation methods. The band gap is shown to be strongly modified by the boundary material. As an application we suggest reflection and guiding of light by patterning the material on top/below the slab.

, , ,

Photonic crystals are periodic dielectric structures that exhibit a photonic band gap, which can be utilized to control and confine light . Photonic crystal slabs are planar structures, i.e., they are periodic in a plane, but have finite thicknesses in the third dimension. The slab acts as a planar waveguide, if the average refractive index of the slab is higher than that of the surroundings. Thus light can be controlled three-dimensionally, when the slab is patterned to form a two-dimensional photonic crystal . Two-dimensional photonic crystal slabs are considerably easier to fabricate than three-dimensional photonic crystals. Analogously, light can be controlled two-dimensionally with one-dimensional photonic crystal slabs . The characteristics of photonic crystal slabs are notably different from the corresponding infinite photonic crystals and thus miniaturized photonic crystal components for, i.e., integrated optics cannot be designed based on the features of infinite structures.

Waveguides are basic components of an integrated optical circuit. Waveguides in two-dimensional photonic crystal slabs are formed by introducing line defects to the periodic structure . Light with a frequency in the band gap of the photonic crystal cannot propagate in the periodic structure and thus follows the line defect. We are proposing a new way to guide light using the influence that the material surrounding the photonic crystal slab has on the band gap. The photonic crystal slab has a band gap for different frequencies depending on the material that is on top/below the slab.

We consider photonic crystal slabs of heights $`h`$ of the order of the period $`P`$. We calculate the band structures of two-dimensional photonic crystal slabs using a pseudo-spectral method . Here, we show calculations for a specific geometry to illustrate the effects of different types of boundary materials above and below a photonic crystal slab. With boundary material we refer to the half-spaces above and below the slab. We show that varying the boundary material results in large changes in the band structures. We suggest to use this effect for guiding of light by patterning the boundary material — an alternative way for light guiding along crystal defects. The idea is the following (also explained in Fig. 2): planar two-dimensional photonic crystal would be fabricated with no defects and the waveguide would be realized by covering a channel (or leaving a channel uncovered) in otherwise uncovered (covered) slab. The difference of the band structures in covered/uncovered areas confines light in the waveguide. We demonstrate the idea with three-dimensional finite difference time domain (FDTD) simulations. For one-dimensional photonic crystal slabs the same effect can be used for reflecting light traveling along a one-dimensional photonic crystal slab by changing the boundary material abruptly. We demonstrate this with a two-dimensional FDTD method. The FDTD methods are described in Ref. .

To demonstrate that boundary material affects the band structure, we calculated band structures of two-dimensional photonic crystal slabs using a pseudo-spectral method . The band structures for a triangular lattice photonic crystal slab with boundary materials air and a dielectric, with a dielectric constant $`\epsilon _b=3`$, are shown in Fig. 1. It can be seen that when the boundary material is air, the frequencies around $`\omega P/(2\pi c)=0.3`$ fall into the band gap for even modes in the M-direction \[Fig. 1 (a)\]. However, when the boundary material has $`\epsilon _b=3`$, there is an even mode at $`\omega P/(2\pi c)=0.2937`$ at the M-point \[Fig. 1 (b)\]. Thus a waveguide can be made for this frequency in the M-direction by covering the photonic crystal slab in the waveguide region and leaving it otherwise uncovered. Light will follow the covered region where an allowed mode exists.

We performed three-dimensional FDTD calculations to study this type of waveguiding using the boundary material. The structures are shown in Fig. 2. We simulated the propagation of a Gaussian pulse with wavelength $`\lambda =1500`$ nm. The period of the photonic crystal lattice was $`P=450`$ nm, so that $`P/\lambda =0.3`$. The spectral full width at half maximum (FWHM) of the pulse was $`\mathrm{\Delta }\omega P/(2\pi c)=0.1`$, so that the main frequencies of the pulse fit into the band gap of Fig. 1 (a). First we calculated the transmission of the pulse along a defect waveguide made by removing holes from the photonic crystal slab \[Fig. 2 (a)\]. Then we calculated the transmission for a photonic crystal slab with no defects \[Fig. 2 (b)\] which was covered by strips of dielectric with $`\epsilon _{wg}=3`$ \[Fig. 2 (c)\]. The width of the strips was $`2P`$. The transmission of a Gaussian pulse for the guiding by boundary material \[Figs. 2 (b) and 2 (c)\] was 40% larger than in the traditional defect waveguide case \[Fig. 2 (a)\]. The energy densities of a Gaussian pulse propagating in the waveguide is presented in Fig. 3.

We studied also 60 and 90 turns in the hexagonal lattice. For the 60 turn (see Fig. 4), it was found that 51-55% of the light follows the waveguide turn. (Throughput of 51% is reached for FWHM $`\mathrm{\Delta }\omega P/(2\pi c)=0.1`$ and 55% for FWHM $`\mathrm{\Delta }\omega P/(2\pi c)=0.056`$. The width of the band gap is $`\mathrm{\Delta }\omega P/(2\pi c)=0.066`$.) In contrast, only a small fraction (about 10%) of the light followed the 90 turn because there is no solution in the K-direction for the used wavelength in the covered photonic crystal ($`P/\lambda =0.3`$). We repeated the simulations for a system where the photonic crystal was replaced by a material with a refractive index corresponding to the average index of the photonic crystal material (n = 2.9493): the amount of guided light was similar. Therefore, for this geometry and material, the effect of the photonic crystal reduces to an average refractive index material in a conventional ridge waveguide configuration. Finally, we repeated the simulations of the 60 and 90 turns for photonic crystal defect waveguides, and the guiding efficiency was considerably smaller (of the order of few percents). Note that the geometries and the coupling into the waveguides were not optimized, thus the results could probably be improved by optimizing the parameters. Furthermore, the comparison to ridge and defect waveguides refers to this specific example and should not be generalized without further study. However, the results demonstrate the in-principle usability of the idea.

The one-dimensional photonic crystal slab is a periodic stack of dielectric rectangular rods which have different dielectric constants. The considered geometry is shown schematically in Fig. 5. Parameters needed to define the structure are explained in the figure caption. Light traverses the structure ($`y`$-direction in Fig. 5), it reflects from each layer and interferes resulting in a band structure. We consider a polarization with field components $`H_x`$, $`E_y`$, and $`E_z`$. As discussed above, the band gap of a photonic crystal slab depends on the dielectric constant of the boundary material. In structures with periodicity in one dimension only, it can be utilized for reflecting light by changing the boundary material abruptly. For the considered structure, light with frequency around $`\omega P/(2\pi c)=0.3`$ has a solution in the photonic crystal when the slab is in air, but falls into the band gap when the slab is sandwiched between dielectric $`\epsilon _b=13`$. We can see from the two-dimensional FDTD simulations that light reflects at the point where the boundary material changes \[see Fig. 6(a)\]. Part of the light is diffracted to the boundary material that has a high dielectric constant. To estimate the efficiency of the reflection, we calculated the energy density of the light that was reflected back along the photonic crystal slab relative to the energy density of the incident pulse. This is illustrated in Fig. 6(b), where the energy density of the reflected pulse is shown as a function of the dielectric constant of the boundary material. It was calculated using the two-dimensional FDTD method. We can see that with increasing dielectric constant of the boundary material more light is reflected along the photonic crystal slab (indicating that for a totally reflecting boundary material such as metal the effect would be optimal). This and other properties of the crystal and its boundaries can be used to optimize the reflection. We have also studied a slab with a refractive index that is the average of that of the photonic crystal slab ($`n=\sqrt{130.2+10.8}`$). The considered boundary materials were dielectric with $`\epsilon _b=13`$ and metal (Au) with $`n=0.053+7.249i`$. In these cases only a small fraction of the light is reflected back along the waveguide. This shows that, for this geometry and choice of materials, the effect of the photonic crystal cannot be reduced to a homogeneous material with a corresponding average refractive index.

Our studies show that in some cases (the 2D periodic structure considered here) the guiding of light by patterning the surrounding material of a photonic crystal slab seems to reduce to the use of a conventional ridge waveguide. However, in other cases (the 1D example), the effect of the photonic crystal clearly cannot be reduced to a homogeneous material with the average refractive index. The situation somewhat resembles the case of microstructured fibers: there, in some cases the guiding can be explained by a conventional step index fiber structure with an effective refractive index, whereas in other cases the guiding is a genuine photonic band gap effect. In both cases, the fibers have turned out to provide interesting applications.

The advantages of the guiding scheme proposed here depend on the type of guiding it actually provides. In all cases (even when the guiding is effectively equivalent to guiding by a homogeneous material ridge waveguide), the advantage over using a ridge waveguide is the integrability with photonic crystal components, e.g. active elements (one may even imagine realizing microcavities by the effect) and additional degrees of freedom in designing the mode properties. The advantage over defect waveguides may be better guiding in some cases (as demonstrated here), easier fabrication or better integrability. In addition to these advantages, when the guiding is a genuine photonic band gap effect, one should be able to realize, for suitable crystal symmetries, sharp bends that are not possible with ridge waveguides. For the hexagonal lattice considered here, it was unfortunately not possible to demonstrate this since the frequencies of the band gaps as well as the frequency of the guided mode differ considerably in the K- and M-directions.

In summary, we have studied the effect of boundary materials and finite height of one- and two-dimensional photonic crystals to the band gap of the photonic crystal. We used a pseudo-spectral method and a finite difference time domain method. The band gap shows dependency on the boundary material. The strong effect of the boundary material on the band structures can be used for novel type of guiding of light in two-dimensional photonic crystals: the channel would be realized by covering a part (having the shape of the channel) in otherwise uncovered slab. We tested the idea by FDTD simulations of pulse propagation in a two-dimensional photonic crystal, and the proposed guiding by boundary material provided twice as effective transmission as the traditional guiding by a defect, for the (unoptimized) parameters used. Whether or not the guiding can be understood by replacing the photonic crystal by a homogenous material with the average refractive index was shown to depend on the choice of geometry and materials.

We thank the Academy of Finland for support (Project Nos. 48845, 53903, 205454). Jerome Moloney and Aramais Zakharian (Arizona Center for Mathematical Sciences, University of Arizona) are acknowledged for providing the original 3D FDTD code and Juuso Olkkonen (VTT Electronics, Finland and Optical Sciences Center, University of Arizona) for development and support of the 3D FDTD code.


text:
# Tunneling into Multiwalled Carbon Nanotubes: Coulomb Blockade and Fano Resonance

Present address: Gordon McKay Laboratory of Applied Science, Harvard University, MA 02138, USA\]

## Abstract

Tunneling spectroscopy measurements of single tunnel junctions formed between multiwalled carbon nanotubes (MWNTs) and a normal metal are reported. Intrinsic Coulomb interactions in the MWNTs give rise to a strong zero-bias suppression of a tunneling density of states (TDOS) that can be fitted numerically to the environmental quantum-fluctuation (EQF) theory. An asymmetric conductance anomaly near zero bias is found at low temperatures and interpreted as Fano resonance in the strong tunneling regime.

Coulomb blockade (CB) has been studied intensely in a multi-junction configuration, in which electron tunnel rates from the environment to a capacitively isolated “island” are blocked by the e-e interaction if the thermal fluctuation is below the charging energy $`E_c=e^2/2C`$ and the quantum fluctuation is suppressed with sufficiently large tunnel resistance $`R_tR_Q=h^2/2e`$. In the case of a single-junction circuit, the understanding of CB is less straightforward. The Coulomb gap, supposedly less significant for the case of low-impedance environments, should be established only if the environmental impendence exceeds $`R_Q`$.

In single-walled carbon nanotubes (SWNTs), the reduced geometry gives rise to strong e-e interaction. Indeed, CB oscillations and evidence of Luttinger liquid (LL) have been observed Bockrath . In contrast to SWNTs, in which only two conductance channels are available for current transport, MWNTs with diameter in the range of $`d`$=20-40 nm have several tens of conductance channels. The energy separation of the quantized subbands, given by $`\mathrm{\Delta }E=\mathrm{}\upsilon _f/d`$, is about 13-26 meV taking the Fermi velocity $`\upsilon _f=8\times 10^5`$ m/s. This value is about an order of magnitude smaller than that of SWNTs. Experiments indicate that MWNTs are considerably hole-doped, thus a large number of subbands, on the order of ten, are occupied. Observations of weak localization Langer , electron phase interference effect Bachtold2  and universal conductance fluctuations Langer  support the view that low frequency conductance in MWNTs is contributed mostly by the outmost graphene shell and is characterized by 2D diffusive transport. In addition to the phase interference effects, a strong e-e interaction has also been observed in MWNTs. Pronounced zero-bias suppression of the TDOS has been observed several times in the tunneling measurements Haruyama ; Tarkiainen ; Bachtold . Moreover, the TDOS shows a power-law, i.e., $`\nu (E)E^\alpha `$, which resembles the case of a LL. It is noteworthy that in the EQF theory, for a single tunnel junction coupled to high-impedance transmission lines, such a scaling behavior is also predicted at the limit of many parallel transmission modes Ingold . The physical origin of these power laws is the linear dispersion of bosonic excitations that are characteristic both for LL, which is a strictly 1D ballistic conductor, and a single tunnel junction connected to a 3D disordered conductor. In the latter case, the quasiparticle tunneling is suppressed at $`Ve/2C`$, therefore the charge is transported with 1D plasmon modes. The fact that MWNTs have many conductance modes, together with the observation of a crossover from power-law to Ohmic behavior at higher voltages Tarkiainen , suggests that the EQF theory is more appropriate to describe the observed TDOS renormalization. However, most of these measurements were done in multi-junction configurations. A single tunnel junction measurement is needed to further clarify this issue.

In this Letter, millimeter-long CVD-grown MWNTs Pan  are measured by a cross-junction method. The MWNTs samples are composed of loosely entangled nanotubes of diameters between 20 to 40 nm that are roughly parallel to each other and up to 2 mm long. The single tunnel junctions are formed by crossing a very thin ($`<1\mu `$m) MWNT bundle with a narrow strip of metal wire fabricated on an insulating substrate. We explored different electrode materials including Au, Cu, Sn, and Al. In the case of Sn and Al, a small magnetic field is applied to suppress the superconducting state below $`T_c`$. No obvious change of device characteristics attributable to the choice of metals is observed. With such a cross-junction configuration the measured conductance is contributed exclusively by the tunnel junction, since the current is passed along one arm of the MWNT/metal and the voltage is measured along the other, non-current carrying arm. Thus the device can be understood as a small number of single junctions in parallel. Despite the simplicity of their fabrication, we find that the devices are very stable and sustain several cooling cycles without apparent change of characteristics. More than 20 samples are measured, all yielding strong zero-bias suppression of the TDOS.

The current $`I(V)`$ and $`GdI/dV`$ are calculated by a golden-rule approach incorporating the environmental influence by $`P(E)`$, the probability for a tunneling electron to lose energy to the environment. $`P(E)`$ can be calculated from its Fourier transform, the phase correlation function $`J(t)`$. In the transmission-line model, $`J(t)`$ is determined by the total environmental impedance $`Z_t(\omega )=(i\omega C+Z^1(\omega ))^1`$, where $`Z(\omega )`$ is the external environmental impedance. We apply the method in Ref. Ingold2  to evaluate $`P(E)`$ from an integral equation without the need of going to the time domain. The I-V characteristics are then calculated for finite temperatures and arbitrary junction impedance. The parameters that appear in the numerical calculation are the damping strength $`\alpha =Z(0)/R_Q`$, the inverse relaxation time $`\omega _{RC}`$, and the quality factor $`Q=\omega _{RC}/\omega _S`$ with the resonance frequency of the undamped circuit $`\omega _S=(LC)^{1/2}`$. It is noted that the quality factor $`Q`$ plays only a minor role as an additional adjustable parameter. Therefore it is not important to the fit, and is always set to unit.

We then try to fit our experimental data with the EQF theory. Since $`E_c`$ and $`R_t`$ are determined by the high-voltage data, the isolation resistance $`R_{iso}=Z(0)`$ is the only adjustable parameter in our fitting. In contrast to the case of Ref. Zheng , where $`R_{iso}`$ is formed by ideal Ohmic and temperature-independent resistors, $`R_{iso}`$ in our case are provided by the resistive impedance of the MWNTs themselves, which should be temperature-dependent. Indeed, from the fitting, we see an evolution of $`\alpha `$ from 0.1 at 20 K to 0.25 at 1-4 K and then $`\alpha `$ saturates (Fig. 1c). Note that dimensionless units are used with $`G`$ normalized by 1/$`R_t`$ and voltage normalized by $`e/2C`$. Therefore the number of the single tunnel junctions has no effect on the fitting, and our measurements of different samples can be directly compared. We find that for different samples with scattered characteristics (see Table I), at low temperatures the exponent $`\alpha `$ reaches a universal value of 0.25-0.35, which agrees with previous multi-junction measurements Tarkiainen ; Bachtold . It yields $`R_{iso}`$ = 3.3-4.6 k$`\mathrm{\Omega }`$, which is roughly a constant for different samples. This coincidence is not accidental but reflects the intrinsic electrodynamic modes of the MWNTs, which can be modelled as an ideal resistive LC transmission line: $`R=(L^{}/C^{})^{1/2}`$ with the kinetic inductance estimated as $`L^{}=R_Q/2N\upsilon _f1`$ nH/$`\mu `$m for $`N`$ 10 - 20 modes and the capacitance $`C^{}`$ 20-30 aF/$`\mu `$m. Therefore, the “environment” with respect to the single-junctions in our devices is provided by MWNTs themselves, not by the external circuits.

As mentioned previously, for an isolated single tunnel junction with many parallel transport modes, the EQF theory predicts a power-law asymptotics: $`G(V,T)/G(0,T)f(V/T)=|\mathrm{\Gamma }[\frac{1}{2}(\alpha +2)+\frac{ieV}{2\pi k_BT}]/\mathrm{\Gamma }[\frac{1}{2}(\alpha +2)]\mathrm{\Gamma }[1+\frac{ieV}{2\pi k_BT}]|^2`$, where $`G(0,T)`$ is the zero-bias conductance, and $`\mathrm{\Gamma }(x)`$ is the Gamma function Zheng . For $`eV/2\pi k_BT1`$, a voltage power law $`G(V,T)V^\alpha `$ is expected. Note that in the above scaling function, exponent $`\alpha `$ is the only adjustable parameter. In Fig. 2, we can observe such scaling behavior. Here $`f(V/T)`$, represented by solid line, is calculated taking $`\alpha =0.27`$. The inset shows the original data in which the dashed line is calculated numerically using the same exponent $`\alpha =0.27`$.

CB in MWNTs is recently described with a microscopic theory in Ref. Egger . In the 2D diffusive regime, the exponent for tunneling into the bulk of a MWNT is: $`\alpha =(R/h\nu _0D)ln(1+\nu _0U_0)`$. Here $`R`$ is the tube radius, $`U_0`$ is the intra-tube Coulomb interaction, $`D=\upsilon _f^2\tau /2`$ is the charge diffusivity and the “bare” DOS $`\nu _0=N/4h\upsilon _f`$ with $`N20`$. If we use $`\upsilon _f\tau 60`$ nm from the magnetoresistance data KangMR , then under the condition of $`U_0/h\upsilon _f1`$ the above formula yields $`\alpha 0.25`$, which agrees with our experiment.

Besides the voltage power law phenomenon, which is characteristic of strong Coulomb interactions, another energy scale - the discrete energy levels due to electron confinement - emerges as temperature decreases, and it modifies the tunneling spectra. Unlike the case of a two-junction configuration, where the electrons are confined by the two contacts if the nanotube is clean enough, in a single-junction configuration the electrons can be confined by disorders when the MWNTs are “dirty”, such that the impedance of the local environment of the junction is larger than $`R_Q`$. A stacking mismatch between adjacent walls and other structural imperfections are possible sources of disorders in MWNTs, resulting in discrete energy levels.

By cooling down the devices to below 1 K, we observe that $`G`$ develops a narrow resonance-like anomaly at very low bias (Fig. 3). The asymmetric anomaly builds up consistently as temperature decreases and even shows a dip structure. The line shape resembles that of a Fano resonance.

Unlike the Coulomb structure, which is caused by static e-e interaction, a Fano resonance arises from an e-e exchange interaction between two interfering scattering channels: a discrete energy level and a continuum band. It has recently been rediscovered in mesoscopic systems such as semiconductor quantum dots, SWNTs, etc. Gores ; Wiel ; Nygard ; Liang . Taking a quantum dot as an example, it can be considered as a gate-confined droplet of electrons with localized states. The coupling of the dot to the leads can be tuned to control the system to enter different transport regimes: If the dot is weakly coupled with the environment, a well-established CB develops. The charge transport is suppressed except for narrow resonances at charge degeneracy points. When the tunnel barriers become more transparent, the dot enters the Kondo regime – i.e., below a characteristic Kondo temperature $`T_K`$, spin-flip co-tunneling events introduce a narrow symmetric TDOS peak at $`E_F`$ that can be interpreted as a discrete level. If the coupling is strong enough, the interference between this discrete level and the conduction continuum gives rise to an asymmetric Fano resonance. In Ref. Gores , such a crossover from a well-established CB through Kondo regime to Fano regime has been clearly observed.

We find that the asymmetric resonance curve of $`G`$ can be fitted by the Fano’s formula: $`G(ϵ+q)^2/(ϵ^2+1)`$. Here $`q`$ is the so-called asymmetry parameter, $`ϵ=(eV\epsilon _0)/(\gamma /2)`$ is the dimensionless detuning from resonance, and $`\gamma `$ is the FWHM of the resonance. As expected, the asymmetric parameter $`q`$ as a measure of the degree of coupling between the discrete state and the continuum increases when the temperature drops.

The Kondo temperature, estimated from the relation $`\gamma =2k_BT_K`$, is consistent for each device at different temperatures, and agrees with the observed FWHM of the resonance. If Kondo physics truly exists, then $`G`$ at voltages near the resonance should show non-monotonic temperature dependence around $`T_K`$ Nygard . Indeed, the $`GV`$ curves in Fig. 3 exhibit such behavior: $`G`$ at the peak position first drops with $`T`$ and then rises up below $`1`$ K (plotted in Fig. 4b). Moreover, we find that a perpendicular magnetic field gives rise to effects that are two-fold (Fig. 4a and 4c): First, the background conductance increases monotonically with the applied $`B`$. Second, the dip seen at zero field gradually disappears and the resonance turns into a nearly symmetric peak at high field, similar to the pattern observed in semiconductor quantum dots Gores .

The magnetic field effect can be explained as follows: First, adding a flux should change the amplitude and/or the phase for the resonant channels and therefore break down the coherent backscattering and increase the forward transmission through the channel. Second, the magnetic field can destroy the interference between the resonant and nonresonant paths, transforming a resonant dip into a peak.

Since Kondo physics is historically interpreted as the interplay between the d-orbitals of magnetic impurities and the conduction continuum, we have to eliminate the possibility that the observed Fano resonance comes from residual traces of the Fe/Si catalyst in the MWNT samples. It has been shown that magnetic impurities in MWNTs cause an enhancement of thermoelectric power Grigorian , which is absent in our careful TEP measurements KangTEP . Furthermore, no Fe signature can be detected in the body of the MWNT bundles within the instrumental resolution in the energy dispersion X-ray and TEM studies. Therefore the observed Fano resonance must be an inherent property of MWNTs, which can be virtually treated as quantum dots strongly coupled to the leads. We noticed that recently Buitelaar et al. have also seen clear traces of CB and Kondo resonance in MWNT SET devices Buitelaar . Their observed Kondo temperature $`T_K=1.2`$ K is in good agreement with our results.

In summary, strong zero-bias suppression of the TDOS is observed in MWNTs’ single tunnel junctions that can be explained well by the EQF theory. The observed exponent $`\alpha `$ 0.25–0.35 is found to be consistent for all the samples. Similar to the case of an open quantum dot, in low-impedance junctions we find that a Fano-resonance-like asymmetric conductance anomaly builts up below mV energy scales. It seems that MWNTs provide us a good laboratory to study the interplay between strong e-e interaction and disorder scattering.

The authors acknowledge fruitful discussions with R. Egger, M. Bockrath, W. Zheng, T. Xiang and G.M. Zhang. This work is supported by the National Key Project for Basic Research, the Knowledge Innovation Program of CAS, and NSFC.


text:
# A Test of a Test for Chaos

## I Introduction

A promising new binary test for chaos has been proposed by Gottwald and Melbourne . The test returns 0 for a non-chaotic system and 1 for a chaotic system. It does not offer any quantitative information and so cannot quantify the degree of irregularity in the way a Lyapunov exponent might. However, it can provide a quick and useful diagnostic for chaos.

Consider a base dynamical system $`\dot{𝐱}=𝐟(𝐱)`$ of any dimension with orbits $`𝐱(t)`$. The proposed method to determine if it is chaotic is based on a Euclidean $`𝐄(2)`$ group extension of the underlying base dynamics . The dynamics is explicitly extended to include two new variables $`(p,q)`$ defined through

$`\dot{p}`$ $`=`$ $`\varphi (𝐱)\mathrm{cos}(\omega _0t),`$ (1)
$`\dot{q}`$ $`=`$ $`\varphi (𝐱)\mathrm{sin}(\omega _0t),`$ (2)
where $`\varphi (𝐱)`$ is any observable of the base dynamics and $`\omega _00`$ is an arbitrary constant frequency that is needed to damp off any linear growth that may be common to both non-chaotic and chaotic orbits. If the observable $`\varphi `$ is drawn from a non-chaotic system, then the $`(p,q)`$ subspace will be bounded. But if the observable $`\varphi `$ is drawn from a chaotic system, then the motion in the $`(p,q)`$ subspace will be Brownian and unbounded. The extended variables diffuse through the subspace as the observables of the base dynamics jump around unpredictably.

Gottwald and Melbourne introduced a binary test for chaos by defining a mean-square displacement

$$M(t)=\underset{T\mathrm{}}{lim}\frac{1}{T}_0^T\left(p(t+\tau )p(\tau )\right)^2𝑑\tau $$
(3)
and characterizing the behaviour of $`M(t)`$ through

$$K=\underset{t\mathrm{}}{lim}\frac{\mathrm{log}M(t)}{\mathrm{log}t}.$$
(4)
If there is no chaos then the motion in $`p(t)`$ is bounded and $`K0`$. If the base dynamics is chaotic then $`p(t)`$ will exhibit Brownian diffusion so that $`\mathrm{\Delta }pt^{1/2}`$ as $`t\mathrm{}`$ and $`K1`$. The $`K`$-test can be used for both continuous and discrete dynamical systems even when the precise underlying base dynamics is unknown .

This test has certain advantages of simplicity and was shown to confirm the results of testing for chaos using Lyapunov exponents for the forced van-der-Pol system . We show here that the $`K`$-test also confirms the trend of the Lyapunov exponents for the simple Hamiltonian system of Henon and Heiles . So the $`K`$-test is applicable to Hamiltonian as well as dissipative systems.

While the test works well as a diagnostic of the transition from regularity to chaos we express some reservations about the use of the test on an orbit by orbit basis. Most importantly, we found the method difficult to interpret for systems that experienced chaotic transients in contrast to the easy interpretation of dissipative systems near a strange attractor, or in highly chaotic Hamiltonian systems. To illustrate this we use the Lorenz model  to show that transient chaos, which can be found in the Lorenz system by other means, goes undetected by the $`K`$-test. Orbits that enter highly chaotic regions of phase space and then depart into regular regions are not Anosov and may not show the diffusive motion in phase space required to yield a $`K1`$ result from the test. These transient systems may not satisfy the criteria for a system within the remit of the $`K`$-test. This signals one of the possible limitations of the K-test.

## II Henon-Heiles

The Henon-Heiles system is an ideal Hamiltonian system on which to evaluate the utility of the $`K`$-test. It provides a simple and well understood model for the motion of stars in a galactic potential as well as for the motion of non-linearly-coupled molecules. The transition from regularity to chaos as the energy is increased is well documented and can be identified with Poincaré surfaces of section, Lyapunov exponents, and the Painlevé property.

The Henon-Heiles Hamiltonian is

$$H=\frac{1}{2}\left(p_x^2+p_y^2+x^2+y^2\right)+x^2y\frac{1}{3}y^3$$
(5)
with equations of motion

$`\dot{x}`$ $`=`$ $`p_x`$ (6)
$`\dot{y}`$ $`=`$ $`p_y`$ (7)
$`\dot{p}_x`$ $`=`$ $`x2xy`$ (8)
$`\dot{p}_y`$ $`=`$ $`yx^2+y^2.`$ (9)
We add the Euclidean extension

$`\dot{p}`$ $`=`$ $`x\mathrm{cos}(\theta )`$ (10)
$`\dot{q}`$ $`=`$ $`x\mathrm{sin}(\theta )`$ (11)
$`\dot{\theta }`$ $`=`$ $`\omega _0`$ (12)
where we have chosen $`\varphi =x`$. We evolve these equations numerically and check that the Hamiltonian remains conserved throughout the simulation. Our results are summarized in Figs. 1 and 2. The orbits selected for display have identical initial conditions up to the energy $`H`$ which was varied from orbit to orbit. The value of $`K(t)`$ is measured as the slope of the line $`ln(M(t))`$ versus $`ln(t)`$ and is plotted in Fig. 1 as a function of the energy $`H`$. The transition from non-chaotic orbits ($`K0`$) to chaotic orbits ($`K1`$) as the energy increases is clearly demonstrated. The transition confirms the break up of tori in the phase space for the base dynamics (and therefore the non-integrability of the base dynamics) as the energy grows. We also confirm the onset of chaotic motion for these orbits with a numerical determination of the principal Lyapunov exponent shown as a function of $`H`$ in Fig. 2.

A few comments on Fig. 1 should be made. The convergence of the value of $`K`$ towards either $`0`$ or $`1`$ can be improved if the numerical simulations are allowed to run far longer and a more precise extrapolation method is used so as to avoid any early transients. Of course, this renders the application more numerically intensive and slower without yielding a very great improvement in convergence. There is also a clear dependence of the motion in the Euclidean subspace on the specific value of $`\omega _0`$ employed. However, the overall trend in $`K`$ is typically the same for various values of $`\omega _0`$.

While Fig. 1 does show the required transition from 0 to 1, it is important to stress how the conservative system differs from the dissipative system. For a dissipative system all trajectories, regardless of initial conditions, are drawn onto the same attractor. Therefore the $`K`$-test should show no dependence on initial conditions. This is not true for a Hamiltonian system. There can be a mixture of regular and irregular orbits and so the test result can depend on initial conditions. Only for a completely ergodic system will $`K1`$ independent of initial conditions. Also, since chaotic transients can often arise in conservative systems, the $`K`$-test may not always give such a crisp transition. Chaotic transients in dissipative systems are discussed further in the next section.

## III Lorenz

The standard Lorenz system demonstrates a transition from no chaos to transient chaos and then further to a full chaotic attractor. The standard system is

$`\dot{X}`$ $`=`$ $`\sigma (XY),`$ (13)
$`\dot{Y}`$ $`=`$ $`XZ+rXY,`$ (14)
$`\dot{Z}`$ $`=`$ $`XYbZ,`$ (15)
with $`\sigma =10`$ and $`b=8/3`$ and we vary the constant $`r`$. The transition to chaos is summarized in Ref. . For $`r<1`$ there is one fixed point at $`X=Y=Z=0`$ and no chaos. For $`1<r\mathrm{}<13.96`$ there are two additional attractors, one at

$`X_R`$ $`=`$ $`(b(r1))^{1/2}`$ (16)
$`Y_R`$ $`=`$ $`(b(r1))^{1/2}`$ (17)
$`Z_R`$ $`=`$ $`r1,`$ (18)
and another at $`(X_L,Y_L,Z_L)=(X_R,Y_R,Z_R),`$ but still no chaos. A non-chaotic orbit drawn onto the simple fixed point $`(X_R,Y_R,Z_R)`$ with $`(\dot{X},\dot{Y},\dot{Z})0`$ is shown in Fig. 3. A closely neighbouring orbit will follow a similar path onto the same attractor since there is no exponential sensitivity to initial conditions, no positive Lyapunov exponent, and no chaos.

However, for values of $`r\mathrm{}>13.96`$ a non-attracting form of transient chaos develops around a homoclinic orbit in phase space, first noted by Kaplan and Yorke . Two orbits undergoing transient chaos are shown in Fig. 4. One orbit winds around a bit before being drawn onto the regular attracting point $`(X_L,Y_L,Z_L)`$ while its very near neighbour ends up on $`(X_R,Y_R,Z_R)`$.

We choose to use the method of fractal basin boundaries to locate regions of transient chaos. An initial slice in phase space is color-coded according to whether the orbit eventually lands on the fixed point $`(X_R,Y_R,Z_R)`$ (black) or the fixed point $`(X_L,Y_L,Z_L)`$ (white). In the absence of chaos all the basins will have smooth boundaries. In the presence of chaos the boundaries become fractal, demonstrating both extreme sensitivity to initial conditions and a chaotic mixing of orbits. Such fractal basin boundaries are illustrated in Fig. 5; the orbits shown in Fig. 4 are drawn from along the fractal. They have nearly identical initial values but divergent outcomes, a characteristic of chaotic transients.

For $`r\mathrm{}>24.06`$ the transient chaos gives way to chaos on attractors which merge into the famous Lorenz strange attractor beyond $`r\mathrm{}>24.74`$ . An orbit that drifts onto the strange Lorenz attractor is shown in Fig. 6.

We apply the $`K`$-test to the Lorenz system and show that it effectively marks the transition from non-chaotic motion to chaos on a strange attractor at $`r24.74`$ where we see $`K`$ rising from $`0`$ to $`1`$. However, the test is unable to pick up the chaotic transient behaviour for values of $`13.96\mathrm{}>r\mathrm{}>24.06`$. Notice that a similar random scan through the Lyapunov exponents in Fig. 8 also misses the transient episodes although short time exponents along the fractal basin boundary can be isolated.

Specific orbits of Fig. 7 can be isolated to illustrate the behaviour in the $`(p,q)`$ subspace explicitly. The orbit at $`r=18`$ corresponds to the right-most orbit in Fig. 4. The orbit winds around before drifting onto the fixed point $`(X_R,Y_R,Z_R)`$. This is reflected in the $`(p,q)`$ subspace by the stray steps taken before the orbit moves onto the smooth, bounded ring of the regular attractor as shown in Fig. 9. This results in a regular oscillation in the mean-square displacement $`M(t)`$ that will eventually die away to give $`K0,`$ as shown in Fig. 7. The $`K`$-test reflects the regularity of the attractor and is unable to detect the subtle chaotic transient in the early motion.

By contrast, an orbit at $`r=30`$ in Fig. 7 is drawn onto the Lorenz strange attractor and its strongly chaotic behaviour is detected by the $`K`$-test. The chaotic motion is well reflected by the Brownian diffusion evidenced in the $`(p,q)`$ subspace shown in Fig. 10. This orbit leads to $`K1`$ as expected.

## IV Summary

We have confirmed that the $`K`$-test provides a simple and easy diagnostic for the transition from regularity to chaos for Hamiltonian as well as dissipative systems. For the test to be effective, an orbit must spend sufficient time on the hyperbolic chaotic attractor in a dissipative system, or be confined to a highly chaotic region of phase space if the system is Hamiltonian. However, for chaotic transients which move into a chaotic region of phase space and then out into a regular region of phase space the motion is not consistently Brownian. Consequently, the test can yield ambiguous and confusing results although a qualitative look over the $`(p,q)`$ subspace can provide guidance as to the transient irregularity of the base dynamics. Additionallly the K-test is not ideal for relativistic settings since it depends on the time coordinate used and so, like the Lyapunov exponents, is not a covariant indicator of chaos. These limitations of the $`K`$-test are not unexpected: no one probe of chaos can suit every scenario. Nor are they fatal to its utility. Rather, we hope that they will help to map out the territory over which this simple test is a reliable guide to the presence of chaos

## Acknowledgements

We are especially grateful to both G. Gottwald and I. Melbourne for answering questions regarding their method. JL thanks N.J. Cornish for useful conversations. JL is supported by a PPARC Advanced Fellowship and an award from NESTA.