text:
# Melting of the vortex lattice in high 𝑇_𝑐 superconductors

## Abstract

The precise measurements of vortex melting point towards a need to develop a quantitative theoretical description of thermal fluctuations in vortex matter. To tackle the difficult problem of melting, the description of both the solid and the liquid phase should reach the precision level well below $`1\%`$. Such a theory in the framework of the Ginzburg - Landau approach is presented. The melting line location is determined and magnetization, specific heat jumps along it are calculated. We find that the magnetization in the liquid is larger than that in the solid by $`1.8\%`$ regardless of the melting temperature, while the specific heat jump is about $`6\%`$ and slowly decreases with temperature. The magnetization curves agree with experimental results on $`YBCO`$ and Monte Carlo simulations.

A magnetic field generates an array of vortices in type II superconductors. The vortices strongly interact with each other forming highly correlated configurations such as the vortex lattice. In high $`T_c`$ cuprates at relatively high temperatures, vortices move and vibrate due to thermal fluctuations to the extent that the lattice can melt becoming a ”vortex liquid”. Several recent remarkable experiments clearly determined that the vortex lattice melting in high $`T_c`$ superconductors is a first order phase transition with magnetization jumps and spikes in specific heat . The magnetization and entropy jumps were measured using diverse techniques such as local Hall probes, SQUID , torque magnetometry and integration of the specific heat spike. Related investigations indicate that, in addition to the spike, there is also a jump in specific heat. In spite of those precise measurements of vortex melting, a quantitative theoretical description of vortex melting is still lacking. We present such a theory in the framework of the Ginzburg - Landau approach in this work.

Thermal fluctuations in vortex matter have attracted great attention since the high Tc superconductors were discovered over a decade ago. In highly anisotropic superconductors like $`BSCCO`$ near the melting point vortices are quite well separated and the system can be approximated by an array of pointlike objects. In less anisotropic ones like $`YBCO`$, near the melting point the vortices overlap and one has to use the Ginzburg - Landau (GL) model. The model is defined by the following free energy:

$$F=d^3x\frac{\mathrm{}^2}{2m_{ab}}\left|𝐃\psi \right|^2+\frac{\mathrm{}^2}{2m_c}|_z\psi |^2a(T)|\psi |^2+\frac{b^{}}{2}|\psi |^4+\frac{\left(𝐁𝐇\right)^2}{8\pi }$$
(1)
where $`𝐀=(By,0)`$ describes magnetic field (considered constant and nonfluctuating, see below) in Landau gauge and covariant derivative is defined by $`𝐃i\frac{2\pi }{\mathrm{\Phi }_0}𝐀,\mathrm{\Phi }_0\frac{hc}{e^{}}`$($`e^{}=2e`$).When $`\kappa =\lambda /\xi `$ is large (greater than $`10`$ ) where $`\lambda `$ magnetic penetration depth and $`\xi `$ coherence length, the magnetic field can be considered constant and nonfluctuating (an excellent approximation in the region studied). Statistical physics is described by the statistical sum $`Z=D\psi D\overline{\psi }\mathrm{exp}\left(\frac{F}{T}\right)`$. It accurately describes thermal fluctuations in the range of relatively large magnetic fields ($`H>>`$ $`H_{c1}`$ $`=100G`$ in $`YBCO`$) and temperatures near $`T_c`$ ($`70K`$ to $`110K`$ in $`YBCO`$ ). Near $`T_c`$, $`a(T)`$ can be approximated bt $`\alpha T_c(1t)`$. In this paper, we will consider only high $`\kappa `$ and temperatures near $`T_c`$. The model is, however, highly nontrivial even within the lowest Landau level (LLL) approximation. In this approximation, which is valid when the magnetic field is high, only lowest Landau level mode is retained and the free energy simplifies (after rescaling) :

$$f=\frac{1}{4\pi \sqrt{2}}d^3x\left[\frac{1}{2}|_z\psi |^2+a_T|\psi |^2+\frac{1}{2}|\psi |^4\right].$$
(2)
The simplified model has only one parameter, the dimensionless scaled temperature $`a_T=\left(\frac{b\omega }{4\pi \sqrt{2}}\right)^{2/3}a_h`$ , where $`\omega =\sqrt{2Gi}\pi ^2t,a_h=\frac{1tb}{2},t=T/T_c,b=B/H_{c2}`$. The dimensionless Ginzburg number $`Gi`$ characterizing the importance of thermal fluctuations is $`\frac{1}{2}\left(\frac{8\pi \kappa ^2\xi T_c\gamma }{\mathrm{\Phi }_0^2}\right)^2`$, and the anisotropy parameter $`\gamma `$ is $`\sqrt{m_c/m_{ab}}`$.

The LLL GL model was studied by a variety of different nonperturbative analytical methods. Among them the density functional, $`1/N`$ , elasticity theory  and others . The model was also studied numerically in both 3D  and 2D  However we will not discuss those approaches in this paper.

While applying the renormalization group (RG) on the one loop level to this model, Brezin, Nelson and Thiaville  found no fixed points of the (functional) RG equations and thus concluded that the transition to the solid phase is not continuous. The RG method therefore cannot provide a quantitative theory of the melting transition. Two perturbative approaches were developed and greatly improved recently to describe the solid phase and liquid phases, respectively. The perturbative approach on the liquid side was pioneered by Ruggeri and Thouless, who developed an expansion in which all the ”bubble” diagrams are resummed. Unfortunately, they found that the series are asymptotic and, although first few terms provide accurate results at very high temperatures, the series become inapplicable for $`a_T`$ less than $`2`$ which is quite far above the melting line (located around $`a_T=10`$). We obtained recently an optimized gaussian series which is convergent rather than asymptotic with radius of convergence of $`a_T5`$, but still above the melting line.

On the solid side, Eilenberger calculated the fluctuations spectrum around Abrikosov’s mean field solution. Maki and Takayama noticed that the vortex lattice phonon modes are softer than those of the usual acoustic phonons in atomic crystals and this leads to infrared (IR) divergences in certain quantities. This was initially interpreted as the destruction of the vortex solid by thermal fluctuations and the perturbation theory was abandoned. However the divergences resemble the ”spurious” IR divergences in the critical phenomena theory. A recent analysis demonstrated that all these IR divergences cancel in physical quantities and the series therefore are reliable. The two loop calculation was performed, so that the LLL GL theory on the solid side is now precise enough even for the description of melting.

However, on the liquid side, a theory for $`a_T5`$ is required. Moreover this theory should be extremely precise since the internal energies of the solid and the liquid near melting differ by a few percents only. Developing such a theory requires a better qualitative understanding of the metastable phases of the model. It is clear that the overheated solid becomes unstable at some finite temperature. It is not clear however whether the overcooled liquid becomes unstable at some finite temperature (like water) or exists all the way down to $`T=0`$ as a metastable state. Despite its limited precision, the gaussian (Hartree - Fock) variational calculation, is usually a very good guide to the qualitative features of the phase diagram. While such a calculation in the liquid was performed quite some time ago, on the solid side a more complicated one sampling inhomogeneous states was performed recently. The results are as follows. The solid state is the stable one below the melting temperature, becomes metastable at somewhat higher temperatures and is destabilized at $`a_T5`$. The liquid state becomes metastable below the melting temperature, however in contrast to the solid, it does not loose metastability all the way down to $`T=0`$ and the excitation energy approaches zero.

Meanwhile, similar qualitative results have been obtained in a different field of physics. A variety of analytical and numerical methods have indicated that liquid (gas) phase of the classical one component Coulomb plasma also exists as a metastable state down to $`T=0`$ with energy gradually approaching that of the Madelung solid and the excitation energy diminishing. We speculate that the same phenomenon would appear to happen in any system of particles interacting via long range repulsive forces. In fact the vortices in the London approximation resemble repelling particles with the force even more long range than the Coulomb. In light of this impetus to consider the above scenario in the vortex matter, we provide both theoretical and phenomenological evidence that the above picture is a valid one.

Assuming the absence of singularities on the liquid branch allows to develop a sufficiently precise theory of the LLL GL model in a vortex liquid (even including an overcooled one) using the Borel - Pade (BP) method at any temperature. After clarifying several issues that prevented its use and acceptance previously we then combine it with the recently developed LLL theory of solids to calculate the melting line and the magnetization and specific heat jumps along it. Early on, Ruggeri and Thouless tried unsuccessfully to calculate the specific heat by using BP because their series was too short. Subsequent attempts to calculate the melting point by using BP also ran into problems. Hikami, Fujita and Larkin (18) attempted to find the melting point by comparing the BP (liquid) energy with the one loop solid energy and, in doing so, obtained the melting temperature $`a_T^m=7`$. However, their one loop solid energy was incorrect and, in any case, the two loop correction is necessary. As demonstrated in the following, the BP energy combined with the correct two loop solid energy computed recently not only gives $`a_T^m=9.5`$, but also allows to obtain a wide range of quantitative predictions within the model.

Now we present the solution of the LLL GL model. The liquid LLL (scaled) effective free energy (of the scaled model defined in eq.(2)) is written as $`f_{liq}=4\epsilon ^{1/2}[1+g\left(x\right)]`$ . The function $`g`$ can be expanded as $`g\left(x\right)=c_nx^n`$ , where the high temperature small parameter $`x=\frac{1}{2}\epsilon ^{3/2}`$ is defined as a solution of the gaussian gap equation, $`\epsilon ^{3/2}a_T\epsilon ^{1/2}4=0`$ for the excitation energy . The coefficients can be found in ref. . We will denote $`g_k(x)`$ by the \[$`k,k1`$\] BP transform of $`g(x)`$ (other BP approximants clearly violate the correct low temperature asymptotics). The BP transform is defined as $`_0^{\mathrm{}}g_k^{}\left(xt\right)\mathrm{exp}\left(t\right)𝑑t`$ where $`g_k^{}`$ is the $`[k,k1]`$ Pade transform of $`_{n=1}^{2k1}\frac{c_nx^n}{n!}`$. For $`k=4`$ the liquid energy already converges. We used in this work $`k=5`$ to achieve the required precision ($`0.1\%`$). The liquid energy completely agrees with the optimized gaussian expansion results above its radius of convergence at $`a_T=5`$. In addition similar results we obtained in the 2D GL model agree with existing Monte Carlo simulations. The solid effective energy to two loops is:

$$f_{sol}=\frac{a_T^2}{2\beta _A}+2.848\left|a_T\right|^{1/2}+\frac{2.4}{a_T}.$$
(3)
Comparing the solid energy to that of the liquid (inset in Fig.1), reveals that $`a_T^m=9.5`$. This is in accord with experimental results. As an example, on Fig.2 we present the fitting of the melting line of fully oxidized $`YBa_2Cu_3O_7`$  which gives $`T_c=88.2,H_{c2}=175.9,Gi=7.0`$ $`10^5,\kappa =50`$. Melting lines of optimally doped untwinned  $`YBa_2Cu_3O_{7\delta }`$ and $`DyBa_2Cu_3O_7`$  are also fitted extremely well. For example, for $`YBa_2Cu_3O_{7\delta }`$ in , the melting line fitting gives $`T_c=93.07,H_{c2}=167.53,Gi=1.9`$ $`10^4,\kappa =48.5`$ (see also ref.  ).

3D Monte Carlo simulations  are not precise enough to provide an accurate melting point since the LLL scaling is violated. One gets different values of $`a_T^m=14.5,13.2,10.9`$ at magnetic fields $`1,2,5`$T respectively. This violation is perhaps due to a small sample size (of order 100 vortices). The situation in 2D is better since the sample size is much larger. We performed similar calculation in the 2D LLL GL model and found that the melting point is $`a_T^m=13.2`$ . It is in good agreement with the MC simulations. Phenomenologically the melting line can be located using Lindemann criterion or its more refined version using Debye - Waller factor. The more refined criterion is required in the case of $`YBCO`$ since vortices are not point-like. Numerical investigation of the Yukawa gas indicated that the Debye - Waller factor $`e^{2W}`$ (a ratio of the structure function at the second Bragg peak at melting to its value at $`T=0`$) is about $`60\%`$. We get using methods of ref.  $`e^{2W}=0.59`$ for $`a_T^m=9.5`$.

The scaled magnetization is defined by $`m\left(a_T\right)=\frac{d}{da_T}f_{eff}\left(a_T\right)`$. At the melting point $`a_T^m=9.5`$ the magnetization jump ratio defined by $`\mathrm{\Delta }M`$ divided by the magnetization at the melting on the solid side is found to be equal to

$$\frac{\mathrm{\Delta }M}{M_s}=\frac{\mathrm{\Delta }m}{m_s}=.018.$$
(4)
This prediction is compared on Fig.2 (the upper inset) with the experimental results on fully oxidized $`YBa_2Cu_3O_7`$ (ref. , rhombs) and optimally doped $`YBa_2Cu_3O_{7\delta }`$ (ref. , stars) . These samples probably have the lowest amount of disorder which is not included in the calculations. From the model, we caculate the specific heat jump ratio at the melting

$$\mathrm{\Delta }c=0.0075\left(\frac{22b+t}{t}\right)^20.20Gi^{1/3}\left(b1t\right)\left(\frac{b}{t^2}\right)^{2/3}.$$
(5)
It is compared on the lower inset on Fig.2 with the experimental values of ref.  (using the fitting parameters given above).

In addition to describing the melting, we present here an example of quantitative results that are obtained using the present approach - the magnetization curves. Our LLL magnetization curve coincides with the LLL Monte Carlo result of ref.  (which is very accurate since the LLL scaling is obeyed) to the precision of MC. However in experiments away from the melting line higher Landau levels (HLL) are no longer negligible. Naively in vortex solid when the distance from the mean field transition line is smaller than the inter - Landau level gap, $`1tb<2b`$, one expects that the higher Landau modes can be neglected. More carefully examining the mean field solution reveals that a weaker condition $`1tb<12b`$ should be used for a validity test of the LLL approximation  in vortex solids. In vortex liquid one has to go beyond the mean field to estimate the HLL contribution . In 3D Lawrie in  calculated the excitation energy in the framework of the gaussian (Hartree - Fock) approximation. The excitation gap is 10 times smaller than inter Landau gap for fields in a wide range around melting line for fields larger than $`0.1T`$ in $`YBCO`$. Therefore in the range of values of the interest in the present paper the LLL contribution should be dominant. Experimentally it is often claimed that one can establish the LLL scaling for fields above $`3`$ T (see, for example, ref. ).

The theoretical expressions we use are the LLL contribution to the magnetization plus the corrections due to HLL calculated in gaussian approximation. The results are compared on Fig.3 with the experimental magnetization curves of ref. . We use the parameters from the fitting of the melting curve (see Fig.2) . The agreement is fair at higher fields, while at low magnetic fields the higher Landau levels theory beyond gaussian approximation is required. The LLL scaling has a limited validity away from melting line.

To summarize, the problem of the quantitative description of melting of the vortex lattice in the framework of the LLL Ginzburg - Landau approach is solved. The results for melting line, magnetization jump and specific heat jump are in good agreement with experiments and MC simulations. This is the first quantitative theory of the first order melting of any kind to our knowledge. We believe that similar methods can be applied to other systems undergoing the first order melting transition.

We are grateful to T. Nishizaki and A. Junod for providing details of their experiments. The work was supported by NSC of Taiwan grant NSC$`\mathrm{\#}`$90-2112-M-009-039.


format: text

text:
# HEP Applications on the Grid Canada Testbed

## I Introduction

One of the motivations for establishing a computational Grid is the ability to exploit unused computing cycles at remote facilities. In Canada, there has been a significant injection of new funding for mid-size capability and capacity computational facilities. One of the conditions of the funding is that these facilities must make 20% of the resources available to outside users. A number of these facilities are associated with high energy physics (HEP) centers or groups. However, the majority are facilities that are shared between many research fields. We would like to exploit these unused cycles for HEP applications such as Monte Carlo simulations. This paper describes the Grid Canada Testbed which consists of resources at both HEP and non-HEP sites. We demonstrate that HEP simulations can be run without installing application-specific software at the remote site. We discuss the results of a series of tests we have run over the Testbed and conclude with a discussion of our plans for the future.

## II Grid Canada Testbed

Grid Canada (www.gridcanada.ca) is an organization established by CANARIE Inc. (www.canarie.ca), who provide the national research and educational network in Canada; the C3.ca Association (www.c3.ca), which is an organization that represents the user community of High Performance Computing (HPC) facilities in Canada; and the National Research Council of Canada (NRC) (www.nrc.ca).

The goal of Grid Canada (GC) is to foster collaboration on Grid activities in Canada. GC is responsible for the Canadian Certificate Authority and has helped coordinate a number of demonstration activities such as the 1 TB high-speed file transfer project between the TRIUMF Laboratory in Vancouver and CERN igrid  and a grid test using a number of HPC (shared memory machines) computers across Canada. In addition, GC has helped to establish the Linux Testbed that will be described in this paper.

The GC Linux Testbed involves researchers at the University of Victoria who are members of the BaBar and ATLAS Collaborations. In addition, personnel at the NRC Institute for Information Technology, CANARIE Inc. and the British Columbia Ministry for Management Services have also made significant contributions to the Testbed.

The Testbed currently includes 12 sites in Canada ranging from a number of sites in British Columbia to those 3000 km away in Ontario. Typically each site has between 1 and 8 computers. All sites are required to run Globus 2.0 or 2.2 globus  on their computers except for the Victoria site which is described in the next paragraph. In addition, sites must also run an OpenAFS client openafs . Latencies between the two most-distant sites vary between 50 and 100 milliseconds. The remote sites access the application-specific software and, for some tests, read the input and write the output data via AFS.

The Victoria Grid testbed is shown in fig 1. It consists of two servers for AFS and Objectivity objy . A head node (grid.phys.uvic.ca) runs the Globus software with a series of processing nodes (grid1, …, grid6) hidden behind a firewall. Jobs are distributed to the processing nodes using Condor condor  although PBS pbs  has been tested. The processing nodes have the ability to access the AFS Server through the firewall.

Access to the software via AFS by many remote machines can cause significant degradation in the CPU utilization to the Server and subsequently reduce the efficiency of the remote processing machines. AFS has a faster, more efficient method of reading clones or read-only volumes. For example, if a volume is read-only, then the client only has to contact the AFS server once to cache the volume. We tested AFS performance using two clones and two replicated volumes on two separate machines. The AFS Server machine (basilisk.phys.uvic.ca) held the master and a cloned volume, with the Objectivity Server machine (objy.phys.uvic.ca) holding a replicated and a cloned volume. This configuration appeared to be sufficient for the size of the Testbed so that no bottlenecks occurred as a result of the AFS Servers.

Initial tests of the ATLAS software using a local machine and a machine 1200 km distant showed that using a cloned volume increased the access to the software at the remote machine. Real execution time at the remote machine dropped by about 40% whereas (as expected) no change was observed at the local machine. Running simultaneous multiple jobs over the Grid did not degrade the AFS performance. We concluded that having cloned volumes significantly increased the AFS performance. We plan to continue to monitor the AFS performance as we continue to scale up the size of the GC Testbed.

Jobs were submitted from Victoria using simple scripts. No resource broker or portal was used however, we are examining the available tools (such as GridPort gridport ) and plan to incorporate them into the Testbed. Monitoring of the GC Testbed was done using the Cacti package cacti  providing network traffic and CPU utilization.

Detailed descriptions of the Victoria Testbed can be found at http://grid.phys.uvic.ca.

## III Description of Tests

HEP simulation applications begin by the generation of events which are collections of 4-vectors representing the trajectories of particles from an initial collision of two incident particles. The 4-vectors are then passed through a program that simulates the response of the experimental detector. The output is in a format that is usually identical to the format used to store the real collision data, but includes additional simulation information. After the detector simulation, the data is passed through an application that reconstructs the data into physical quantities such as tracks in the tracking chambers and energy in the calorimeters. It is now common practice to inject background events (either real or simulated) on top of the simulated physics events.

We chose to examine the performance of the GC Testbed using the simulations of the BaBar and ATLAS experiments. The main difference between the ATLAS and BaBar simulation applications is the manner in which they store their data. ATLAS uses conventional files in a Zebra format developed at CERN while BaBar uses the commercial Objectivity database. In the following subsections we describe the results of running the BaBar and ATLAS simulations on the GC Testbed. It is worth pointing out that neither the BaBar or the ATLAS simulation are optimized for operation over a wide-area Grid.

### III.1 Tests using the BaBar Monte Carlo Simulation

The BaBar Collaboration studies electron-positron collisions at the SLAC facility. In addition to collecting data from the experiment, the collaboration generates a significant amount of Monte Carlo simulated data in order to model the response of the detector and to help assess efficiencies and systematic errors in the physics studies.

The simulation application is run in 3 phases for event generation, detector simulation and reconstruction. Data is written to Objectivity by the event generation application. The detector simulation reads the output of the event generation and writes output that is then used by the reconstruction. The reconstruction application writes out data that can be used for physics studies.

The BaBar application accesses the database multiple times per event. The handshaking required between the database and application combined with the large latency resulted in relatively poor utilization of the remote processors. For example, we observed that CPU utilization was approximately 10% at sites 3000 km away from the Objectivity Server.

A series of 1 day tests were run involving many of the GC Testbed sites. A total of approximately 50,000 events were generated per test. The BaBar software was accessed from the AFS Server. The network traffic to the Objectivity and AFS Servers is shown in figure 2. Approximately 0.5 MB/s of data was constantly read and written to the database. AFS traffic was generally only a few KB/s with a peak during the start of the application.

Our initial tests showed that one of the areas that had to be addressed was the issuing and releasing of database locks. For example, a remote application locks the entire federation during a creation of a container. If a remote machine lost network connectivity while a global database lock was present then this resulted in the entire Testbed being halted as none of the sites would be able to write to the database. As a result we developed a series of monitoring and cleanup utilities that would eliminate any problems associated with the database locks.

We consider the test a success in that we could run the BaBar application on non-HEP resources. It is clear that the efficiency at the remote sites was low and we believe a number of modifications to the BaBar application could improve the throughput. However, given the plan by BaBar and the LHC experiments to move away from the Objectivity database, it was decided not to pursue this application further.

### III.2 Tests using the ATLAS Data Challenge software

The ATLAS Collaboration will study high energy proton-proton collisions at the Large Hadron Collider at the CERN Laboratory. ATLAS is expected to produce large volumes of data starting in 2007. The ATLAS Data Challenge (DC) is a project to test and develop the ability to exploit remote computing resources for Monte Carlo Simulation. For our tests, we used the application that generates background (called pile-up) events. Input and output data are standard files in a Zebra format. Future ATLAS DC will use the Root data formatroot ; such a change will make no difference to the way in which the current test was performed.

We tested three modes of operation on the testbed:

1. Both input data and output data transferred using AFS
2. Input data transferred using GridFTP and output data transferred using AFS
3. Both input and output data transferred using GridFTP

All software and log files use AFS.

It has been recognized that AFS was not an ideal protocol for transfering large data files, however, we felt it was important to quantify the differences. We found AFS to extremely slow for reading input data and resulted in CPU utilizations being less than 5%. AFS caches its output; as a result we tested a mode where the input was sent via GridFTP to the remote site and let AFS control the output. The results of this mode of operation were more encouraging, however the performance of AFS is comparable to single stream FTP and using GridFTP is a much more effective solution.

The mode of operation we chose was one where the input files were tarred and compressed, and sent via GridFTP to the remote site. At the remote site, the files would be untarred and uncompressed, and the job executed. At the end of the job, the output files would be tarred and compressed, and then transferred via GridFTP back to the Victoria disk server. This was found to be a successful mode of operation. Access to software via AFS proved to be reliable and efficient. We chose an application that would run extremely quickly on the processing nodes (30 minutes on a 2 GHz CPU). Compression, transfer and uncompression was approximately 20 minutes for a 2 GB file. Although the overall utilization of a remote CPU could be considered to be low if one includes the compression and transfer times, a more realistic application, such as simulation or reconstruction, would take hours rather than minutes to process data files of the same size. In this situation, we believe the CPU utilization would significantly improve.

## IV Summary

We have demonstrated that HEP applications can be run on generic computing resources with relatively high efficiency. With (modest) design changes in the applications we expect that further gains in efficiency would be possible. Further tests are planned using a sites with significantly larger resources with the aim of establishing a production level Grid resource in Canada.

###### Acknowledgements.

We would like to thank those people who made available resources at their sites. The support of the Natural Sciences and Engineering Research Council of Canada, the National Research Council of Canada, the C3.ca Association and CANARIE Inc. is aknowledged.


text:
# Untitled Document

COMMISSIONS 27 AND 42 OF THE IAU

INFORMATION BULLETIN ON VARIABLE STARS

Number xxxx

Konkoly Observatory Budapest xx xxx 2003 HU ISSN 0374 – 0676

NEW R CRB-TYPE STAR HadV98

KATO, TAICHI1; KATSUMI, HASEDA2

1 Dept. of Astronomy, Kyoto University, Kyoto 606-8502, Japan, e-mail: tkato@kusastro.kyoto-u.ac.jp

2 Variable Star Observers League in Japan (VSOLJ), 2-7-10 Fujimidai, Toyohashi City, Aichi 441-8135, Japan, e-mail: khaseda@mx1.tees.ne.jp

HadV98 = CoD $``$2212017 = GSC 6825.253 is a variable star discovered by K. Haseda (Haseda and Kato 2001). The variable had been almost constant until 2001 April, when it suddenly started fading. The object was suspected to be an R CrB-type star, from the time of the discovery alert, based on its light behavior and its Tycho-2 moderate color ($`BV`$1.0). This star has been very recently confirmed to be a genuine R CrB-type variable from spectroscopy (Hasselbach et al. 2002).

We studied the variability of this object using Haseda’s discovery material combined with the ASAS-3 survey data (Pojmanski 2002). The light curve constructed from Haseda’s observations and the ASAS-3 public data is presented in Figure 1. Haseda’s observations were performed with a $`D`$ = 10 cm f/4.0 telephoto lens and unfiltered T-Max 400 emulsions. The passband of observations covers the range of 400–650 nm. The photographic magnitudes were estimated using GSC stars. A 0.08 mag systematic correction was added to Haseda’s photographic observations, in order to best match the ASAS-3 $`V`$ magnitudes. The recorded range of variability in $`V`$ was 10.7–13.4, but the true minimum seems to have been fainter.

The light curve shows a rapid decline and a slower recovery, which are characteristic to R CrB-type optical light curves (Clayton 1996). The 2001–2002 fading was composed of at least two distinct minima (Figure 2), one lasting in 2001 May – 2001 July (JD 2452030–2452130) and another minimum suddenly starting from 2001 October (JD 2452190–), whose entire picture was not well observed because of the solar conjunction. Such an occurrence of closely spaced two minima are also sometimes seen in other R CrB-type fadings and other dust-forming objects (Kato et al. 2002).

We also searched for a pulsation signal, which is frequently seen in many R CrB stars, most notably recorded in RY Sgr and V854 Cen (cf. Clayton 1996), using the ASAS-3 observations at maximum between JD 2452697 and 2452762 (20 measurements). No remarkable pulsation was found with an amplitude larger than 0.1 mag. This finding is consistent with the lack of remarkable variability in Haseda’s observations before 2001. HadV98 apparently belongs to a group of R CrB stars showing less distinct pulsations at maximum.

We are grateful to G. Pojmanski for making the ASAS-3 survey data publicly available, and generously allowing us for unlimited usage. This work is partly supported by a grant-in-aid (13640239, 15037205) from the Japanese Ministry of Education, Culture, Sports, Science and Technology.

References:

Clayton, G. C., 1996, PASP, 108, 225

Haseda, K., Kato, T., 2001, vsnet-alert 6226
(http://www.kusastro.kyoto-u.ac.jp/vsnet/Mail/alert6000/msg00226.html)

Hasselbach, E., Clayton, G. C., Smith, P. S, 2002, AAS Meeting 201, 17.11

Kato, T., Haseda, K., Takamizawa, K., Yamaoka, H., 2002, A&A, 393, 69L

Pojmanski, G. 2002, Acta Astronomica, 52, 397


text:
# Controlled coupling of selected single quantum dots with high-Q microdisk microcavities

## Conclusions

A technique to position single quantum dots in microdisks has been developed. Optical studies have confirmed the presence of only one QD in the microdisk. One major application of this technique is the possibility to compare theoretical and experimental measurements on the lifetime modifications, as the coupling of the dot to the field depends strongly on the dot position. Further studies should enable us to reach the strong coupling regime for a single dot in a microdisk.

## Acknowledgments

Part of this work was performed within the Nanometer Structure Consortium in Lund, Sweden and was supported by NFR, TFR, NUTEK and SSF. One of us (VZ) acknowledges funding from the European Union’s Marie Curie fellowship.


text:
# Pressure induced transition from a spin glass to an itinerant ferromagnet in half doped manganite Ln0.5Ba0.5MnO3 (Ln=Sm and Nd) with quenched disorder

## Abstract

The effect of quenched disorder on the multiphase competition has been investigated by examining the pressure phase diagram of half doped manganite Ln0.5Ba0.5MnO3 (Ln = Sm and Nd) with A-site disorders. Sm0.5Ba0.5MnO3, a spin glass insulator at ambient pressure, switches to a ferromagnetic metal with increasing pressure, followed by a rapid increase of the ferromagnetic transition temperature TC. The rapid increase of TC was confirmed also for Nd0.5Ba0.5MnO3. These observations indicate that the unusual suppression of the multicritical phase boundary in the A-site disordered system, previously observed as a function of the averaged A-site ionic radius, is essentially controlled by the pressure and hence the band width. The effect of quenched disorder is therefore much enhanced with approaching the multicritical region.

One of the hallmarks of strongly correlated electron physics is the multiphase competition produced by a delicate interplay among charge, spin and orbital degrees of freedom ref1 ; ref2 . An appealing example of such multicriticality can be seen in perovskite manganites. Competing interactions/orders inherent in the manganites, such as double-exchange ferromagnetism vs. super-exchange antiferromagnetism and charge-orbital order vs. metallic state, tend to produce the multicritical state where external stimuli occasionally cause the dramatic phase conversion, e.g., between a metal and an insulator, or between a ferromagnet and an antiferromagnet ref3 . The colossal magnetoresistance (CMR) phenomena ref4  are believed to be a consequence of such phase competitions.

Recently, the vital role of quenched disorder in the critical region of manganites has been attracting considerable interest. The quenched disorder ref5  may result in a phase separation ref6  into the competing two ordered phases on various time-scales and length-scales, which substantially modifies the criticality and hence the response to the external field. The first suggestion came from the observation of relaxor-like behavior induced by Cr-doping on the Mn sites in the charge/orbital ordered (CO/OO) state ref7 ; ref8 . This is highly likely associated with the ferromagnetic (FM) metallic clusters coexisting with the background CO/OO state. The effect of quenched disorder on the A-site was subsequently examined, which shows a drastic modification of the phase diagram as a function of the average A-site ionic radius.

The half-doped manganites Ln0.5Ba0.5MnO3 (Ln = Lanthanides) are one of the ideal systems to explore the physics of quenched disorder. The perovskite manganites in this chemical composition have two possible forms of the crystal structure depending on synthetic condition ref9 . One is with the A-sites randomly occupied by Ln3+ and Ba2+ ions. The melt-quenched sample shows the complete solid solution of Ln and Ba ions on the A-sites with a simple cubic structure as the average structure. The other is the A-site ordered perovskite with the alternate stack of LnO and BaO sheets along the c-axis with intervening MnO2 sheets. This is due to the large difference in the ionic radii of Ln3+ and Ba2+. The respective MnO2 sheets in this tetragonal form are free from random potential which would otherwise arise from the Coulomb potential or local strain from the random occupation of Ln/Ba on the A-sites.

The phase diagram ref10 ; ref11  of ordered Ln1/2Ba1/2MnO3 as a function of ionic radius is schematically shown by the dotted line in Fig. 1. The increase of ionic radius results in the reduction of Mn-O bond bending and hence the increases of the band width. With relatively small Ln ions from Dy to Sm, the CE-type charge ordered state with diagonal orbital stripes develops at low temperatures. With larger Ln ions from Nd to La, on the other hand, the ferromagnetic metal phase and, eventually, the A-type antiferromagnetic phase with ordering of $`d_{x^2y^2}`$ orbitals are stabilized upon cooling the sample. These phases meet with each other at around Ln = Nd, forming a well defined multicritical point. When the solid-solution of Ln3+ and Ba2+ is formed, however, this multicritical region is drastically altered ref11 . The solid line in Fig. 1 indicates that the phase diagram of A -site disordered Ln0.5Ba0.5MnO3. It is clear that the phase transitions near the multicritical point are suppressed significantly. In contrast to the ordered Ln0.5Ba0.5MnO3, the charge and orbital ordered phase goes away and alternatively a spin glass phase without long-range order in any charge and orbital sectors emerges. With increasing the ionic radius, the ferromagnetic metal phase takes over the spin glass phase around Nd. The A-type antiferromagnetic phase, however, does not show up at low temperatures. Besides, as compared with the A-site ordered systems, the ferromagnetic phase extends to the smaller ionic radius and the ferromagnetic transition temperature, TC, is appreciably suppressed near the critical region.

The contrast between the A-site ordered and the disordered systems demonstrates a substantial influence of quenched disorder on the multiphase competition in the manganites. However, it should be emphasized here that the control parameter of the disordered system is not only the bandwidth but also the degree of disorder. By changing the Ln ion, the degree of disorder is inherently modified associated with the change in the ionic radius. It is difficult to identify which factor plays a dominant role in the unusual phase diagram of the disordered system. For example, the rapid decrease of TC in the ferromagnetic phase with approaching the critical region can be interpreted as an enhanced influence of quenched disorder in the critical region. Alternately, it may be understood in terms of the increased disorder with decreasing the ionic radius of Ln ions from those comparable with large Ba ions.

To separate the two contributions arising from the disorder and the bandwidth, we have employed hydrostatic pressure to modify the bandwidth alone, while keeping the degree of disorder constant. We have selected two A-site disordered compounds, Sm0.5Ba0.5MnO3 and Nd0.5Ba0.5MnO3, and applied hydrostatic pressures up to 10 GPa. We found a pressure induced spin glass insulator to ferromagnetic metal transition, which can be mapped onto the phase diagram obtained by examining the solid-solution effect. This clearly demonstrates that the unusual phase diagram of disordered Ln0.5Ba0.5MnO3 arises from the bandwidth.

The A-site solid-solutions of Ln0.5Ba0.5MnO3 (Ln = Sm and Nd) were grown by the floating-zone method in a single crystalline form. The starting polycrystalline rods were prepared by a standard solid-state reaction. The x-ray powder diffraction patterns of the resultant crystalline rods indicate the formation of single phase Ln0.5Ba0.5MnO3 without any evidence for the A-site ordering. The resistivity and the magnetization data of disordered Sm0.5Ba0.5MnO3 and Nd0.5Ba0.5MnO3, at ambient pressure, are shown in the right panel of Fig. 1. The resistivity of disordered Sm0.5Ba0.5MnO3 shows an insulating behavior down to the lowest temperature measured and no evidence for the long-range charge ordering is observed. Corresponding to the phase line TSG in Fig. 1, we observe a spin glass transition at 50 K in the magnetization but no appreciable anomaly at TSG is observed in the resistivity. In the Nd0.5Ba0.5MnO3 sample, the resistivity shows a very rapid decrease below 140 K, which corresponds to the transition to a ferromagnetic metal phase. Note that the Curie temperature here is reduced to almost a half of the cation ordered Nd1/2Ba1/2MnO3 ref10 ; ref11 . Hydrostatic pressures up to 10 GPa were applied to those two crystals by a cubic-anvil-type pressure apparatus ref12 . The electrical resistivity and the AC susceptibility measurements were conducted within the pressure cell down to T = 4.2 K.

In the top panel of Fig. 2, the results of resistivity measurements on Sm0.5Ba0.5MnO3 under various pressures were summarized. Up to 4 GPa, no appreciable change was observed except for the decrease of the magnitude of resistivity. The resistivity shows an insulating behavior down to 4.2 K, indicating that the system remains in the spin glass phase. Above 4 GPa, in contrast, a substantial change was observed at low temperatures. The resistivity shows a rapid decrease upon cooling, which is a characteristic behavior of the transition to a ferromagnetic metal. This clearly indicates that the system has experienced a pressure induced spin glass insulator to a ferromagnetic metal transition. Above the critical pressure, the transition temperature, TC, goes up very rapidly with increasing pressure. These behaviors can be visually summarized as a phase diagram in Fig. 3. We immediately notice that the phase changes as a function of pressure are essentially the reproduction of those with the average ionic radius on the A-site.

It is now clear that the dominant controlling factor of the unusual phase diagram in disordered Ln0.5Ba0.5MnO3 is the bandwidth. Then, the pressure dependence of transition temperatures in disordered Nd0.5Ba0.5MnO3 should be connected to that of Sm0.5Ba0.5MnO3 by shifting the pressure axis. The bottom panel of Fig. 2 indicates the evolution of the temperature dependent resistivity and the AC susceptibility with increasing pressure for the A-site disordered Nd0.5Ba0.5MnO3. TC defined as a temperature of resistivity anomaly shows a rapid increase with increasing pressure. At high pressures, the resistivity anomaly tends to become weaker, which make the determination of TC ambiguous. In the AC susceptibility, however, we are able to see the well defined anomalies associated with the ferromagnetic transition (indicated by arrows), which provides us with a solid estimate of TC. Thus obtained pressure dependence of TC is plotted in Fig. 3. By shifting by 10 GPa, the TC-P curve of disordered Nd0.5Ba0.5MnO3 seems to be smoothly connected to that of Sm0.5Ba0.5MnO3. This provides a further evidence for the dominance of pressure as the controlling parameter.

Another important point seen in Fig. 2 is that TC of Nd0.5Ba0.5MnO3 tends to saturate with increasing pressure. In the saturation region, TC recovers to almost 300 K, which is close to those observed for ordered Ln1/2Ba1/2MnO3 with Ln = Nd and Pr. These observations strongly suggest that the suppression of TC is much reduced once away from the critical region. Note again that the degree of disorder is kept constant in this experiment. We therefore conclude that the effect of quenched disorder is substantially enhanced when coupled with the (multi)criticality.

Recently, it was pointed out that the quenched disorder triggers a phase separation into ferromagnetic and charge ordered domains ref13 . In their scenario, the formation of a clustered state with randomly oriented ferromagnetic clusters results in the decreases of the clean limit TC to a much lower TC in the ferromagnetic region. This accounts for the anomalous suppression of TC with approaching the critical region ref14 . More recent works claim that the state below the clean limit TC is not statically clustered but homogeneous state with substantial critical charge/lattice fluctuations induced by quenched disorder ref15 . Present pressure measurement alone cannot specify whether the phase separation is static or dynamic.

In summary, we have observed a pressure-induced switching of a spin glass into a ferromagnet in a half-doped manganite Sm0.5Ba0.5MnO3 with quenched disorder, which is followed by a rapid increase of Curie temperature upon pressure. The remarkably rapid increase of Curie temperature was confirmed also in itinerant ferromagnet Nd0.5Ba0.5MnO3 with a slightly larger ionic radius. Thus obtained pressure phase diagram can be mapped onto the unusual phase diagram of disordered Ln0.5Ba0.5MnO3 as a function of ionic radius, indicating that the critical control parameter is the bandwidth under the presence of quenched disorder. It is now firmly established that the effect of quenched disorder is much more enhanced in the critical region, resulting in a stronger suppression of phase transition with approaching the multicritical point. The high pressure measurements have proved themselves as a remarkably useful probe to explore the physics of quenched disorder in multi-phase competing systems.

The authors would like to thank Y. Motome, N. Nagaosa and Y. Tomioka for enlightening discussion. This work is partly supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology, Japan and New Energy and Industrial Technology Development Organization (NEDO), Japan.


text:
# Efficient mixed-force first-principles molecular dynamics

## Abstract

We present an efficient method to mix well converged ab initio forces with simpler and faster ones in molecular dynamics. While the cheap forces are evaluated every time step, the converged ones correct the trajectory only every $`n`$ time steps. For convenience, both types of forces are calculated with the same basic scheme, using density functional theory, norm conserving pseudopotentials, and a basis set of numerical atomic orbitals. The cheap forces are evaluated with a short-range minimal basis set and the non-selfconsistent Harris functional. Since these evaluations are hundreds of times faster than those of the converged forces, they add a neglegible cost, and the boost in computational efficiency is approximately a factor $`n`$. Our results indicate that one can use values of $`n`$ of up to 10, without affecting significantly the calculated structural and dynamical magnitudes.

Molecular dynamics (MD) is a fundamental tool in atomistic materials simulation Allen and Tildesley (1987). A majority of practicioners have used classical, semiempirical interatomic potentials. This is necessary for the large sizes and long times required to simulate many processes of enormous scientific and technological interest, from materials deformation and fracture Bulatov et al. (1998) to protein folding Duan and Kollman (1998). A large effort has been devoted to develop interatomic potentials for many types of systems Voter (1996). However, the quantitative reliability of such potentials in situations of bond formation and breaking is highly questionable. In such cases, it is imperative to use the much more expensive ab initio MD methods Car and Parrinello (1985); Payne et al. (1992), generally limited to a few hundred atoms and a few tens of picoseconds. Thus, it is essential to find methods that accelerate the integration of the dynamical equations, thus allowing for longer simulations.

In classical dynamics, one of such methods Tuckerman et al. (1992) uses multiple time scales to integrate the equations of motion for systems with both fast and slow dynamical degrees of freedom. The same method can be used to compute separately the hard, short-ranged forces from the soft, long-ranged ones. De Vita and Car De Vita and Car (1998) have proposed to adapt ‘on the fly’ the parameters of a classical potential using sporadic or periodic evaluations of ab initio forces. In this work, drawing ideas of those previous works, we propose a simple method to speed up dramatically ab initio MD. In principle, it could be implemented by combining classical and ab initio forces. However, such an approach would still require to develop a suitable classical force field for every new system with different interactions. Therefore, instead we take advantage of the fact that, while standard density functional forces require the simultaneous convergence of many parameters, much lower values of those parameters can still yield quite reasonable forces. Thus, by reducing drastically the size of the basis set, the Brillouin zone sampling, or the number of selfconsistency iterations, it is possible to reduce the computer time by enormous factors, and still obtain forces which are considerably more reliable than those of classical interatomic potentials.

To test our scheme, we have chosen the SIESTA method Ordejón et al. (1996); Soler et al. (2002), which is specially well suited to span the range from ‘quick and dirty’ calculations to fully converged ones. It uses density functional theory Kohn and Sham (1965) (DFT), norm-conserving pseudopotentials Troullier and Martins (1991) and a basis set of numerical atomic orbitals of strictly finite range Sankey and Niklewski (1989); Anglada et al. (2002). To calculate converged forces, we might typically use the generalized gradient approximation (GGA) to exchange and correlation (with spin polarization if required), double-$`\zeta `$ polarized (DZP) basis orbitals with a relatively long range, fine integration grids in real and reciprocal space, and a well converged selfconsistency between density and potential. For the cheap forces we may save on many different parameters, depending on the system and the properties studied. Thus, we may use the local density approximation (LDA), a minimal single-$`\zeta `$ basis set with short range, a coarser integration grid in real space, just the $`\mathrm{\Gamma }`$ point in reciprocal space, and the non-selfconsistent Harris functional Harris (1985). All together, the cheap forces are typically hundreds of times faster to compute than the converged ones, and therefore they add a neglegible cost to the overall calculation, thus making unneccessary to resort to classical force fields.

As usual, we use the Born-Oppenheimer approximation and we treat the nuclei as classical particles, subject to the Hellmann-Feynman forces (including all Pulay corrections). The equations of motion are solved with the standard velocity-Verlet algorithm Allen and Tildesley (1987), what ensures the time reversibility of the trajectories Tuckerman et al. (1992). The atomic forces at time $`t`$ are defined as $`𝐅_{fast}(t)+\mathrm{\Delta }𝐅(t)`$, where

$$\mathrm{\Delta }𝐅(t)=\{\begin{array}{cc}n(𝐅_{conv}(t)𝐅_{fast}(t))\hfill & \text{if }(t/\mathrm{\Delta }tmodn)=0\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(1)
Thus, the expensive converged forces $`𝐅_{conv}`$ need to be evaluated only once every $`n`$ time steps $`\mathrm{\Delta }t`$. In those ‘correction steps’, the trajectories generated by the cheap (fast) forces $`𝐅_{fast}`$ are corrected by applying a force ‘kick’ equal to the difference between the converged and fast forces at that time, multiplied by $`n`$. The factor $`n`$ accounts for the concentration of the continuous force correction $`𝐅_{conv}(t)𝐅_{fast}(t)`$ in one out of every $`n`$ steps. The method of Ref. Batcho and Schlick, 2001, based on the position-Verlet algorithm, was reported to have a better numerical stability in response to the correction ‘kicks’. The efficiency of that method in the present context will be studied in future works.

Figure 1 shows schematically the positions, velocities and forces of a particle moving in one dimension, generated with our mixed-force algorithm.

For simplicity, we take the converged force and the initial velocity equal to zero, so that the correct converged position is also zero at all times. The mixed-force trajectory, for a constant negative fast force, shows periodic force kicks that change discontinuously the velocity and invert the trajectory at the correction steps.

We have applied this method to simulate a system of 64 silicon atoms at an average temperature of $`2000`$ K and an average pressure close to zero. This high temperature was intentionally chosen to test the method under specially stringent conditions, whith high kinetic energies and frequent formation and breaking of bonds. The simulations were performed with the SIESTA program Soler et al. (2002) but standard Hamiltonian diagonalizations were used instead of order-$`N`$ methods Ordejon et al. (1993); Kim et al. (1995), because of the metallic character of liquid silicon. For the cheap forces we use the Harris functional, a minimal basis set with a range of 3.5 and 4.0 Bohr for $`s`$ and $`p`$ orbitals, a real-space integration grid with a plane wave cutoff of 40 Ry, and only the $`\mathrm{\Gamma }`$ $`k`$-point. For the converged forces, we use the self-consistent Kohn-Sham functional in the LDA, a double-$`\zeta `$ polarized (DZP) basis set with a range of 5.4, 6.5, and 3.8 for $`s,p`$ and $`d`$ basis orbitals, a real-space grid with a 80 Ry plane wave cutoff, and only the $`\mathrm{\Gamma }`$ $`k`$-point. The forces are corrected according to equation (1) every $`n`$ time steps, with $`\mathrm{\Delta }t=1`$ fs.

Figure 2 compares the magnitudes of the fast and converged forces, and of their difference.

It can be seen that the latter is a relativelly small and smooth correction, what explains why it may be evaluated and applied less frequently.

Figure 3 represents the divergence of the trajectories, generated with different values of $`n`$, away from the reference converged trajectory (which corresponds to $`n=1`$ in Eq. (1)).

It can be seen that even the $`n=10`$ trajectory diverges much more slowly than that obtained purely from the fast forces (labelled $`n=\mathrm{}`$, i.e. with no force corrections).

Figure 4 shows the total energy as a function of time.

The energy conservation is considerably worse in the selfconsistent converged-force trajectory ($`n=1`$) than in the Harris-force trajectory ($`n=\mathrm{}`$). This probably reflects larger effects of charge sloshings and analytic discontinuities in the forces, due to frequent level crossings in this highly disordered system. However, it is important to notice that the energy conservation in the mixed-force trajectory ($`n=10`$) is similar to that in the converged trajectory.

Despite the high simulation temperature, the charge transfer in elemental liquid silicon may be expected to be considerably smaller than in an ionic system, making the non-selfconsistent Harris functional specially adequate. In fact, we have seen that structural magnitudes like the bond length and bond angle distributions are not very different using the Kohn-Sham and Harris functionals. Therefore, we have also studied a more challenging system, liquid silica, using 72 atoms at a high average temperature of 5500 K and a low density of 0.42 g/cm3, typical of porous silica aerogels Rahmani et al. (2001). The distributions of bond legths and angles, presented in Figures 5 and 6, are indeed very different using the two functionals ($`n=1`$ and $`n=\mathrm{}`$).

Despite this, the mixed-force method, with up to $`n=10`$, yields essentially the same distributions as the converged Kohn-Sham trajectory. Similar results, to be presented elsewhere, were obtained for an even more ionic system, liquid magnesium oxide, with 54 atoms at 6500 K and 30 GPa.

It might be expected that dynamical magnitudes are more sensitive than thermodynamic averages to changes in how the MD trajectories are obtained. Figure 7 shows the velocity autocorrelation function for the three systems studied, as a function of the interval $`n`$ between force corrections.

As expected, the trajectory of the non-selfconsistent Harris functional is reasonably accurate only for elemental liquid silicon. But, in every case, the mixed-force method, with up to $`n=10`$, yields essentially the same velocity autocorrelations as the converged Kohn-Sham trajectories. We have also calculated self-difussion coefficients from the average quadratic distances traversed as a funtion of time. Thus, for liquid silicon we obtain, respectively, $`(2.4\pm 0.1)`$, $`(2.5\pm 0.1)`$, $`(2.6\pm 0.1)`$, $`(2.6\pm 0.1)`$, and $`(2.0\pm 0.1)\times 10^4`$ $`\text{cm}^2/\text{s}`$, for $`n=1,5,10,20`$ and $`\mathrm{}`$. Again, the mixed-force value, even with $`n=20`$, is the same, within the statistical error, as that of the converged trajectory. This indicates that dynamical and kinetic magnitudes, as well as structural or thermodynamic averages, are well reproduced even with quite large values of the boost factor $`n`$.

In conclusion, we have presented a new method to greatly accelerate ab initio molecular dynamics simulations by combining cheap force evaluations with accurate converged ones. Our results show that the method is very robust with respect to the reduced accuracy of the cheap forces. Although the acceleration factor will undoubtly depend on the system simulated, our present results indicate that factors of 10 can be expected in most cases.

We thank specially Gabriel Fabricius for discussions, for help with the liquid silicon parameterization, and for sharing with us his data-processing programs. We also acknowledge useful discussions with Emilio Artacho. This work has been supported by the Fundación Ramón Areces and by Spain’s MCyT grant BFM2000-1312.


text:
# Eclipse observations of high-frequency oscillations in active region coronal loops

## 1 Introduction

The coronal heating mechanism is the subject of a great deal of debate. With a temperature of more than a million degrees, the corona is several orders of magnitude hotter than the photosphere and chromosphere, thus ruling out the possibility of heating via thermal conduction. Popular theories which attempt to explain coronal heating can be broadly grouped into two categories (see the review article by Priest & Schrijver 1999). One possibility is that a large number of magnetic reconnections, followed by current dissipation, result in frequent micro- or nano-flare activity (Parker 1988). The other theory argues that the heating is dominated by the damping of magnetohydrodynamic (MHD) waves – either propagating from the lower solar atmosphere or induced in active regions by reconnection – through ion viscosity and electrical resistivity (first introduced by Hollweg 1981). MHD waves have two very distinct extremes: magnetoacoustic (divided into slow and fast mode) and ‘pure’ Alfvén (divided into compressional and non-compressional) waves. Magnetoacoustic waves cause pressure variations in the coronal plasma as they propagate. By contrast, Alfvén waves will either be transverse and non-compressional, propagating parallel to the magnetic field, or compressional and modifying the magnetic flux density perpendicular to the field. The differences between the two main categories of waves are expected to be observable, since Alfvén waves cause only Doppler shifts in observed lines, whereas the magnetoacoustic waves are expected to cause intensity variations as well. The latter should be more readily observable since the intensity normally varies with the square of the electron density.

In the past, a number of authors have reported intensity, velocity and line width fluctuations in the corona. Koutchmy et al. (1983) used coronagraph observations to find evidence of fluctuation in the velocity measurements of the coronal line of Fe xiv at 5303 Å with periods of 300s, 80s and 43s. Their observations were limited by sky fluctuations, which led to their suggestion of using satellite telescopes or observations during solar total eclipse. Furthermore Pasachoff & Landman (1984), using observations made during a total eclipse, detected intensity fluctuations with frequencies in the range 0.5–2.0 Hz (periods of 0.5–2.0 s). Since then, more detections of possible MHD oscillations have been published and Aschwanden et al. (1999) produced a catalogue of all detected periodicities across the spectrum from 0.01 to 1000 s. Subsequently, a number of authors have reported further detections of coronal oscillations. Singh et al. (1997) observed fast magnetosonic disturbances with frequencies of 0.2 Hz during the 1995 total solar eclipse. Cowsiket et al. (1999) applied similar techniques to detect oscillations with frequencies in the range of 10-20 mHz during the 1998 total solar eclipse, while Sakurai et al. (2002) used spectroscopic data to detect Doppler velocity with a frequency of $``$57 mHz. More recently, total eclipse observations were published by Pasachoff et al. (2002), who found some evidence for waves with frequencies in the range of 0.75-1.0 Hz.

Porter et al. (1994a, b) have used numerical methods to simulate the damping of energy from slow- and fast-mode MHD waves. They concluded that slow-mode waves can deposit enough energy to heat the corona under certain conditions, and for periodicities $`\tau `$100 s; for fast-mode waves the upper limit to periodicities is $``$1 s. In a bid to detect such high-frequency oscillations, Phillips et al. (2000, hereafter P00) developed the Solar Eclipse Coronal Imaging System (SECIS), an instrument capable of making high-cadence observations of solar eclipses. Williams et al. (2001, 2002, hereafter W01 and W02 respectively) reported SECIS detections of oscillations with frequencies around 0.1 Hz, which might provide enough energy to heat the corona efficiently. Here we report further instances of such high-frequency oscillations in the data as those presented by W01 and W02.

## 2 Data analysis & Results

### 2.1 Observations

Observations of the solar corona in Fe xiv 5303 Å (formed at $``$2.0 $`\times `$106 K during the August 1999 eclipse were taken using the SECIS instrument (for a detailed description of the instrument see P00 and W01), with a sampling rate of 44 frames per second. The field of view of the instrument was $`0.5\times 0.5`$ deg2 with a pixel size 4.0 arcsec pixel-1. To increase the signal-to-noise, for each pixel the intensities of the eight adjacent pixels were added to it; i.e. as in W01 and W02 we have summed over a 3$`\times `$3 pixel2 area. Two active region loops in NOAA AR 8651, observed both by the Solar and Heliospheric Observatory (SoHO) and through an Fe xiv (5303 Å) filter by SECIS, are highlighted in Fig. 1. Both loops were chosen as they are among the most well isolated in this active region and can clearly be seen in data from both SECIS and the EUV Imaging Telescope (EIT) on board SoHO. The latter observations were used to confirm the positions of these loops using several reference features.

### 2.2 Wavelet analysis

As in W01 and W02, we chose to analyse the data using a continuous wavelet analysis (for further details on this technique see Torrence & Compo 1998, hereafter TC98). Although Fourier analysis is overall more widely used, wavelet analysis has recently gained popularity, due to its ability to detect oscillations localised both in time and frequency. If an oscillation only lasts for a small fraction of the time series duration, Fourier analysis will be unable to detect it, whereas wavelet analysis is sensitive to even transient oscillatory signals. As the oscillations we hope to detect are relatively short (of the order of a few seconds) and our data extend over a period of around 40 s, wavelet analysis is ideal. Other authors (e.g. Gallagher et al. 1999; Ireland et al. 1999; Banerjee et al. 2000) have applied the same technique to detecting solar coronal oscillations. Additionally, the wavelet technique has consistently been used by W01 and W02, and a comparison with the results from these publications may help to draw interesting conclusions.

A Morlet wavelet was used for the analysis of our data, with

$$\psi (\eta )=\pi ^{1/4}\mathrm{exp}(i\omega _0\eta )\mathrm{exp}(\frac{\eta ^2}{2}),$$
(1)
where $`\eta =t/s`$ is the dimensionless time parameter, $`t`$ is the time, $`s`$ the scale of the wavelet (i.e. its duration), $`\omega _0=s\omega `$ is the dimensionless frequency parameter, and $`\pi ^{1/4}`$ is a normalization term (see TC98).

### 2.3 Detections

Fig. 2 contains a time integrated (40 s) SECIS image of approximately the same region as Fig. 1. Since the lower parts of the corona are covered by the Moon, only the apexes of the loops identified in EIT Fe xii (195 Å) are observed by SECIS. The solid lines indicate the positions of the loops marked in Fig. 1 and the pixels marked with squares the areas that show intensity oscillation lasting three periods or longer. All twenty regions containing detected oscillations are labelled from A to T. Importantly, all the detections were made outside both bright loops toward the side where the corona is more tenuous. Although the areas of the corona that appear to host intensity oscillations are outside the visible loops, this has not led us to believe that these waves perturb outside coronal loops. It could simply be that the loops they travel through are low emission-measure structures that are too faint to be visible above the background. Furthermore, one may observe that in the case of the left-hand loop, the detections exactly coincide with a shift of the line which highlights the loop by one pixel to the left and two pixels down. The dashed line of Fig. 2 illustrates the close correlation between the shifted line of the loop and the pixels that have been detected with long intensity oscillations. Also, phase analysis of points A–T reveals no correlation between the phase of the oscillations detected in each point in space or time. This leads us to believe that this is not a simple case of a standing or travelling wave.

Figs. 3–6 contain the results of the wavelet transform of a sample of points A to T in Fig. 2. The remaining figures are available from the first author upon request. Each figure is divided into three sections. Part (a) shows the time series generated by the Fe xiv (5303 Å) line filter (for details about the filter and the data reduction see W01). Part (b) is the power density wavelet transform of the time series in (a) : the brighter an area, the greater the power at the given time and frequency. The vertical axis is logarithmically scaled in frequency space while the horizontal axis is time on a linear scale. The time axis in (b) exactly coincides with that of the time series in (a), with the hatched region of (b) marking the cone-of-influence (COI). Everything inside the COI should be treated with suspicion, since any detections in this area may be influenced by edge effects in the wavelet transform calculation. For example a detection of a 0.5 Hz oscillation 0.5 s before the end of the time series is unreliable as there is not sufficient time for the oscillation in the wavelet packet to finish. For a more detailed discussion of the problem see TC98. The contours of panel (b) indicate where detected power exceeds the 99% confidence level, i.e., there is a 1% chance of the detection being due to Poisson noise.

Panel (c) contains the global wavelet spectrum, which is the wavelet analogue of the standard Fourier transform. It is produced by summing the power density wavelet transform over the whole time series, while the dotted line running along the frequency axis is the global significance level (again summed over time) at the same value (99%) as the contours in panel (b). The horizontal dot-dashed line near the bottom of panel (c) marks the bottom of the COI and all detections below this frequency should be discarded. Since the third panel provides no time information for the detections, the rest of the COI cannot be defined.

At this point it should be emphasized that there are two additional basic differences between panel (c) and the ‘traditional’ Fourier transform. Firstly, the spectrum is much smoother than the discrete Fourier power spectrum. Secondly, due to the wavelet transform’s use of ‘Heisenberg boxes’ (Mallat 1998), the better defined the transform is in time, the less well defined it will be in frequency.

Table 1 contains all the frequencies detected in each of the points from A to T. The length of these oscillations is also included in units of periodicities (i.e., a duration of three means that this particular oscillation lasted for three oscillatory periods at this frequency).

## 3 Discussion

Although W01 and W02 have published similar detections from SECIS observations, the present work is by far the largest number of detections published to date. All twenty points presented here (Table 1) have passed a number of selection criteria. These include those used by W01 and W02 and the most important are summarised below:

* The frequencies of the detections are distinct from known instrumental frequencies (see W01).
* The contours of panel (b) were chosen at a 99% confidence level. Only oscillations within those contours were considered.
* All the reported detections lasted for at least three periods, so as to rule out rapid increases or decreases in the signal. When the duration of the oscillations was calculated, any portion within the COI was discarded.

As with all high-cadence systems producing data analysed for oscillations, the introduction of instrumental frequencies and the effect of noise is always a concern. To address the first limitation, wavelet analysis was deployed for parts of the image covered by the moon’s disk and the very faint parts of the corona. Any non-localised instrumental variations in pixel intensity should have affected these areas as much as the pixels of the active region. With the exception of those instrumental frequencies that are already known and discussed by W01, no other frequencies were detected. The other well-know cause of false detections is noise. To limit the possibility of a false detection because of noise, we chose to ignore any detections with frequencies above 1 Hz. It is widely accepted (for example Starck & Murtagh (2002) and references therein), that Gaussian or Poisson noise only affects frequencies of the same order as the sampling rate of the time series. As the sampling rate of the SECIS 1999 observations was 44 frames per second, the reported frequencies should be fairly unaffected by Gaussian or Poisson noise (which includes types of noise such as the CCD readout).

The two coronal loops of AR 8651 analysed here show a significant number of oscillations. Although only twenty oscillations were detected lasting three periods or longer, several tens of wave signatures were found which last between two and three periods at a $`>99`$% confidence level. This number is much larger than that expected considering the previous work on SECIS data (W01 and W02). Moreover, all previous detections were from the interior of bright coronal loops, while all of our oscillations are detected in fainter loops within the same active region, toward the tenuous part of the corona. Both of the above peculiarities can be explained by introducing a different physical mechanism to that suggested in W02. Zaqarashvili & Roberts (2002, 2003) have suggested a swing wave-wave interaction mechanism which may cause energy transformation from fast magnetoacoustic waves propagating across a magnetic field to Alfvén waves propagating along the field. They argue that, for a given medium density and magnetic field, the energy of fast waves can be converted to Alfvén waves with a basic harmonic at half the wavelength of the fast-mode wave. The above mechanism (also called swing absorption) may provide a possible explanation for the intensity oscillations reported here. In this scenario, Alfvén waves created in the upper photosphere (as described by Zaqarashvili & Roberts 2002, 2003) propagate along magnetic field lines adjacent to the coronal loops in Figs. 1 and 2. In the case of the first loop (points A, B, C, J, K, L, M, N of Fig. 2), a magnetic field line was running almost parallel to the left bright loop of Fig. 1, producing the alignment highlighted in Fig. 2. The rest of the points with detected oscillations (D–I and O–T) belong to one or another of the magnetic field lines in the same active region.

Litwin & Rosner (1993) proposed a multi-thread model, where many tiny loops with different physical parameters but in steady state equilibrium are superimposed forming the observed coronal loop structures. Aschwanden et al. (2000) used this model to explain the nonuniform heating of coronal loops observed by the Transition Region And Coronal Explorer (TRACE). Several such small threads may form thinner and fainter loops outside the bright structures appearing on Fig. 1 and 2. As the above oscillations travel through low emission-measure loops, it is easier for us to detect more oscillations than through high emission-measure loops (such as those studied by W01 and W02). Therefore the conditions under which those MHD waves propagate along points A–T are significantly different from those of the oscillations reported by W01 and W02. In their case the emission from the propagating waves is much stronger and as they and travel through a large number of threads, they stand significantly above the background emission of the loop. This is more clearly seen in W02, where the event that caused the oscillations also cause them to propagate with the same phase, enabling W02 to calculate the velocity of the perturbation by the phase difference of the traveling wave across the loop. In contrast to those detections, the area outlined by points A–T contains a smaller amount of threads with lower emission-measure, therefore relatively weaker MHD oscillations can be detected. This is supported by the results of the phase analysis that reveal no correlation in the phase of the 20 points with detected oscillations. The mechanism that produced the weaker waves (compared to those reported by W01 and W02) is less likely to produce them with the same phase. As the single-phase MHD waves are relatively rare events in the solar corona, we would expect them to appear in more extreme conditions (such as the release of large amounts of energy at the foot points of the loops) while the propagation of weaker oscillations through a smaller number of threads is more likely to take place under phase mixing conditions.

Cooper et al. (2003) suggest a possible mechanism that explains in some detail the detection of intensity oscillations as a line-of-sight effect of entirely incompressible MHD waves. In this model, when observed at an angle $`\theta `$ to the direction of propagation, the wave-induced deformation in a coronal loop causes intensity variations. This is because the amount of optically thin emitting plasma along the line of sight changes as a function of time. Cooper et al. (2003) find that the observed amplitude of the intensity oscillation can vary as a function not only of the true intensity of the oscillation and the angle between the propagation direction and the line-of-sight, but also of the wavelength of the perturbation. Furthermore, the observed frequency also varies as a function of the angle $`\theta `$, meaning that the detected periods listed in Table 1 may simply be higher harmonics of the true values.

Using the Cooper et al. (2003) and Zaqarashvili & Roberts (2002) results we were able to provide a satisfactory explanation of how the detected incompressible MHD waves were created in the vicinity of an active region in the photosphere, transmitted through low emission-measure loops to the lower corona and then detected as intensity oscillations by our imaging system. The large number (comparing to previous work) of detections and the alignment of some of these can be explained as the low emission that comes from the tenuous plasma makes any intensity oscillations more apparent, since the oscillating material makes for a higher percentage of the detected intensity.

###### Acknowledgements.

The authors would like to thank K.J.H. Phillips for his collaboration on the SECIS project. ACK acknowledges funding by the Leverhume Trust via grant F00203/A. DRW acknowledges a CAST studentship funded by the Department of Employment & Learning and the Rutherford Appleton Laboratory. JMCA acknowledges CAST studentship funded by DEL and QUB. ACK & DRW would like to thank Valery Nakariakov and Temury Zaqarashvili for useful discussions.


text:
# Effect of Measurement on the Periodicity of the Coulomb Staircase of a Superconducting Box

(5 November 2003)

## Abstract

We report on the effect of the back-action of a Single Cooper Pair Transistor electrometer (E) on the measurement of charge on the island of a superconducting box (B). The charge is e-periodic in the gate bias of B when E is operated near voltages $`2\mathrm{\Delta }/e`$ or $`4\mathrm{\Delta }/e`$. We show that this is due to quasiparticle poisoning of B at a rate proportional to the number of quasiparticle tunneling events in E per second. We are able to eliminate this back action and recover 2e charge periodicity using a new measurement method based on switching current modulation of E.

The superconducting box, e.g. a small island of Al film very weakly coupled to the outside circuitry by Josephson junctions, has shown considerable promise as a qubit for quantum information processing where the two states can be represented by superpositions of 0 or 1 excess Cooper pairs in the box Shnirman ; Averin ; Nakamura ; Vion . Measurement of the quantum state of this so-called charge qubit without inducing unwanted decoherence is a significant problem as is quasiparticle poisoning, i.e. the introduction of an unpaired electron (quasiparticle) into the box. At temperatures of 10 mK or so, where experiments are commonly done, the number of quasiparticles should, in principle, be negligible. However, such quasiparticle poisoning, due perhaps to the measurement process itself, is commonly observed. A manifestation of this is seen in the so called Coulomb staircase. When a charge, $`q_{g,B}`$, is capacitively induced on the box, one expects Cooper pairs to tunnel resonantly into or out of the box at $`q_{g,B}=n_oe`$, where $`n_o`$ is an odd integer, to maintain the lowest energy charge state of the box. This results in the Coulomb staircase of the charge in the box $`Q_B(q_{g,B})`$ with period $`2e`$ in $`q_{g,B}`$. On the other hand, if there are quasiparticles in the system, then maintaining the lowest charging energy state also leads to quasiparticle tunneling. This gives rise to splitting of the steps in the Coulomb staircase, which shifts toward e-periodicity as the number of quasiparticles increases Lafarge . As a result, the lowest energy state of the box at $`q_{g,B}=n_oe`$ no longer corresponds to a resonant state of the Cooper pair tunneling. For the box qubit, this means that relaxation does not bring the system back to its computational ground state at its operating point. Since the ability to prepare the initial state of the qubit is an absolutely necessary condition for quantum computing, quasiparticle poisoning has been a serious road block for groups working to build a charge qubit Flees ; vanderWal ; Lehnert . Solving this problem in charge qubits is essential.

The purpose of this paper is to investigate the effects of measurement (i.e. back-action) on the measured charge in the box and to develop approaches to minimize these effects. For this study, we use two capacitively coupled Single Cooper Pair Transistors (SCPT)s, one of which acts as an electrometer (E) and the other as a superconducting box (B). The latter gives a good representation of the box and at the same time allows us to study quasiparticle poisoning effects without needing to operate E. The SCPT electrometer can be operated in several different modes when measuring the charge of B Schoelkopf ; Zorin . Commonly, the charge measurements of B are done by operating E in the voltage modulation mode (VM). In this mode, E is biased at a sufficiently high voltage that quasiparticles are generated, so it effectively functions as a SET where the source drain voltage is modulated by $`q_{g,E}`$ with a period of $`e`$. However, it is known that the switching current of a SCPT, i.e. the current at which it switches hystereticaly from the low voltage or phase-diffusion branch to $`eV>2\mathrm{\Delta }`$, is also charge sensitive Flees  and can be used for charge measurement. We refer to this as the switching current mode (SW). The SW mode of operation has been analyzed Cottet , but until now no measurements of the charge on the island of a box using a SW mode of the electrometer have been reported. We present the results of the measurements of the island charge in B by E operated in either the VM or SW mode. The results demonstrate that measurement-induced poisoning, which leads to an e-periodic Coulomb staircase when using E in the VM mode, can be eliminated in the SW mode.

The parameters of the sample (Fig. 1.a) are as follows: E has a normal state resistance $`R_{n,E}=61.4`$ k$`\mathrm{\Omega }`$, charging energy $`E_{c,E}=33`$ $`\mu `$eV and Josephson energy $`E_{J,E}=21.2`$ $`\mu `$eV. For B these parameters are $`R_{n,B}=63.0`$ k$`\mathrm{\Omega }`$, $`E_{c,B}=47`$ $`\mu `$eV and $`E_{J,B}=20.6`$ $`\mu `$eV. $`E_J`$ is the average of the Josephson energies of two junctions as determined from the values of $`R_n`$ and the superconducting gap $`\mathrm{\Delta }200`$ $`\mu `$eV by using the Ambegaokar-Baratoff formula. $`E_c`$ is determined from the amplitude of maximum voltage modulation of the devices. The coupling capacitance between E and B is determined to be 80 aF from their measured coupling of 4.8%. All devices are made using standard two angle shadow evaporations without having normal metal quasiparticle traps close to junctions. The sample is placed in microwave tight copper can located on a temperature regulated stage of a dilution refrigerator having a base temperature of 6 mK. All the measurement leads are filtered by low temperature microwave filters Danthesis  that are thermally anchored to the mixing chamber.

Figure 1.b shows the charge on the island of B measured with E in the VM mode. During this measurement, the source and drain leads of B are at a common potential with respect to its gate. Figure 1.c presents the switching current modulation of B measured with the bias current through E, $`I_E`$, set equal to zero. In this and following measurements, the bias current of B, $`I_B`$, is ramped at a rate of $`140`$ nA/s. As can be seen, the switching current modulation of B is 2e-periodic as expected at low temperature, but the charge of B, measured by the electrometer, is e-periodic. Similar dependences of $`Q_B`$ on $`q_{g,B}`$ were measured with E biased in either of its voltage sensitive regions, i.e. near the gap where $`V_E4\mathrm{\Delta }/e`$ or near the Josephson-quasiparticle peak where $`V_E2\mathrm{\Delta }/e`$. To determine if the e-periodicity of $`Q_B(q_{q,B})`$ is due to the back-action of E on B, the quasiparticle poisoning rate of B, $`\gamma _B`$, was measured for a range of bias conditions of E. In addition to this, we studied how the biases of two other SCPTs located on the same chip, but coupled more weakly to B, affected $`\gamma _B`$.

The details of the technique for determining $`\gamma _B`$ from switching current distributions have been reported previously Mannik . Briefly, as $`I_B`$ is linearly ramped in time, the switching current histogram of B, with its gate biased near $`q_{q,B}=n_oe`$, exhibits two peaks if the number of quasiparticles on the island changes during the measurement of histogram. One peak, $`I_{even}`$, which is close to the maximum of the switching-current characteristic $`I_{sw,B}(q_{g,B})`$, occurs when the island has an even number of electrons (even state) and the other -much lower current- peak, $`I_{odd}`$, is near the predicted switching current minimum at $`q_{g,B}=0`$ when one quasiparticle occupies the island (odd state) (Fig. 2). If $`I_B>I_{odd}`$ and B has not switched, it must be in the even state. The entry of a quasiparticle onto the island effectively changes $`q_{g,B}=e`$ to $`q_{g,B}=0`$ which for $`I_B>I_{odd}`$ will cause B to switch rapidly to the running state, giving the time of the poisoning and thus the quasiparticle poisoning rate $`\gamma _B`$. Previous studies have shown that $`\gamma _B`$ is independent of $`I_B`$ in the region between the peaks, giving an exponential decay of the even state. Several of these histograms with increasing $`\gamma _B`$ as $`V_E`$ and $`I_E`$ increase are shown in Fig. 2. The relatively low bandwidth of our filters limits these measurement to $`\gamma _B<10`$ ms-1.

First we determine $`\gamma _B`$ when all the electrodes of the electrometer are disconnected from the measurement circuitry and grounded. In this case we still observe a small residual rate $`\gamma _B^0=0.06`$ ms-1 for $`q_{g,B}n_oe`$. This non-zero rate can be caused, e.g., by the presence of impurity levels in the superconducting gap of Al Lafarge . For the present discussion it is clear that this small residual rate is not related to the back-action of E. $`\gamma _B`$ is unchanged if E is biased on its supercurrent branch or when E is biased on its return current branch at low voltage $`V_E<200`$ $`\mu `$V $`\mathrm{\Delta }/e`$ as shown in Fig. 3. For $`V_E>200`$ $`\mu `$V, $`\gamma _B`$ decreases rapidly, becoming too short to measure for slightly higher voltages. The modulation characteristics of the current $`I_E`$ with gate voltage also changes at this point from being 2e-periodic for $`V_E<200`$ $`\mu `$V to e-periodic for $`V_E>200`$ $`\mu `$V. Similar cross-overs from e to 2e periodicity at a voltage $`\mathrm{\Delta }/e`$ have previously been seen in related systems (see e.g. Bibow ).

In order to study how $`\gamma _B`$ depends on the voltage and current of E through its entire operating range, and in particular near $`V_E=2\mathrm{\Delta }/e`$ and $`V_E=4\mathrm{\Delta }/e`$, we use two other SCPTs fabricated on the same chip but much more weakly coupled to B. Since these two SCPTs, which we call E1 and E2, are much more weakly coupled to B, they allow us to measure $`\gamma _B`$ for bias currents, $`I_{E1,2}`$, that are several orders of magnitude higher than is possible using E. The island of SCPT E1 is located 96 $`\mu `$m from the island of B and has $`R_n=61.4`$ k$`\mathrm{\Omega }`$, while the corresponding parameters for E2 are 143 $`\mu `$m and $`R_n=54.2`$ k$`\mathrm{\Omega }`$. For these devices, we can measure the rate $`\gamma _B`$ up to voltages $`4\mathrm{\Delta }/e`$. Again, we see a small initial increase of $`\gamma _B`$ at $`V_E200`$ $`\mu `$V and then sharp increases at voltages $`V_E2\mathrm{\Delta }/e`$ and $`4\mathrm{\Delta }/e`$. These voltages correspond approximately to the Josephson-quasiparticle tunneling and sequential quasiparticle tunneling thresholds in a SCPT and are accompanied by sharp increases in $`I_{E1,2}`$. The inset of Fig. 3 shows the rate $`\gamma _B`$ as a function of the currents in E1 and E2 above these two thresholds. To test the hypothesis that $`\gamma _B`$ is proportional to the total number of quasiparticle tunneling events, the currents for the electrometer voltages $`V_E>4\mathrm{\Delta }/e`$ are multiplied by two. This is done since, for a given current at voltages $`V_E>4\mathrm{\Delta }/e`$, there are twice as many quasiparticle tunneling events through the junctions of E as for $`V_E2\mathrm{\Delta }/e`$. As one can see, this scaling collapses all of the data from each electrometer to a common line.

From these data one can conclude that the quasiparticle current of electrometer is the source of back-action noise leading to quasiparticle poisoning and an e-periodic Coulomb staircase of B when measured by E. Further, the quasiparticle generation in B is proportional to the total number of quasiparticle tunneling events per second through E. This relationship could indicate that the back-action of E results from the shot noise of tunneling quasiparticles. On the other hand, this back-action could also be the result of the recombination of quasiparticles in E into pairs. This quasiparticle recombination produces phonons and to a smaller extent photons of energy $`2\mathrm{\Delta }`$. These phonons/photons, which propagate from E to B without energy relaxation could generate quasiparticles in B. Determining the details of the interaction between E and B will require further work. However, it is interesting to note that this sort of recombination noise would likely be suppressed by having normal metal leads close to the junctions. This may provide an explanation of previous results Bouchiat  in which it was possible to observe a 2e-periodic Coulomb staircase using a VM electrometer measuring a superconducting box, which had normal metal “quasiparticle traps” close to the junctions.

Quasiparticle poisoning of B due to the measurement of its charge can, in principle, be eliminated by the operation of E in the SW mode, where its voltage remains well below $`\mathrm{\Delta }/e`$ until the measurement is made. The SW mode of operation of the electrometer is illustrated in the inset of Fig. 4, which shows its switching current vs. gate charge transfer-function. In general, $`q_{g,E}`$ contains a component proportional to $`Q_B`$, in addition to the externally applied bias. So, to measure changes in $`Q_B`$ as a function of $`q_{g,B}`$, we bias E near point “0”, where the switching current is very sensitive to variations of external charge ( inset of Fig. 4), and record several hundred switching events of E for each value of $`q_{g,B}`$. These switching data are then corrected for the measured nonlinearity in the transfer-function of E and averaged. The top curve in Fig. 4. shows the result obtained using this procedure.

While these data are no longer strictly e-periodic, as in Fig. 1b, the Coulomb staircase in Fig. 4 (top line) still shows the split steps around $`q_{g,B}=n_oe`$ characteristic of quasiparticle poisoning. We can see that this is consistent with the measured residual rate, $`\gamma _B^0`$, as follows. The edges of these steps should occur for an energy difference between the even and odd states, $`\mathrm{\Delta }E(q_{g,B})`$, such that the rates for even to odd ($`\gamma _B^0`$) and odd to even ($`\overline{\gamma _B^0}`$) transitions of the island are equal. This ratio is given by $`\gamma _B^0/\overline{\gamma _B^0}=D_e/D_o\mathrm{exp}(\mathrm{\Delta }E(q_{g,B})/kT)`$ Schon , where $`D_e`$ is the sum of the quasiparticle densities on the leads and the island in the even state. For the odd state, this sum is $`D_o=D_e+V_i^1`$, where $`V_i`$ is the volume of the island. The second term in $`D_o`$ accounts for the extra quasiparticle occupying the island in the odd state, which increases the density by $`V_i^1100`$ $`\mu `$m-3. $`D_e`$ can be estimated from the measured residual poisoning rate $`\gamma _B^0`$, the normal density of states of the Al film and the junction resistances, giving $`D_e410^3`$ $`\mu `$m-3. Taking the electron temperature of B to be 15 mK, which is reasonable since B is completely passive in these measurements, gives a length for the short step of $`0.76e`$ in agreement with the data shown in Fig. 4.

The effects which result from the residual rate $`\gamma _B^0`$ can be greatly reduced by flushing the quasiparticle from the island of B before each measurement. As one possible approach to prepare the even parity state, we apply a voltage pulse $`V_B`$ across B just prior to each measurement. The amplitude of $`V_B`$ is chosen such that $`2E_{c,B}<V_Be<\mathrm{\Delta }`$, in order to release the quasiparticle from the electrostatic potential of the island but yet not to generate any new quasiparticles by the pulse explanation . Switching histograms of B show that this procedure prepares the even state with a probability of about 85%. Immediately after B is flushed, the measurement ramp of $`I_E`$ begins. The result of this procedure is shown in the bottom curve in Fig. 4. As one can see, the quasiparticle-induced splitting of the step at $`q_{q,B}=n_oe`$ is no longer apparent. However, individual histograms still show about 30% of the switching events in E near $`q_{q,B}=n_oe`$ are from the poisoned state of B. This is consistent with the imperfect preparation of the initial state and additional poisoning with rate $`\gamma _B^0`$ in the finite time between flushing and the switching event. Thus we see that, with the effects of residual poisoning greatly reduced by the flushing, the measurement of $`Q_B`$ by E in the SW mode gives results consistent with the 2e-periodic switching-current distribution of B.

In conclusion, our measurements clearly show that operation of a SCPT electrometer at voltages $`V_E>200`$ $`\mu `$V causes a substantial generation of quasiparticles in the circuit of the superconducting box leading to an e-periodic Coulomb staircase. The rate of quasiparticle poisoning in the box depends linearly on the total number of quasiparticle tunneling events per second through the junctions of the electrometer. To overcome this back-action from the electrometer, we operate it in a mode which uses switching-current modulation for charge detection. Using this mode of operation, we are able to recover the 2e-periodic Coulomb staircase of the SCPT in the box configuration, which is expected both theoretically and from the 2e-periodicity of its switching current.

The authors thank D. V. Averin, J. R. Friedman and K. K. Likharev for useful discussions and W. Chen and V. V. Kuznetsov for technical assistance. Work is supported in part by AFOSR grant No. F49620010001.