text:
# n-Alkyl Thiol Head Group Interactions with the Au(111) Surface

## Abstract

State-of-the-art first principles calculations based on density functional theory were performed on CH3(CH2)n-1S-Au(111) systems. We show that the adsorption site of methylthiolate at low coverage on the Au(111) surface is the fcc site, not the hcp site as has been recently reported. Further, we report results for chain length dependency and the electronic structure of the system.

The interaction between organic materials and solid surfaces has been extensively studied because of the broad range of industrial applications. Self-assembled monolayers (SAMs) hold special interest, because the presence of the thiol group greatly strengthens the molecule-surface interactions, inducing order in the layer. SAMs have important potential applications in industry, such as sensors, transducers, detectors, packaging and insulating layers for integrated circuits, functionalization of surfaces, thin coatings for electrodes, and corrosion inhibition.

The long-chain alkane thiols \[CH3(CH2)n-1SH, or Cn\] form SAMs on the Au(111) surface. Their simplicity, highly ordered structures, and chemical stability make these systems ideal for study with a variety of techniques including atomic force microscopy, infrared spectroscopy, high-resolution electron-energy-loss spectroscopy, grazing X-ray diffraction, scanning tunneling microscopy (STM), scanning probe microscopy, low-energy electron diffraction (LEED), He atom diffraction, and theory. Despite extensive studies of this system, their are many controversies regarding its structure.

Two competing structural models, the “standard model” and the “sulfur-pairing model”, have considerable support. In the standard model, alkanethiol molecules exhibit a hexagonal ($`\sqrt{3}\times \sqrt{3}`$)R30 lattice. Alkanethiolates occupy three-fold hollow sites of the Au (111) surface with S-S spacings of 4.99 Å, and the molecular axes are tilted by 30-35 with respect to the surface normal. The sulfur-pairing model involves a $`c`$(4 $`\times `$ 2) superlattice of the hexagonal lattice, where alkanethiolates form sulfur head group dimers with S-S spacings of 2.2 Å. Furthermore, other structural models have been reported, including $`n\times \sqrt{3}`$ unit cells ($`8n10`$). There is even controversy among proponents of the sulfur pairing model, with some researchers suggesting top and hollow-bridge positions for the sulfurs, while others propose a hollow and hollow-bridge dimer.

In order to understand the structure fully one must include the effect of coverage, chain length, and temperature. However, because the differences between these structures center around the thiol head-group locations, theoretical investigation of the thiol-gold chemisorption bond is vital. Since interactions with the metal surface appear to dominate the molecular arrangement of SAM systems, the interaction of a methylthiolate (CH3S) with the Au (111) surface presents a prototype system to study this entire class of systems.

Both LEED and ultrahigh-vacuum cryogenic STM studies show a ($`\sqrt{3}\times \sqrt{3}`$)R30 pattern of CH3S monolayers at room temperature, with head-groups at three-fold hollow sites, while at 110 K a disordered hexagonal ($`3\times 4`$) lattice has been reported. The S-S distances in the latter study (3.3 Å) are much longer than the 2.2 Å and much shorter than 4.99 Å from the sulfur-pairing and standard models. Two theoretical studies have concluded that the adsorption site for a single thiol molecule on the gold (111) surface is the hcp hollow site. These results also support the standard model. However, despite the overall agreement of the results, there are discrepancies regarding the S-Au distance and the tilt of the molecule in these two independent studies.

The main objective of our study is to identify the adsorption site for thiols on gold (111), and to explain the structural properties. We have performed first-principles calculations based on density-functional theory (DFT) to study the adsorption of methylthiolate on the Au (111) surface. Our results show that a single CH3S molecule adsorbed on the Au (111) surface chemisorbs on the fcc site. We also find that further increase of the hydrocarbon chain does not affect the adsorption site. All of the first-principles calculations have been performed with the plane-wave pseudopotential code “dacapo” .

In our DFT calculations, the wave functions are expanded in a plane-wave basis set, the electron-ion interactions are described by ultra-soft pseudopotentials (USPP) , and the generalized gradient approximation (GGA)  for the exchange-correlation functional has been used. The Kohn-Sham equations are solved self-consistently, using a Pulay density-mixing scheme  to update the electronic density between iterations. The occupation numbers are updated using a recently developed technique based on minimization of the free-energy functional . A finite electronic temperature is used, in order to reduce the number of $`k`$ points needed, and all total energies are then extrapolated to zero electronic temperature. For the optimization of atomic structure, a damped molecular-dynamics (DMD) method has been used.

When using the slab geometry in the plane wave pseudopotential method, errors most often arise from the pseudopotential or from lack of convergence with respect to the number of $`k`$-points, the cutoff energy, and the system size. When generating pseudopotentials it is of great importance to preserve the eigenvalues for all relevant atomic configurations, not just in the reference configuration. However, in order to have correct lattice constants for metals it is also of importance to preserve the charge density in the tail region for all relevant atomic configurations . For H (1s1), C (2s22p2), S (3s23p3.53d0.5), and Au (5d9.56s16p0.5) the USPPs were generated using cutoff radii (rc) of 0.60, 1.24, 1.45, and 2.00 bohr, respectively, and the maximum transferability errors were less than 2 mRy. All potentials are tested and compared to experiment and to all-electron calculation.

We find that although the $`3d`$ orbitals of S are nearly unpopulated in simple molecules, such as S2, the $`d`$ channel has great impact on the behavior of the pseudopotential. The 3$`d`$ orbitals give S enhanced polarizability and enable a great variety of bonding configurations through hybridization. A comparison between two USPPs for S with and without the $`d`$-projector (transferability errors less than 2 mRy in both potentials) show binding energies (Eb of S2) of 2.51 and 3.08 eV/atom (all-electron and experimental values are 2.49 eV/atom  and 2.37 eV/atom ), respectively.

In order to reduce calculation errors due to finite slab thickness, different slab geometries (3-8 layers thick) have been tested using different numbers of $`k`$-points and vacuum layers (see Fig. 1). A slab consisting of 6+7 Au+vacuum layers, a cutoff energy of 30 Ry, and a 4$`\times `$4$`\times `$1 grid of special $`k`$-points is found to give converged results (convergence within a few meV). Table I shows the equilibrium parameters for the clean Au (111) surface and the free CH3S and CH3SH molecules compared to experiment and other calculations.

The first step toward a detailed understanding of the SAM structure is to investigate the adsorption site of the simplest Cn (CH3S) on the surface. Therefore, we calculated the total energy of a single CH3S molecule in a (2$`\times `$2) surface unit cell, for a coverage of 0.25. As shown in Fig. 3a, the thiolate was moved from top to fcc, then over the bridge site to the hcp site using different molecular orientations. The supercell consisted of 6 Au(111) layers and 7 layers of vacuum. The three topmost Au layers and the molecule were optimized using the DMD method until the total force for the system was less than 0.01 eV/Å. Figure 4 shows calculated potential energy surfaces for two different molecular orientations along the diffusion path. The dashed curve is for the orientation used in Ref.  and the solid line is the minimum-energy path. Our results show that the staggered configuration is preferred over the configuration used in previous studies. The results of surface buckling and other optimization effects for the preferred orientation (solid line) are summarized in Table II and shown schematically in Fig. 3b.

The calculated chemisorption energy difference between fcc and hcp, $`\mathrm{\Delta }`$Ehcp-fcc = 0.10 eV, shows that a single CH3S molecule adsorbed on the gold surface prefers the fcc site over hcp. This is in contrast with results of Refs.  and , which are based on cluster calculations and classical MD. In the cluster calculations, only two gold layers were used. As a result, the surface energy is far from converged, and the effect of gold atoms in the third layer is not included. As can be seen from Fig. 1, the gold (111) surface energy is converged only when six or more gold layers are used. Furthermore, the energy difference between having the molecule at the fcc or hcp sites, ($`\mathrm{\Delta }`$Ehcp-fcc) increases from 0.04 to 0.10 eV as the slab thickness is increased from four to six gold layers. This indicates the importance of slab size and suggests that the use of a two-layer gold slab is insufficient for accurate results.

Regarding the effect of hydrocarbon chain length, we performed calculations using two and three carbons in the chain. The results show no changes in the preferred adsorption site (fcc), and $`\mathrm{\Delta }`$Ehcp-fcc is 0.15 eV and 0.18 eV for two- and three-carbon thiols, respectively. These calculations indicate that the preference for fcc becomes stronger with increasing number of carbons in the chain.

The induced charge density (Fig. 4) shows an increase of charge between the Au and S atoms, when a thiol molecule is adsorbed on the fcc site. The metallic electrons (via Pauli repulsion) shift this region of enhancement, so that it is farther from the surface than expected.

In summary, we have presented state of the art DFT calculations of thiol molecules adsorbed on the Au (111) surface at low coverage. This study demonstrates that the fcc site is the preferred location for a single thiol, and the electronic changes which accompany chemisorption are elucidated. Furthermore, the importance of having a minimum of six gold layers in the model, and of including $`d`$ orbitals in the sulfur pseudopotential were highlighted.

The authors wish to acknowledge G. Scoles for his insightful comments concerning this subject. This work was supported by NSF grant DMR 97-02514 and the Air Force Office of Scientific Research, Air Force Materiel Command, USAF, under grant number F49620-00-1-0170. AMR would like to thank the Alfred P. Sloan Foundation for support. Computational support was provided by the National Center for Supercomputing Applications and the San Diego Supercomputer Center.


format: text

text:
# On the Dynamics of the Evolution of HIV Infection

## Abstract

We use a cellular automata model to study the evolution of HIV infection and the onset of AIDS. The model takes into account the global features of the immune response to any pathogen, the fast mutation rate of the HIV and a fair amount of spatial localization. Our results reproduce quite well the three-phase pattern observed in T cell and virus counts of infected patients, namely, the primary response, the clinical latency period and the onset of AIDS. We have also found that the infected cells may organize themselves into special spatial structures since the primary infection, leading to a decrease on the concentration of uninfected cells. Our results suggest that these cell aggregations, which can be associated to syncytia, leads to AIDS.

The human immunodeficiency virus (HIV), which causes AIDS (the acquired immunodeficiency syndrome) has been the subject of most intense studies, that encompass diverse fields of scientific research, ranging from the natural to the social sciences. Major progress has been achieved by medical and biological researchers in understanding the genetic code of the virus, the role of its high mutation rate, the virus-host interaction, the apparent failure of the immune system to control and eliminate the virus, and the onset of AIDS. Nevertheless, the mechanisms by which HIV causes AIDS still remain unexplained.

The immune response to any virus pathogen (virus, bacteria, etc.) is generated by a complex web of interactions, involving the cooperative and collective behavior of different types of white blood cells (monocytes, T- and B-cells). The time scale involved to develop a specific immune response may vary from days to weeks. The complex dynamic of HIV infection and the ensuing onset of AIDS involves a wide range of time scales . The primary infection exhibits the same characteristics as any other viral infection: a dramatic increase of the virus population during the first $`26`$ weeks, followed by a sharp decline, due to the action of the immune system. For HIV, however, instead of being completely eliminated after the primary infection, a low virus burden is detected for a long asymptomatic time: the clinical latency period. This period may vary from two to ten (or more) years, depending on the patient. Besides the low concentration of HIV detected during this period, a gradual deterioration of the immune system is manifested by the reduction of CD4+T-cell populations in the peripheral blood. The third phase of the disease is achieved when the concentration of the T-cells is lower than a critical value ($`30\%`$), leading to the development of AIDS. As a consequence, the patient normally dies from opportunistic diseases, which would be usually controlled by the immune system. This common pattern observed during the course of the HIV infection  is depicted in Fig. 1, which shows the plasma viremia titer and the CD4+T cell counts in the peripheral blood as functions of time.

Several theories have been proposed to explain why and how the virus remain (albeit in low concentrations) in the organism after the primary immune response, and the causes of the decline of T-cell counts, leading to the onset of AIDS. So far, none of them has provided, a complete explanation for the entire process. The most common assumption made to explain the permanence of the virus is that of fast mutation rate of the HIV, due to its rapid replication and its high degree of variability . This assumption plays a central role in our theory as well.

The challenge to understand the dynamics of the HIV infection has also attracted the attention of mathematicians and physicists that have been working on mathematical models  to describe different aspects of the dynamics of the host-parasite interaction. Although many of these models have contributed to the understanding of various aspects of the development of the disease, they fail to describe either the changes in the amount of virus detected in the body (blood, tissues and other body fluids) during the entire process or the two time scales observed in the course of the HIV infection, i.e., the primary response which runs from days to weeks and the clinical latency period that may vary from months to years. From the dynamical point of view, these different time scales may be related to two distinct kinds of interactions: one local and fast; and the other long-ranged and slow.

Most of the mathematical models proposed so far use population equations, treating the infected patient as a homogenous entity. Therefore local interactions and spatial inhomogeneities, caused by localization of the initial of immune response in lymphoid organs, are not  taken into account. We believe that this feature, which is a natural ingredient of our model, are of central importance.

Experimental evidences  support that the lymphoid tissue is a major reservoir of HIV infection in vivo. Moreover a snapshot of the distribution of cells among the different compartments of the immune system will show only a small fraction ($`24\%`$) of the cells circulating in blood and lymph, while the majority is found in the lymphoid organs . Paradoxically, the process of mobilization and activation of immune cells directed against the virus by the immune system, which occurs in the lymphoid micro-environment in this case, provides a milieu that contributes to the virus spread. .

Such spatial structure emerges naturally in discrete models, based on cellular automata that were shown  to describe well cooperative and collective patterns in experimentally observed immune response . Therefore we model the course of the HIV infection by a cellular automaton. Our model takes into account the main features of the immune response to any pathogen, the high mutation rate of the HIV and a fair amount of spatial localization that may occur in the lymphoid tissues. Our aim is to test whether the combination of these hypotheses can explain the three-phase dynamics observed on the course of the HIV infection (see Fig. 1). The results obtained by simulations of our model are shown in Fig. 2 and, as far as we know, this is the first time that the complex dynamics of the HIV infection process has been so faithfully reproduced by a theoretical model.

We model the immune system cells in the lymphoid tissues by a two-dimensional cellular automaton. To each lattice site we associate a cell, that may be a target for the HIV (as for instance, a CD4+T cell or a monocyte). Each cell can be in one of four states: (a) healthy; (b) infected-A1, corresponding to an infected cell that is free to spread the infection; (c) infected-A2, an infected cell against which a specific immune response has already developed and, finally (d) dead, an infected cell that was eliminated by the immune response.

The initial configuration is mainly composed of healthy cells, with a small fraction, $`p_{HIV}`$, of infected-A1 cells, representing the initial contamination by the HIV. In one time step the entire lattice is updated in a synchronized parallel way, according to the rules described below. The updated state of a cell depends on the states of its eight nearest neighbors.

Rule 1: Update of a healthy cell: (a): if it has at least one infected-A1 neighbor, it becomes infected-A1. (b): if it has no infected-A1 neighbor but does have at least $`R`$ ($`2<R<8`$) infected-A2 neighbors, it becomes infected-A1 (c): otherwise it stays healthy.

Rule 1a mimics the spread of the HIV infection by contact, before the immune system had developed its specific response against the virus. Rule 1b represents the fact that infected-A2 cells may, before dying, contaminate a healthy cell if their concentration is above some threshold.

Rule 2: An infected-A1 cell becomes infected-A2 after $`\tau `$ time steps.

An infected-A2 cell is one against which the immune response has developed and hence its ability to spread the infection is reduced. Here $`\tau `$ represents the time required for the immune system to develop a specific response to kill an infected cell. Such a time delay is required for each infected cell since in our model we view each new infected cell as carrying a different lineage (strain) of the virus. This is the way we incorporate the fast mutation rate of the virus in our model. When a healthy cell is infected, the virus uses the cell’s DNA in order to transcribe its RNA and replicate. During each transcription an error may occur, producing, on the average, one mutation per generation and hence a new strain of the virus is produced .

Rule 3: Infected-A2 cells become dead cells.

This rule simulates the depletion of the infected cells by the immune response.

Rule 4: (a) Dead cells can be replaced by healthy cells with probability $`p_{repl}`$ in the next time step, (or remain dead with probability $`1p_{repl}`$). (b) Each new healthy cell introduced may be replaced by an infected-A1 with probability $`p_{infec}`$.

Rule 4a describes the bone marrow replenishment of the depleted cells, mimicking the high ability of the immune system to recover from the immunosupression generated by infection. As a consequence, it will also mimic some diffusion of the cells in the tissue. Rule 4b simulates the entrance of new infected cells in the system that could come from other compartments of the immune system.

This four-state automaton allows transitions from each of its four states into the next one in a cyclic manner. We performed simulations of the model on square lattices with periodic boundary conditions of $`N=L^2`$ sites, with $`L`$ ranging from $`300`$ up to $`1000`$. All the parameters adopted in the simulations were based on experimental data. The initial concentration of HIV, $`p_{HIV}=0.05`$, was chosen based on the observation that the order of one in $`10^2`$ or $`10^3`$ T cells harbor viral DNA during the primary infection . We adopted $`p_{infec}=10^5`$ due to the fact that only one in $`10^4`$ to $`10^5`$ cells in the peripheral blood of infected individual express viral proteins. Since the immune system has a high ability to replenish the depleted cells, we used $`p_{repl}=0.99`$, although smaller probabilities could also be considered, representing different individuals. Since the time delay parameter ($`\tau `$) may vary from $`2`$ to $`6`$ weeks, we chose $`\tau =4`$. Each of our time steps corresponds to one week.

We present in Fig. 2 the evolution of the densities of healthy cells, infected cells (considering both A1 and A2 types) and dead cells, obtained in simulations of our model . We show results averaged over one hundred simulations (that used different initial configurations) and the corresponding standard deviations (shown as error bars). There is excellent qualitative agreement between our results for the density of healthy and infected cells and the time evolution of the number of CD4+ T cells in the peripheral blood and the plasma viremia titer showed in Fig. 1. The model reproduces the two time scales observed in the dynamics of the HIV infection: a short one associated with the primary infection and a long one related to the clinical latency period and the onset of AIDS. The small error bars related to the first stage of the infection indicated that its dynamics is insensitive to the initial configuration. However, the large error bars during the latency period suggest the presence of some mechanism regulating this phase of infection.

The results presented in Fig. 2 are global properties, i.e. average quantities over the entire system. Our model allows also a closer look at local behavior, which in fact, may provide the clue for understanding also the global properties. From the analysis of the spatial configurations generated by the model in various individual simulations, we noticed that the slow dynamics observed in the latency period is related to the emergence, after completion of the primary response, of some special spatial structures of infected cells. These growing special structures spread the infection in such way that they slowly commit more and more healthy cells, segregating and trapping uninfected cells.

In Fig. 3 we present four “snapshots” of typical configurations obtained during one particular simulation. Starting form an initial configuration mostly composed of healthy cells (blue) with a random distribution of infected-A1 cells (yellow), in subsequent time steps each individual infected-A1 cell generates a pulse of infected cells, of width $`(\tau +1)`$, propagating in all directions. Whenever the average distance $`l_1`$ between individual infected cells in the initial configuration is less than or equal $`(2\tau +1)`$, the independent pulses achieve a maximum coverage of the lattice. For $`\tau =4`$ the maximum coverage occurs after five weeks as shown in Fig. 3a. Note that the distribution of dead (red) cells corresponds to the initial configuration of infected cells. After that the concentration of infected cells decreases to a minimal value at $`2(\tau +1)`$ time steps, establishing the end of the primary infection phase ($`10`$ weeks, in this case).

In the following time steps the presence of infected cells will be dictated by $`p_{infec}`$, according to rule 4b. When new infected cells are introduced, they may generate two different kinds of structures. The simplest case corresponds, as discussed above, to a wave of infected cells propagating in all directions, but in this case since $`p_{infec}<<p_{HIV}`$ the average distance $`l_2`$ between the new infected cells, introduced with probability $`p_{infec}`$ at random locations, is much larger than $`(2\tau +1)`$. Examples of such structure are shown in the bottom right part of Figs. 3b and 3c. The spread of infection generated by these structures alone takes a longer time to reach nearly complete lattice coverage, and after that the infection vanishes. It is very rare, however, to find only these simple structures. A second type of structure also occurs, due to an interplay between rules 1b, 4a and 4b. These special structures are generated by ”sources” of infected cells. These sources appear, for instance, when a new infected-A1 cell, introduced by rules 4a and 4b, is surrounded by at least four dead cells. At every every ($`\tau +3`$) time steps they launch a propagating wave front of infected cells with width ($`\tau +1`$). The period of the wave fronts corresponds to the period of time the infected-A1 cell (source) remains in the same state ($`\tau `$) plus the three time steps to transition to types infected-A2, dead and healthy, in this sequence. Figs 3b and 3c (upper left) show the growing of such structures for two subsequent periods, corresponding respectively to $`18`$ and $`25`$ weeks. As these structures grow, the number of infected cells increases and the concentration of healthy cells decreases. Moreover they segregate and trap uninfected cells between two consecutive waves of infected cells. These special pattern of cell aggregations explain the slow decrease of T cells in infected patients during the latency period. These growing structures may eventually cover the entire lattice, with final densities of different cells evolving towards a steady state with average concentrations: $`(\tau +1)/(\tau +3)`$ for all the infected cells, and $`1/(\tau +3)`$ for both the dead and healthy cells. Note that the average density for healthy cells at the steady state is always lower than the observed critical value($`30\%`$) of CD4+ T cells counts related to the onset of AIDS in infected patients, i.e, the system’s breakdown. Therefore the dynamics of real experimental data is related to the transient behavior of our model and not to its steady state. The long time scale, associated with this transient period that follows the primary infection, during which the model relaxes towards its steady state, corresponds to the clinical latency period.

Analysis of the source distribution at any given time shows that the latency period depends on $`l_3`$, the average distance between sources and, consequently, on the probability of occurrence of sources ($`p_{repl}.p_{infec}.p_{HIV}`$). Actually, since new sources can be released at any time step, the length of this transient time depends on the spatio-temporal average of the distance between sources.

In this work we have shown that our cellular automaton model reproduces quite well the three-stage dynamics observed in the course of the HIV infection, and the different time scales observed in this process. The short time scale, characteristic of the primary infection, increases when $`\tau `$ is increased or when $`p_{HIV}`$ decreases. The long time scale, responsible for the clinical latency period and the onset of AIDS, is associated with the emergence of special structures, that increase slowly the number of infected cells and confine healthy cells. This special pattern formation depends on the value of the parameter $`R`$ of rule 1b, $`p_{repl}`$ and $`p_{infec}`$. Smaller values of $`p_{repl}`$ enhance the presence of dead cells and favor the formation of the special growing structures. We have also performed a careful investigation of the parameter space, and found that three-phase dynamics is reproducible for a wide range of the parameters. The complete study of the parameter space, the detalied discussion of the necessary conditions to generate the special spatial structures and the role they play spreading the infection, is under preparation and will be published elsewhere.

The growing structures of infected cells may be associated with syncytia, the aggregation of cells observed experimentally, suggesting according to our results that they are responsible for the depletion of T cells leading to AIDS. These results actually corroborate some previous suggestions that syncytia are responsible for the permanence of the virus in the system, based on the analysis of the similarities between HIV infection and other diseases .

Finally we emphasize that our cellular automaton model is, as far as we know, the first model which succeeded to reproduce the general features of the complex dynamics of the HIV infection, as observed in infected patients. This work further substantiates claims made in previous studies , that discrete models may be useful to describe emergent properties of complex biological systems, and to understand the mechanisms underlying its dynamical behavior. The reason for our success in describing the three-stage dynamics, whereas the other approaches fail, is that they do not take into account the spatial localization effects that play a major role on the course of the infection. A detailed study of the mechanisms underlying the dynamics described by our model may lead to new differential equation approaches, more suitable to describe the kind of dynamics observed in the course of HIV infection.

We thank Dr. M. Curié Cabral, I. Procaccia, S. Bocalletti, and J-P. Eckmann for enlightening discussions and to Borko Stǒsic for assistance with color graphics. RMZS’ thanks J.F. Fontanari for the hospitality during her visit to IFSC-USP, and FAPESP (project 99/09999-1) for supporting her stay. This work was partially supported by the Brazilian Agencies CNPq, CAPES, FINEP (under the grant PRONEX 94.76.0004/97) and FAPERJ.


text:
# On the structure and stability of germanium nanoparticles

## Abstract

In order to control and tailor the properties of nanodots, it is essential to separate the effects of quantum confinement from those due to the surface, and to gain insight into the influence of preparation conditions on the dot physical properties. We address these issues for the case of small Ge clusters (1-3 nm), using a combination of empirical and first-principles molecular dynamics techniques. Our results show that over a wide temperature range the diamond structure is more stable than tetragonal, ST12-like structures for clusters containing more than 50 atoms; however, the magnitude of the energy difference between the two geometries is strongly dependent on the surface properties. Based on our structural data, we propose a mechanism which may be responsible for the formation of metastable ST12 clusters in vapor deposition experiments, by cold quenching of amorphous nanoparticles with unsaturated, reconstructed surfaces.

In semiconductor nanoparticles quantum confinement leads to an increase of the optical gap compared to the bulk value and thus opens new possibilities for controlling photoluminescence effects, with narrow emission spectra tunable over a wide range of wavelengths. These properties make semiconductor dots attractive for many applications including photovoltaics, lasers and infrared dyes. Furthermore, their brightness, low toxicity and the ability to use a single excitation wavelength make them good alternatives to organic dyes for biological labelling, but their low water solubility has limited biological applications. However, recent experiments have shown that using specific coatings, the surface of selected semiconductor nanodots can be tailored to enhance the chemical interaction with a biological sample and the water solubility.

Understanding the influence of surface reconstruction and passivation on the ground state properties of semiconductor nanodots is a key prerequisite not only in designing biological applications, but also for controlling deposition of nanoparticles on surfaces and aggregation of multiple dots into new structures. In order to tailor the properties of nanodots, it is important to separate the effects of quantum confinement from those due to the surface, and to gain insight into the mechanism by which preparation conditions can influence the dot atomic structure and thus its optical properties.

Here we address these issues for the case of small Ge dots (1-3 nm), whose atomic structure is the most controversial amongst those of group IV and II-VI semiconductors. While some preparation techniques, including chemical methods, yield diamond-like Ge dots irrespective of size, several experiments suggest a structural transition, as the dot diameter becomes smaller than 4-5 nm. In particular, some experiments using vapor deposition techniques indicate a change from a cubic diamond (DIA) to a tetragonal structure, possibly ST12, in contrast to the behavior found for Si and other II-VI dots. In the bulk, the ST12 phase is only obtained from high-pressure experiments . The ST12 gap is direct and about 1.5 eV, as opposed to an indirect gap of 0.7 eV in the cubic phase. The relationship between a possible structural transition as a function of nanoparticle size and cluster preparation conditions, as well as the influence of a structural transition on the dot optical properties are as yet unknown.

Using first principles calculations we have studied the structure and stability of Ge clusters with diameters smaller than approximately 2.5 nm. We have considered both cubic diamond and tetragonal ST12 structures with H-terminated and bare reconstructed surfaces. Summarizing our results, we found that over a wide temperature range, the DIA structure is more stable than ST12-like structures for all clusters with more than 40 atoms (about 1 nm in diameter). The energy difference between the two geometries monotonically increases as a function of size, for both H-terminated and bare reconstructed surfaces. However the magnitude of this energy difference is strongly dependent on the reconstruction and on H-passivation. Our calculations also suggest that surface effects may be responsible for the formation of metastable ST12 clusters from amorphous nanoparticles in vapor deposition experiments: when the dot diameter becomes smaller than $``$ 2-3 nm, a structural transition from amorphous to ST12-like nanodots may be driven by the pressure on the inner core of the dot arising from the reconstructed surface.

In our calculations, Ge nanocrystallites were represented by free-standing clusters  in a large supercell . The total energy of the dots was computed using Density Functional Theory (DFT) in the Local Density Approximation (LDA), using iterative optimization techniques. The electronic wave functions were expanded in plane waves, with an energy cutoff of 11 Ry, and non-local pseudopotentials were used to represent the interaction between the electrons and ionic cores. We considered nanoclusters with spherical shapes, with a number of Ge atoms ranging from 28 to 300. In all cases we considered sizes which allowed us to use the same number of Ge and H atoms for both DIA and ST12 structures, in order to have direct total energy comparisons .

The energy difference between DIA and ST12 Ge clusters with H-saturated surfaces is plotted as a function of cluster size in Fig. 1 (dotted line). Our results show that H-passivated dots with a DIA-like structure are more stable than those with a ST12-like structure for all sizes, with the energy difference increasing as a function of the nanoparticle diameter. The average volume per atom ($`V`$) of the cluster, shown in Fig. 2, is reduced in comparison to the volume per atom ($`V_0`$) corresponding to the bulk first neighbor distances for both DIA and ST12 geometries; for example, the ratio $`V/V_0`$ for DIA (ST12) dots varies from 0.97 (0.96) to 0.99 (0.98) when going from a 45 to a 145 atom cluster. Over the same range of sizes, the effective pressure acting on the cluster core, as evaluated using bulk moduli data, decreases from about 2 to 1 GPa. The magnitude of the DIA-ST12 energy difference, as well as the contraction of the average cluster atomic volume are significantly modified in the presence of unsaturated, reconstructed surfaces, as described below.

In order to find reconstructed geometries for the cluster surfaces we used a combination of empirical and DFT-LDA techniques. Formation of facets is expected for clusters larger than 5 nm, but facets are unlikely in spherical nanoclusters with 1-5 nm diameters. In addition clusters prepared by deposition on surfaces usually exhibit disordered, defected interfaces. We therefore chose to determine surface reconstructed geometries using an annealing procedure. First, using molecular dynamics with a Tersoff potential , we melted the cluster surface by heating it up to 2000 K for 0.08 ns; the temperature was then slowly decreased to zero over 0.2-0.5 ns. During this phase of the calculation, the crystalline core of the nanoclusters was kept frozen, and the shape conserved by confining the system in a spherical cavity. The annealing/quenching series was repeated 2-4 times . The final structure was fully relaxed within DFT-LDA, all of the atoms being allowed to move. During these relaxations, energy gains varied from 310 (275) meV/atom for 95 atom clusters to 88 (81) meV/atom for 300 atom clusters, for dots with a DIA (ST12)-like core structure. Fig. 3 shows the surface structure and the crystalline core for a selected cluster (Ge190).

As shown in Fig. 1 (dotted versus dashed lines), we observed a strong reduction of the energy difference \[$`\mathrm{\Delta }E(N)`$\] between DIA and ST12 clusters when the H-passivated surface is replaced by one which is unsaturated and reconstructed. In the absence of H atoms, and given the spherical shape of the clusters, the computed values of $`\mathrm{\Delta }E(N)`$ can be fitted by separating surface and bulk contributions: $`\mathrm{\Delta }E(N)=N\mathrm{\Delta }ϵ+(36\pi )^{\frac{1}{3}}N^{\frac{2}{3}}\mathrm{\Delta }\gamma `$. Here $`\mathrm{\Delta }ϵ=87`$ meV/at and $`\mathrm{\Delta }\gamma =+59`$ meV/at are the volume and surface contributions to the energy difference $`\mathrm{\Delta }E(N)`$, respectively. The value of $`\mathrm{\Delta }ϵ`$ is smaller than the calculated energy difference between solid DIA and ST12 (for the bulk energy difference we find 130 meV/atom, in good agreement with previous calculations). As indicated by the value of $`\mathrm{\Delta }\gamma =+59`$ meV/at, surface energy is smaller for ST12 clusters and its sign is opposite to that of the volume contribution to $`\mathrm{\Delta }E(N)`$. This is largely responsible for the reduction of $`\mathrm{\Delta }E`$ in dots with reconstructed surfaces, compared to H-saturated clusters. Our results indicate that the ST12 structure should be more stable than DIA for $`N40`$; however for such small sizes, Ge clusters are not expected to exhibit bulk-like geometries, but rather to form complex, non spherical shapes.

As mentioned above, the average atomic volume of a dot with a bare reconstructed surface is reduced in comparison to that of a cluster with a hydrogenated surface . As plotted in Fig. 2, the ratio V/V0 is 0.93 (0.947) for ST12 (DIA) Ge190 clusters. The volume reduction remains strong for our biggest Ge300 clusters, with 0.954 (0.967) for ST12 (DIA). These last values are to be compared, e.g., with the value of 0.982 (0.988) for ST12 (DIA) in the case of the smaller Ge145H108 dot, which has roughly the same number of atoms belonging to the crystalline dot core as Ge300. The average volume reduction of the cluster amounts to an effective pressure on the crystalline core of about 4 and 2.3 GPa, for Ge190 and Ge300, respectively. The effective pressure is slightly higher for ST12 than for DIA geometries (see Fig. 2).

An analysis of the reconstructed surfaces reveals disordered structures in all cases, as expected from a fast quench from a liquid state. For the larger clusters, the bond angles range approximately from 63  to 144 , and the average bond length is close to the calculated first neighbor-distances in amorphous Ge (2.46 Å), i.e. 2 % larger than in the crystalline DIA structure. In general, we observed a strong reduction of the undercoordinated surface atoms after reconstruction, due to atomic dimerization. This effect is stronger for ST12 than for diamond, with ST12 reconstructed nanoclusters exhibiting approximately 20% fewer dangling bonds. This circumstance is due to the smaller size of the ST12 crystalline core and to the broader distribution of bond angles in the bulk ST12 structure, both of which provide greater freedom in the rearrangement of surface atoms. For one cluster size, we studied the effect on the dot stability of H-passivating the dangling bonds of its reconstructed surface. We observed an increase of $`\mathrm{\Delta }E`$, the DIA structure being favored over ST12 (the full relaxation of selected hydrogenated reconstructed surfaces were carried out within DFT-LDA).

As a final step in our study of the stability of Ge nanodots, we estimated the effect of temperature on $`\mathrm{\Delta }E(N)`$, for clusters with reconstructed surfaces, by computing free energy differences in the harmonic approximation. The vibrational free energy $`F_{vib}=_{i=1}^{3N6}\left[\frac{\mathrm{}\omega _i}{2}+k_BT\mathrm{ln}\left(1\mathrm{exp}(\frac{\mathrm{}\omega _i}{k_BT})\right)\right]`$ was determined by computing the vibrational frequencies $`\omega _i`$ using the Tersoff potential. Although not as accurate as total energy differences obtained within DFT-LDA, $`F_{vib}`$ can be used to estimate finite temperature effects as a function of size . Our results, shown in Fig. 1 by open dots, indicate that energy differences between DIA and ST12 are slightly reduced. However temperature effects do not invert the relative stability of the two structures. For examples, at $`N=145`$ a temperature greater than 1180 K, i.e. close to the melting point of Ge, would be required for the reconstructed ST12 cluster to be more stable than DIA.

Our total energy calculations have shown that cubic diamond is the most stable structure of both H-terminated and bare reconstructed Ge clusters in the 1-3 nm size range, despite the importance of surface reconstruction effects. Ge dots with the ST12 structure are metastable and it is interesting to investigate whether there exist experimental conditions which might give rise to metastable ST12 clusters. For example, in vapor deposition or sputtering experiments, amorphous Ge nanoparticles are initially present and annealing treatments are usually required for crystallization to occur. It is therefore relevant to understand whether metastable ST12 nanoparticles can be quenched from amorphous dots.

Based on our calculations, the cores of ST12 and DIA dots with both H-terminated and reconstructed surfaces are compressed. The effective pressure on the dot cores is much larger in the presence of reconstructed surfaces. These results suggest that pressure effects may play a role in quenching metastable ST12 clusters from amorphous nanoparticles. In order to address this issue, we have first investigated the amorphous (a-Ge) to ST12 transition in bulk Ge. Fig. 4 shows the total energy of diamond, ST12 and a-Ge as a function of volume, at T=0, as obtained from our calculations. Both a-Ge and ST12 are metastable, with the amorphous phase being slightly lower in energy than the ST12 crystal. A pressure Pt of 1.5 GPa is required to induce an a-Ge to ST12 transition. Whether such a transition actually occurs depends on the height of the barrier between the two structures and on the temperature. We have not attempted to compute the a-Ge to ST12 energy barrier; however, phenomenological arguments suggest that it should be lower than that between a-Ge and cubic diamond, at temperatures typical of, e.g. dot deposition experiments. Indeed, ST12 is a weakly ordered crystal, which has been used to model a-Ge: it has 12 atoms per unit cell and a space group with few symmetry operations . Most importantly, unlike diamond, the ST12 crystalline network exhibits 7- and 5-fold atom rings, similar to a-Ge. It is therefore conceivable that a transition between a-Ge and ST12 be possible at relatively low T, when a-Ge is under a pressure of 1.5 GPa or higher. Similar considerations for the case of nanodots suggest that the pressure exerted by reconstructed surfaces on amorphous nanoparticle cores initially present in vapor deposition experiments may be large enough to induce a transition from amorphous to ST12 metastable nanoparticles. An extrapolation of the calculated effective pressures on dot cores (see Fig. 2) suggest that for dots with bare reconstructed surfaces and a diameter smaller than 2.5-3 nm, the pressure on the crystalline core is larger than the pressure required in the bulk to induce an a-Ge to ST12 transition (2.5-3 GPa, as compared to 1.5 GPa). Therefore for dots with diameters smaller than 2.5-3.0 nm prepared in vapor deposition experiments, an a-Ge to ST12 transition induced by an effective surface pressure may be possible. On the contrary, the pressure exerted on the core of H-passivate clusters is equal to or smaller than the bulk transition pressure even for clusters with 70-100 atoms (i.e. with a diameter less than 1.5 nm).

The LDA used in our work does not permit quantitative evaluations of optical gaps for Ge dots. However, it is interesting to note that the difference between the energy of the highest occupied orbital (HOMO) and lowest unoccupied orbital (LUMO) as obtained for H-passivated clusters is larger for DIA-like geometries than for ST12-like geometries. While the HOMO position in energy is similar in DIA and ST12 clusters, the LUMO position is higher for DIA than for ST12 clusters. If the same trend was to be confirmed for quasi-particle energies, then a measure of the optical gap of small Ge dots could be a way to discriminate between DIA and ST12 geometries. Work is in progress to go beyond the LDA and estimate optical gaps.

In conclusion, we have shown that Ge clusters with the diamond structure are more stable than tetragonal ST12 dots over a wide temperature range, irrespective of the cluster size, for dot diameters larger than $``$ 1.5 nm. We have proposed a mechanism which may be responsible for the formation of metastable ST12 clusters in vapor deposition experiments, by cold quenching of amorphous nanoparticles with unsaturated, reconstructed surfaces. The pressure exerted on the nanoparticle core by the surface can induce an amorphous to ST12 transition, for clusters with diameters smaller than 2.5-3.0 nm. This may explain why different types of structures are seen in experiments using chemical preparation methods versus physical vapor deposition methods. According to our calculations, chemical methods should always yield diamond structures, consistent with the results of H. W. Lee et al.. Our study indicates that quantum confinement as well as surface effects are both key features in understanding the physical properties of small semiconductor dots. By tuning the surface properties with, for example, a particular choice of surfactant or by otherwise controlling the surface reconstruction, the pressure exerted on the dot core can be modified and used to tailor the atomic structure of the dot and indirectly the electronic properties.

Useful discussions with L. Terminello, C. Bostedt, A. Van Buuren, H. Lee, S. Bastea and E. Schwegler are gratefully acknowledged. This work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory, Office of Basic Energy Sciences, Division of Materials Science, under contract No. W-7405-Eng-48.


text:
# Comment on “Determination of Interlayer Diffusion Parameters for Ag/Ag(111)”

In a recent Letter  Roos and Tringides (RT) collected evidence in support of the hypothesis that the preexponential factor $`\nu _s`$ in the interlayer jump rate on Ag(111) exceeds the corresponding quantity $`\nu _t`$ for inlayer transport by about two orders of magnitude. As part of the argument they present a simplified analysis of an experiment carried out by Bromann et al. (BBRK) , in which second layer nucleation on top of predeposited Ag islands on Ag(111) was investigated. Here I point out a mistake in the analysis of RT, and show that the correct application of their idea does not allow to conclude that $`\nu _s/\nu _t>1`$. I then discuss more broadly our present understanding of interlayer diffusion on Ag(111) in view of a recent reanalysis  of the BBRK experiment.

RT estimate the probability $`f`$ that an atom deposited during the second dose descends from the island within the time $`\tau `$ between successive deposition events. They write $`f=(\lambda /d)p`$, where $`\lambda =\sqrt{D_t\tau }`$ is the diffusion length, $`d`$ the island diameter and $`p`$ the probability that an atom poised at an edge site jumps down from the island. The ratio $`\lambda /d`$ is referred to by RT as “the number of edge interrogations”. In fact this number is much larger than $`\lambda /d`$. Provided that $`\lambda d`$, which is the case of interest here, it is given by the number of diffusion jumps $`D_t\tau =\lambda ^2`$ multiplied by the fraction $`L/A`$ of edge sites among all sites on the island; $`L`$ is the island perimeter and $`Ad^2`$ the island area in units of the lattice constant, so $`(L/A)D_t\tau L(\lambda /d)^2\lambda /d`$.

To obtain the correct expression for $`f`$, note that the probability that an atom on the island descends in a small time interval $`dt`$ is $`(L/A)D_s(1f)dt`$, where $`D_s=pD_t=(\nu _s/\nu _t)D_t\mathrm{exp}(\mathrm{\Delta }E_s/kT)`$ is the interlayer hopping rate. Integrating up to time $`t=\tau `$ then yields

$$f=1\mathrm{exp}[(L/A)D_s\tau ],$$
(1)
which corrects Eq.(1) of RT. Using the numbers given by RT, for the case of islands of radius 30 $`\mathrm{\AA }`$ at temperature $`T=130`$ K, we obtain $`f=1\mathrm{exp}[130(\nu _s/\nu _t)]`$ for $`\mathrm{\Delta }E_s=0.13`$ eV. Thus $`f=1`$ provided $`\nu _s/\nu _t0.1`$.

In addition to the BBRK experiment, RT base their conclusion on the analysis of two other growth situations involving interlayer transport. In all three cases they work at a single temperature, which implies that information about $`\nu _s`$ can be extracted only if $`\mathrm{\Delta }E_s`$ is known. The value of $`\mathrm{\Delta }E_s`$ used by RT was obtained in  by analyzing the dependence of the fraction of islands with second layer nuclei (a quantity somewhat similar to $`1f`$) as a function of the radius of predeposited islands, at two different temperatures $`T=120`$ K and 130 K. The analysis was based on a theory of second layer nucleation due to Tersoff et al. , which has recently been shown to be quantitatively incorrect . Other groups  have estimated $`\mathrm{\Delta }E_s`$ assuming that the adatom density at second layer nucleation is comparable to that at which first layer islands nucleate, which is generally not true .

A reanalysis  of the BBRK data using the correct expression for the rate of second layer nucleation yields a very large step edge barrier $`\mathrm{\Delta }E_s0.32`$ eV, accompanied by a prefactor $`\nu _s8\times 10^{19}`$ s<sup>-1</sup>, which would imply $`\nu _s/\nu _t4\times 10^8`$. As such a large preexponential factor seems hard to justify physically, this suggests that a single pair of diffusion parameters is insufficient to describe interlayer transport on Ag(111), and that additional second layer nucleation experiments at variable temperature are called for. A temperature-dependent measurement of interlayer transport based on the decay of vacancy islands was performed by Morgenstern et al., who estimate $`\mathrm{\Delta }E_s=0.13`$ eV and $`\nu _s/\nu _t=10^{0.6\pm 0.5}<1`$ . The discrepancy between these numbers and the (correctly analyzed) BBRK experiment indicates that, despite considerable effort, interlayer diffusion for Ag(111) remains an open problem.

Useful correspondence with Michael Tringides is gratefully acknowledged. This work was supported by DFG within SFB 237.


text:
# Simulation of Flux Lines with Columnar Pins: Bose Glass and Entangled Liquids

The seminal paper of Abrikosov in 1957  showed, within a mean field calculation, that for fields below a critical field $`H_{c2}(T)`$ a type II superconductor allows magnetic flux to penetrate in the form of a triangular array of vortices each carrying a unit flux quantum $`\varphi _0=hc/2e`$. This picture is modified considerably in the case of high $`T_c`$ cuprates which are quasi two-dimensional and have high transition temperatures . The most dramatic consequence of the enhanced thermal fluctuations on the Abrikosov lattice is its melting into a vortex liquid by a first order transition as indicated in resistivity, neutron scattering, specific heat, SQUID magnetometry and local Hall probe  measurements. Peak effect  measurements on NbSe2 have also found a reentrant phase diagram in the H-T plane . Theoretical calculations have primarily been confined to density functional approaches on the layered system of pancake vortices  valid for BSCCO, lattice simulations of frustrated 3D XY-models , Monte Carlo or Langevin dynamics of the vortices , and determination of the melting line using Lindemann criterion . Simulations of vortices, treated as line or polymer-like objects, have been done using an analogy with interacting bosons  by variational and path integral techniques .

The effect of pins on the structure and dynamics of the vortex lattice adds another dimension to vortex behavior. The original work of Larkin and Ovchinnikov  predicted that for $`d<4`$, point pins destroy the long range order of the vortex lattice and break it up into small ordered coherence volumes. From the point of view of practical applications where the aim is to minimize the motion of vortices in order to reduce dissipation, it is found that columnar pins along the vortex direction are extremely effective in pinning vortices . There have been very few simulations attempting to study the phases and phase transitions in the presence of columnar pins. The ones based on 3D XY models assume an underlying numerical grid which affects the vortex structures for the small system sizes and high vortex densities that have been studied .

With this motivation, in this Letter we present the first simulations of the structure and thermodynamics of up to 1000 vortices with and without strong columnar pins using a continuum path integral quantum Monte Carlo (PIMC) simulation of the corresponding interacting bosons in (2+1)D. There is no bias of a superimposed lattice and no prior input of the nature of wave functions or the type of the transition. We extend previous simulations  to include columnar pins on considerably larger lattices.

Our main findings are as follows– In the pure system, we see the reentrant nature of the melting line in agreement with the predictions of Nelson . We find a triangular solid phase at low temperatures $`T`$ which melts into an entangled liquid phase at higher $`T`$ by a first order transition with a finite entropy jump. The defect density, and structure factor at the first reciprocal lattice vector, $`S(𝐆)`$, also show sharp discontinuities at the transition. In the solid phase, $`S(𝐆)𝒪(N)`$ and drops sharply to $`𝒪(1)`$ in the liquid phase. In the presence of 10% columnar pins, the vortex lattice transforms into a Bose glass phase at low $`T`$. The Bose glass consists of patches of ordered regions with only short range positional and orientational order. This phase is qualitatively different from the Bragg glass phase with point pins. With increasing temperature, there is a transition from a Bose glass to an entangled ‘defected liquid’ which is considerably smeared.

Consider a system of $`N`$ flux lines in 3D characterized by a bending energy $`\stackrel{~}{ϵ_1}=\{\varphi _0/(4\pi \lambda )\}^2`$  where $`\lambda `$ is the London penetration depth. The inter-vortex interaction is approximated as $`V(r_{ij})=ϵ_0K_0(r_{ij}/\lambda )`$ where $`ϵ_0=\varphi _0^2/(8\pi ^2\lambda ^2)`$ is the scale of the inter-vortex potential, $`K_0`$ is the modified Bessel function of the first kind and $`𝐫_i(z)`$’s are the 2D position vectors of the flux lines in a plane perpendicular to their length at a height $`z`$. $`K_0(r)\mathrm{ln}(r)`$ as $`r0`$ and is valid in the high vortex density limit, the regime studied in Ref. . In this paper, we fix the density on the lower part of the reentrant melting curve  and therefore retain the full $`K_0(r)`$ form of the potential.

The partition function of the flux lines is mapped onto the world lines or Feynman paths of bosons  described by the Hamiltonian

$$H=\mathrm{\Lambda }^2\underset{i=1}{\overset{N}{}}_i^2+\underset{i<j}{}K_0(r_{ij})$$
(1)
The bosonic path integrals are in two spatial dimensions and one imaginary time direction, which is equivalent to the thickness of the sample along the direction of the magnetic field. All lengths are measured in units of $`\lambda `$ and energies in units of $`ϵ_0`$. The dimensionless de Boer parameter $`\mathrm{\Lambda }=k_BT/(2\stackrel{~}{ϵ_1}\lambda ^2ϵ_0)^{1/2}`$ is a measure of quantum fluctuations (the relative strength of the kinetic energy vs the interaction energy) of the boson system, or equivalently a measure of temperature $`T`$ in the vortex system. Once the mapping is established, we do a PIMC simulation  of interacting bosons at an inverse boson temperature $`\beta =1/T_b`$ . Periodic boundary conditions are used in all directions and permutations between bosons are included. If $`T_b<\mathrm{\Lambda }^2\rho `$, where $`\rho =N\lambda ^2/(\mathrm{Area})`$ is the dimensionless density of the bosons, exchanges between the particles are possible and can generate entangled vortex lines.

In Fig. 1(a) we show the structure of 1000 vortices at low temperature $`\mathrm{\Lambda }=0.045`$ where the vortex system shows a beautiful triangular lattice with small excursions of the vortex lines about their equilibrium positions seen in Fig. 2(a). For $`\mathrm{\Lambda }>0.064`$ the vortices melt into an entangled liquid as seen from the actual configurations of the vortex lines in Fig. 2(b). Concomitantly, the superfluid density $`\rho _s`$ of the corresponding boson system jumps from zero in the solid phase to a finite value in the liquid phase. The entanglement correlation length , which is the typical length along the $`z`$-direction for vortex lines to cross, is on the order of the sample size in the solid and drops to $`4a_0`$, where $`a_0`$ is the planar lattice constant, in the liquid phase.

The degree of positional order is quantified by the structure factor $`S(𝐤=𝐆)`$ at the first reciprocal lattice vector of the triangular lattice. As seen in Fig. 3, with increasing temperature there is a sharp transition from a crystalline to a liquid phase at $`\mathrm{\Lambda }_m0.064`$. In a temperature-driven first order transition, the free energy is continuous across the transition but the internal energy and entropy are discontinuous. In analogy in the mapped boson problem at $`T=0`$, the total energy is continuous but there are sharp jumps in the potential and kinetic energies between $`\mathrm{\Lambda }=0.062`$ and $`0.064`$. The melting transition is located by the intersection of the polynomial fits to the energy in the liquid and solid phases at $`\mathrm{\Lambda }_m=0.064`$ in good agreement with Ref. .

In the crystal $`S(𝐤)`$ shows Bragg peaks with $`S(𝐆)𝒪(N)`$, where $`N`$ is the number of vortices, while above $`\mathrm{\Lambda }_m`$, $`S(𝐤)`$ is liquid-like with $`S(𝐆)𝒪(1)`$. By fitting the structure factor data to the form $`S(𝐆)=1+(N1)\mathrm{exp}(G^2u^2)`$ for a finite system we have calculated the Lindemann number $`c_L=u^2^{1/2}/a_0=0.24`$ at the transition just within the solid phase, in good agreement with other works, where $`u^2`$ is the mean squared displacement of the particles from their lattice positions .

In addition to the positional correlations, the orientational order for six-fold symmetry is described by $`\mathrm{\Psi }_6(𝐫_i)=(1/6)_je^{i6\theta _{ij}}`$, where the sum over $`j`$ runs over all the nearest neighbors of $`i`$ in a triangulation picture, and $`\theta _{ij}`$ is the angle the bond between particles at $`𝐫_i`$ and $`𝐫_j`$ makes with an arbitrary reference axis. The orientational correlation function defined as $`g_6(r)=\mathrm{\Psi }_6(r)\mathrm{\Psi }_6^{}(0)`$ also shows long range order in the clean system at low $`T`$. We find that the correlation lengths for positional and orientational order can reliably be extracted only from data on large system sizes $`1000`$ vortices and we find that both types of order vanish at $`\mathrm{\Lambda }_m`$.

A measure of imperfections in the triangular arrangement of the flux lines is the number of topologically-defected non six-fold coordinated sites in the Delaunay triangulation of imaginary time slices of the system. This defect density is small at low temperatures in the clean system and shows a sharp jump across the melting transition (Fig. 4).

Effect of Columnar Pins:

In the boson picture, randomly placed potential wells are correlated in the imaginary-time direction and correspond to columnar pins in the vortex representation. We model the columnar pin as an attractive delta function which traps a vortex all along its length. In the boson picture the trapped vortex corresponds to a fixed classical particle. The vortex phase in the presence of columnar pins has been previously referred to as a Bose glass phase . Our simulations provide the first detailed description of the structure and properties of this so-called Bose glass. In Fig. 1(b) we show the configuration of the 1000 vortex system at $`\mathrm{\Lambda }=0.045`$ with $`10\%`$ of the vortices, those trapped by columnar pins, held fixed. The pins are found to seed or nucleate an ordered crystalline patch, however, since the pins are randomly located, the different patch orientations are not commensurate and frustrates the translational and orientational order of the lattice. Columnar pinning dramatically reduces $`S(𝐆)`$ from around 250 in the pure 1000 vortex system to about 6 in the presence of pins at $`\mathrm{\Lambda }=0.045`$ in the Bose glass phase (Fig. 5(a)). In addition, $`S(k)`$ has a jagged peak, with shoulders and split peaks, which are characteristic of glasses. The Bose glass is rather different from a structural glass or a frozen liquid which has a peak height of $`S(𝐆)1.5`$. At higher temperatures the Bose glass melts into a defected liquid phase where the unpinned vortices are entangled (Fig. 2(c)). In the parameter region investigated, we do not find any evidence for a disentangled liquid either with or without pins. The Bose glass phase that we find in these simulations is in a different limit from that studied by Nelson and Vinokur  where the number of columnar pins is much larger than the number of vortices.

The topological defects generated by the columnar pins in the vortex lattice destroy the positional and orientational order within few lattice spacings (Fig. 5(b)) in contrast to a Bragg glass  in the presence of weak point pins which retains long range orientational order and algebraic positional order.

Columnar pinning has a marked effect on the first order melting transition. The kinetic and potential energies vary smoothly through the transition. The sharp drop in $`S(𝐆)`$ in Fig. 3 also disappears. The sharp jump in defect density is replaced by a much smoother increase through the transition as seen in Fig. 4. The local variations in the melting temperature caused by the quenched random potential below a certain critical dimension can cause a sharp first order transition to get rounded off in presence of quenched randomness . It is possible that at lower pin densities, the first order jumps in various quantities persist with reduced magnitudes. Further simulations are in progress to study the evolution of the sharp first order transition with increasing defect density.

P. Sen would like to acknowledge B. Militzer, G. Bauer and E. Draeger for help with the code. We would also like to thank S. Banerjee, S. Bhattacharya, C. Dasgupta, T. Giamarchi, A. Grover, G. Menon, and A. Paramekanti for useful discussions. We acknowledge the support of the NCSA, UIUC for computational resources and software and the NSF grant no. DMR 98-02373.

e-mail: ntrivedi@tifr.res.in

FIGURE CAPTIONS


text:
# DETECTION OF COSMIC SHEAR WITH THE WILLIAM HERSCHEL TELESCOPE

## 1 Introduction

Understanding the large-scale distribution of matter in the universe continues to be a major issue in modern cosmology. Weak gravitational lensing promises to be a particularly effective method for determining properties of large-scale structure, since it provides direct information concerning the total mass distribution, independently of its state and nature.

The images of distant field galaxies obtained at a telescope are slightly coherently distorted, due to weak lensing by large-scale structure. With extensive measurements of this shear on various scales, one would obtain a direct measure of the power spectrum of density fluctuations along the line of sight.

However, the first stage in such a programme is the detection of the cosmic shear signal; this itself is a challenge, because the rms shear amplitude is small - a few percent on arcminute scales. Recently four papers describing the detection of this effect have been released (Wittman et al 2000, van Waerbeke et al 2000, Bacon et al 2000a, Kaiser et al 2000), presenting mutually consistent results with careful analysis of systematic effects.

Here we overview the detection of cosmic shear obtained with the 4.2m William Herschel Telescope, fully discussed in Bacon, Refregier & Ellis (2000a). The current paper describes the survey strategy, and discusses how the data are analysed to overcome convincingly the contribution of systematic effects to the shear. Simulations used to check our methodology are explained, and the cosmological implications of our results are discussed.

## 2 Observations

The goal of our survey is to obtain a homogeneous sample of deep fields, chosen to be on random lines of sight separated by $`>5^{}`$ in order to sample independent structures. The galactic latitude of the fields was tuned to afford $`200`$ stars within the field of view, necessary to correct for anisotropic PSF systematics.

We carried out deep $`R`$-band imaging on 14 such fields with the Prime Focus camera on the WHT. This has an 8’$`\times `$16’ field of view and pixel size 0.237”. The fields were exposed for one hour in $`R`$; median seeing was 0.81”, having excluded exposures with seeing $`>`$ 1.2”. We reach a magnitude depth of $`R_{\mathrm{median}}=25.2`$, with $`R_{\mathrm{median}}`$ of 23.4 for our selected sample, $`z_{\mathrm{median}}0.8`$, and a number density $`N=14.3`$ arcmin<sup>-2</sup>.

## 3 Analysis of Systematic Effects

The aim of our analysis is to remove carefully systematic effects from the galaxies’ measured ellipticities, leading to unbiased measures of the small ($``$1%) mean shear components for each field. We wish to measure the mean shear in 8’$`\times `$8’ cells, for increased shear signal and to allow cross-correlation tests.

We used the Kaiser, Squires & Broadhurst (1995) method (KSB) implemented by Kaiser’s imcat software for object detection and shape measurement. A detection signal-to-noise $`\nu >15`$ limit was imposed on our usable galaxy catalogue to remove correction systematics found at low signal-to-noise level.

The shear induced by telescope optics must be dealt with; by calculating the telescope distortion using objects’ relative positions in several dithers, we show that this is negligible (shear due to telescope $`<0.003`$ everywhere; see figure 1).

Next, the PSF anisotropy from e.g. tracking errors must be removed (see figure 2). Before correction this effect induces an rms stellar ellipticity $`e=0.07`$, but after subtracting a fitted 2-dimensional cubic to stellar ellipticities, the residual stellar rms is a mere $`e=1.4\times 10^3`$. We correct the galaxy ellipticities following Luppino & Kaiser (1997), using the stellar fit model and responsivities to smear measured for the galaxies.

Finally, galaxies are corrected for isotropic smear, i.e. the fact that smaller galaxies’ shapes have been more affected by seeing-induced circularisation. Again we follow Luppino & Kaiser’s (1997) method; this results in estimates for the mean shear components in each 8’$`\times `$8’ cell.

## 4 Shear Measurements on Simulated Data

In order to check our correction of systematics, and to calibrate KSB-measured shear to real shear, we constructed simulated WHT fields on which to carry out the above shear analysis (Bacon et al 2000b); one can apply a chosen shear to a field, and test its recovery by our algorithm. By creating a joint probability model of the magnitude - number density - ellipticity - radius distribution of galaxies in Ebbels’ (1998) HST Groth Strip survey catalogue, we were able to draw out statistically similar simulated catalogues for shearing and analysis.

The catalogues were visualised with IRAF artdata; telescope-specific pixelisation, throughput, anisotropic PSF, poisson and readout noise were added. A null set of 20 simulated fields without shear were created, together with a further 30 fields with an rms shear of 1.5%. The KSB analysis described above was carried out on each field.

## 5 Results

We shall now compare the detection results for the simulated and real data. The left-hand panel of figure 3 shows the mean shear components found for our 1.5% rms shear simulations (30 cells). The inner circle represents the variance that would be expected from noise alone; one can see that there is an excess variance $`\sigma _{\mathrm{lens}}^2`$, which turns out to be significant; $`\sigma _{\mathrm{lens}}^2=(0.013)^2\pm (0.006)^2`$, to be compared with our input $`\sigma _{\mathrm{lens}}^2=(0.015)^2`$. The fact that the method has detected the simulated signal with the correct amplitude is encouraging.

The middle panel of figure 3 shows the mean shear components for 20 simulations with no shear added. Note that here the variance is accounted for by noise alone; as expected, there is no excess variance.

The shear results for our observed fields are shown in the right-hand panel of figure 3. Again we see an excess variance; with a thorough statistical analysis we find this to be significant, with $`\sigma _{\mathrm{lens}}^2=(0.016)^2\pm (0.008)^2\pm (0.005)^2`$, the errors being due to noise and uncertainty on any remaining systematics respectively. This corresponds to a 3.4$`\sigma `$ detection of cosmic shear.

By comparing our $`\sigma _{\mathrm{lens}}^2=(0.016)^2\pm (0.012)^2`$ (error now includes cosmic variance) with that expected for popular cosmological models, we find that COBE-normalised SCDM is ruled out at 3$`\sigma `$ level, whereas cluster-normalised $`\tau `$-, $`\mathrm{\Lambda }`$-, and OCDM are highly consistent with the data.

Our results also afford us a measure of $`\sigma _8`$ for a given cosmological model. For instance, for $`\mathrm{\Lambda }`$CDM with $`\mathrm{\Omega }_m`$ = 0.3, we obtain $`\sigma _8=1.47\pm 0.51`$, consistent with cluster abundance determinations, $`\sigma _8=1.13\pm 0.19`$ (Viana & Liddle 1996).

With increased numbers of fields in the future, cosmic shear variance measurements will afford very precise estimates of $`\sigma _8`$, while skewness measurements of the distortion field will provide an independent estimate of $`\mathrm{\Omega }_m`$.

## 6 Conclusion

Evidence for the detection of shear arising from large-scale structure has been presented based on an analysis of 14 fields obtained at the William Herschel Telescope. Particular attention has been paid to questions of systematic correction and testing of measurement method by simulations. Prospects are now bright for measuring with greater accuracy the amplitude of the cosmic shear signal. The key uncertainties to overcome are noise and cosmic variance, i.e. more independent lines of sight will lead to a better estimate of the shear amplitude. Future cosmic shear surveys will consequently provide powerful constraints on cosmological parameters.

## Acknowledgments

Particular thanks go to Nick Kaiser for providing the Imcat software, and to Douglas Clowe for valuable discussions as to its use. We also thank Roger Blandford, Chris Benn, Andrew Firth, Mike Irwin, Konrad Kuijken, Peter Schneider, Andrew Liddle, Yannick Mellier, Roberto Maoli and Jason Rhodes for useful discussions. Max Pettini kindly provided us with one of the WHT fields. This work was performed within the European TMR lensing network.

## References


text:
# Charge profile in vortices

\[

## Abstract

The electric charge density in the vortex lattice of superconductors is studied within the Ginzburg-Landau theory. We show that the electrostatic potential $`\phi `$ is proportional to the GL function, $`\phi |\psi |^2|\psi _{\mathrm{}}|^2`$. Numerical results for the triangular vortex lattice are presented.

\] Abrikosov vortices in type-II superconductor consist of magnetic flux encircled by a rotating condensate of Cooper pairs. The inertial and Lorentz forces both act in the centrifugal direction and thus cause a small depletion of the charges in the vortex core. The electrostatic potential generated by this depleted charge supplies a centripetal force which balances the centrifugal force .

While the depleted charge occupies a region comparable with the vortex core, i.e. a region extending over the Ginzburg-Landau (GL) coherence length $`\xi `$, the screening charge extends over a region extending over half the London penetration depth $`\lambda `$. This charge distribution results in an electrostatic field which may be observed either by nuclear quadrupole resonance  or by surface scanning as proposed recently by Blatter et al.  or by the Torricelli-Bernoulli effect . Besides, the charge of the vortex core is expected to contribute to the forces acting on the vortex. For example, Feigel’man et al.  speculate that this charge can explain the still puzzling sign reversal of the Hall voltage.

The charge in the asymptotic region away from the vortex core has been studied by LeBlanc . In this region one can benefit from the simplicity of the London theory. Within the London theory, the local balance of forces acting on the Cooper pair may be written in form of the Newton equation,

$$m^{}\dot{𝐯}=e^{}\phi +e^{}𝐯\times 𝐁.$$
(1)
Here the mass $`m^{}=2m`$ and charge $`e^{}=2e`$ are twice the values of a single electron. Unlike in common applications of the Newton equation (1) here the velocity is given by the London relation, $`m𝐯=e𝐀`$, and Eq. (1) thus specifies the electric field $`\phi `$ which maintains this motion. Using the identity $`\dot{𝐯}=(𝐯)𝐯=\frac{1}{2}\left(𝐯^2\right)𝐯\times \times 𝐯`$ and the London approximation $`e𝐁=m\times 𝐯`$, one can show that Eq. (1) is solved by the Bernoulli potential,

$$\phi =\frac{mv^2}{2e}.$$
(2)
The quadratic dependence on the velocity shows that the electric potential $`\phi `$ decreases over a characteristic length $`\lambda /2`$. The charge density $`\rho =\epsilon ^2\phi `$ behaves similarly.

Obviously, the Bernoulli potential (2) is not suited for the core region, where the local London approximation fails. Moreover, close to the core the velocity diverges, $`v\frac{e\mathrm{\Phi }_0}{2\pi m}\frac{1}{r}`$, where $`\mathrm{\Phi }_0`$ is the quantum of flux and $`r`$ is the distance from the vortex axis. According to Eq. (2) and $`\rho =\epsilon ^2\phi `$, the total screening charge follows this divergence, $`_r^{\mathrm{}}dr\rho r^2`$. The charge of the core estimated from the London expression (2) thus crucially depends on the estimate of the core diameter and the profile of the charge is incorrect in the core region.

In the core region, a non-local approach like the Ginzburg-Landau theory is required . Like the London approach, the Ginzburg-Landau theory does not include the electrostatic potential explicitly. Here we show how this potential can be determined from its effect on the motion of Cooper pairs. The Bernoulli potential then has to be replaced by the potential

$$\phi =\frac{\beta }{2e}\left(|\psi |^2|\psi _{\mathrm{}}|^2\right),$$
(3)
with $`\beta `$ and $`\psi _{\mathrm{}}`$ having the standard meaning in the Ginzburg-Landau theory .

The proof of (3) parallels the proof of the Bernoulli potential. The Cooper pairs move according to the Schrödinger equation,

$$\frac{1}{2m^{}}\left(i\mathrm{}+e^{}𝐀\right)^2\psi +e^{}\phi \psi =0,$$
(4)
which replaces the force balance equation (1).

The electrostatic potential $`\phi `$ also causes a mismatch between the local Fermi level and the global chemical potential that results in a local suppression of Cooper pairs. Briefly, the charge densities of superconducting and normal electrons have a ratio different from the unperturbed state, $`\rho _\mathrm{s}:\rho _\mathrm{n}\rho _\mathrm{s}^{\mathrm{}}:\rho _\mathrm{n}^{\mathrm{}}`$, where the superscript $`\mathrm{}`$ denotes values far from the perturbation. Due to this mismatch, the share of normal electrons is locally increased. In the linear approximation one has

$$\rho _\mathrm{n}=\rho _\mathrm{n}^{\mathrm{}}\frac{e^2}{D}\phi .$$
(5)
Taking the gradient, one can read (5) as balance between the electric force and a kind of osmotic pressure of normal electrons, $`e𝐄=D(\rho __\mathrm{n}/e)`$. In this sense, equation (5) describes the response of normal electrons to the electric field, importance of which has already been noticed by London . The “diffusion coefficient” $`D`$ is identified below.

Finally, the electrostatic potential obeys the Poisson equation,

$$\epsilon ^2\phi =\rho =\rho _\mathrm{s}+\rho _\mathrm{n}+\rho _{\mathrm{lattice}}.$$
(6)
With the wave function normalized to the charge density of the condensate, $`\rho _\mathrm{s}=e^{}|\psi |^2`$, we can use the condition of neutrality, $`e^{}|\psi _{\mathrm{}}|^2+\rho _\mathrm{n}^{\mathrm{}}+\rho _{\mathrm{lattice}}=0`$, to eliminate the normal density from the Poisson equation,

$$\epsilon ^2\phi =e^{}\left(|\psi |^2|\psi _{\mathrm{}}|^2\right)\frac{e^2}{D}\phi .$$
(7)
The set of equations (4) and (7) is closed by the Maxwell equation, $`\times \times 𝐀=\mu 𝐣`$, with the superconducting current $`𝐣=\frac{e^{}}{m^{}}\mathrm{Re}[\overline{\psi }(i\mathrm{}+e^{}𝐀)\psi ]`$.

Let us suppose that the system is nearly neutral so that the left hand side of (7) may be disregarded. In this case, the potential can be determined from the suppression of the superconducting density,

$$\phi \frac{2D}{e}\left(|\psi |^2|\psi _{\mathrm{}}|^2\right).$$
(8)
Substitution of (8) into the Schrödinder equation (4) yields a closed equation for $`\psi `$,

$$\frac{1}{2m^{}}\left(i\mathrm{}+e^{}𝐀\right)^2\psi +4D|\psi |^2\psi 4D|\psi _{\mathrm{}}|^2\psi =0.$$
(9)
One can see that for $`4D=\beta `$ and the standard GL relation, $`\alpha =\beta |\psi _{\mathrm{}}|^2`$, Eq. (9) turns into the GL equation,

$$\frac{1}{2m^{}}\left(i\mathrm{}+e^{}𝐀\right)^2\psi +\beta |\psi |^2\psi +\alpha \psi =0.$$
(10)
With $`D=\beta /4`$, Eq. (8) coincides with Eq. (3).

Now we show that the charge neutrality is a good approximation in the general case. A simple estimate shows that the deviation from charge neutrality, $`\rho =\epsilon ^2\phi `$, is much smaller than the increase of the normal electron charge density $`4e^2\beta ^1\phi `$. Indeed, the potential varies on the scale of the coherence length, $`^2\phi \xi ^2\phi `$, therefore the ratio of these two charge densities behaves as

$$\frac{\epsilon ^2\phi }{4e^2\beta ^1\phi }\frac{\beta \rho _\mathrm{s}^{\mathrm{}}\kappa ^2}{emc^2}\frac{k_\mathrm{B}^2T_c^2\kappa ^2}{ϵ_Fmc^2}\frac{|\rho _\mathrm{s}^{\mathrm{}}|}{|\rho _{\mathrm{lattice}}|}\frac{\mathrm{}^2\lambda ^2}{m^2c^2\xi ^4}.$$
(11)
Here we have used the BCS value of $`\beta `$ , see Eq. (13), and introduced the GL parameter $`\kappa =\lambda /\xi `$. Even for extreme values of the high-$`T_c`$ materials \[$`\kappa 100`$, $`T_c90`$ K, the Fermi energy $`ϵ_F0.08`$ eV, $`m4m_e`$, $`\epsilon =4\epsilon _0`$, $`\rho _\mathrm{s}^{\mathrm{}}\rho _{\mathrm{lattice}}`$\] the estimate (11) gives about $`10^5`$. For conventional superconductors the ratio (11) is typically about $`10^9`$. Accordingly, the deviation from the charge neutrality may be disregarded when one solves for the wave function. With $`\epsilon ^2\phi 0`$ the Poisson equation (7) simplifies to Eq. (3).

Now we compare the above result with alternative approaches. First, far from the core one can make the local approximation of the Ginzburg-Landau equation (10)

$$\frac{m^{}v^2}{2}\psi +\beta |\psi |^2\psi +\alpha \psi =0,$$
(12)
where the velocity is given by the current, $`𝐣=e^{}𝐯|\psi |^2`$. One can see that in this limit the electrostatic potential found from (12) and (8) equals the Bernoulli potential (2). We note that the incompressible superconducting liquid discussed by London  corresponds to the limit $`\kappa \mathrm{}`$ of the GL theory . In this limit, the local approximation is exact.

Second, formula (3) is similar to the formula proposed by Khomskii and Kusmartsev , $`\phi _{\mathrm{KK}}=\mathrm{\Delta }^2/(4eϵ_F)`$. It turns out, however, that from the BCS-Gor’kov relations ,

$$\beta =\frac{24(\pi k_\mathrm{B}T_c)^2}{7\zeta (3)ϵ_FN},\psi =\frac{\sqrt{7\zeta (3)N}}{4\pi k_\mathrm{B}T_c}\mathrm{\Delta },$$
(13)
one obtains

$$\frac{\beta }{2e}|\psi |^2=\frac{3\mathrm{\Delta }^2}{4eϵ_F},$$
(14)
i.e., our value is three times larger than $`\phi _{\mathrm{KK}}`$. The similarity of (14) with the potential of Khomskii and Kusmartsev is, however, accidental. While (14) follows from the intertial and Lorentz forces, the potential $`\phi _{\mathrm{KK}}`$ results from the condensation energy of the pairing. We have neglected the condensation energy in our treatment.

In the rest of this paper we show actual charge profiles of an isolated vortex and of the vortex lattice. We use known solutions of the Ginzburg-Landau equation and evaluate the electrostatic potential $`\phi `$ from Eq. (3). The deviation of the charge density from neutrality then follows from the Poisson equation, $`\rho =\epsilon ^2\phi `$.

We start the discussion using the isolated vortex model suggested by Clem ,

$$|\psi |^2=|\psi _{\mathrm{}}|^2\frac{r^2}{r^2+\xi _\mathrm{v}^2},$$
(15)
where $`\xi _\mathrm{v}\xi `$ is given by the minimum of the free energy as $`\xi _\mathrm{v}=\xi \sqrt{2}\sqrt{1K_0^2(\xi _\mathrm{v}/\lambda )/K_1^2(\xi _\mathrm{v}/\lambda )}`$, where $`K_0`$ and $`K_1`$ are modified Bessell functions. This wave function fails at long distances, $`r\lambda `$, where the correct behavior follows the local approximation, $`|\psi |^2=|\psi _{\mathrm{}}|^2(1+\frac{m^{}v^2}{2\alpha })`$, but it provides a good approximation in the vicinity of the core, $`r\xi `$, where the dominant charge contribution is located. The electric field, $`𝐄=\phi `$, obtained from (15) and (3) is directed radially,

$$𝐄=\frac{\beta }{e}|\psi _{\mathrm{}}|^2\frac{\xi _\mathrm{v}^2}{(r^2+\xi _\mathrm{v}^2)^2}𝐫.$$
(16)
The charge density corresponding to (15),

$`\rho `$ $`=`$ $`2\epsilon {\displaystyle \frac{\beta }{e}}|\psi _{\mathrm{}}|^2{\displaystyle \frac{\xi _\mathrm{v}^2(\xi _\mathrm{v}^2r^2)}{(r^2+\xi _\mathrm{v}^2)^3}}`$ (17)
$`=`$ $`{\displaystyle \frac{a_\mathrm{B}e}{\pi }}{\displaystyle \frac{\xi _\mathrm{v}^2}{2\xi ^2}}{\displaystyle \frac{\xi _\mathrm{v}^2r^2}{(r^2+\xi _\mathrm{v}^2)^3}},`$ (18)
is depleted in the core, $`r<\xi _\mathrm{v}`$, and enhanced outside the core, $`r>\xi _\mathrm{v}`$. Here, $`a_\mathrm{B}=\frac{4\pi \epsilon \mathrm{}^2}{me^2}`$ is the Bohr radius in the crystal.

The second form of (18) can be compared with the result of Blatter et al. , who found a more complicated factor, $`\rho _{\mathrm{Bl}}\rho \times \pi ^2\frac{dT_c}{d\mathrm{ln}ϵ_F}`$. We note that Blatter et al. approached the problem from the microscopic side using the BCS relation between the total charge and the chemical potential. This charge is moreover screened in the spirit of Jakeman and Pike  within the Thomas-Fermi approximation. Since $`\frac{dT_c}{d\mathrm{ln}ϵ_F}\mathrm{ln}(\mathrm{}\omega _D/T_c)`$ 1 to 10, their numerical results are close to ours.

The charge of the vortex core per unit length, $`Q=2\pi _0^{\xi _\mathrm{v}}r\rho dr=2\pi \epsilon E(\xi _\mathrm{v})\xi _\mathrm{v}`$, results as

$$Q=\pi \epsilon \frac{\beta }{2e}|\psi _{\mathrm{}}|^2=\frac{\pi \epsilon \alpha }{2e}=\frac{a_\mathrm{B}e}{4\xi ^2}.$$
(19)
In a single pancake (a vortex of length $`d=1.1710^9`$ m) in $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ (YBCO) at 77 K, the depleted charge is $`Qd=0.004`$ e. We have used $`a_\mathrm{B}=0.05`$ nm and $`\xi =1.91`$ nm, where the coherence lenght follows from $`\xi =\lambda /\kappa `$ with $`\lambda =191`$ nm and $`\kappa =100`$. With the quasiparticle screening in the spirit of van Vijfeijken and Staas, this value will be reduced by $`\rho _\mathrm{s}/\rho =2|\psi _{\mathrm{}}|^2/N1(T/T_c)^4=0.5`$ for $`T=77`$ K.

Note that the vortex charge (19) has a value independent of the core diameter $`\xi _\mathrm{v}`$. A similar value, $`Q=1.28\times a_\mathrm{B}e/(4\xi ^2)`$, results for a different model function, $`|\psi |=|\psi _{\mathrm{}}|\mathrm{tanh}(r/\xi )`$, recommended in . This seems to contradict the expectation from the local approximation, where the core diameter plays a dominant role. In the local approach, the velocity close to the core reads $`v\frac{e\mathrm{\Phi }_0}{2\pi mr}`$. From the Bernoulli potential (2) and the Poisson equation one finds that inside the cylinder of radius $`r`$ the local approximation predicts a charge

$$Q_{\mathrm{la}}=\frac{e\epsilon \mathrm{\Phi }_0^2}{2\pi mr^2}.$$
(20)
With the core radius equal to the coherence length, $`r=\xi =\mathrm{}/\sqrt{2m^{}\alpha }`$, the classical estimate (20) reproduces the nonlocal result (19). Although the exact agreement of Eq. (20) with Eq. (19) is rather accidental, it shows that classical estimates used in the literature are reasonable.

An example of the electrostatic potential for the triangular vortex lattice is presented in Fig. 1. One can see that in the core region the potential has a shape similar to the Clem model. The maximum amplitude, $`\phi _{\mathrm{max}}=1.3`$ mV, is identical to the maximum potential obtained from the Clem model since the GL function $`|\psi |^2`$ reaches zero at the vortex center in both cases. The GL function used in Fig. 1 has been evaluated by the iteration method described in .

The density of charge corresponding to the potential in Fig. 1 is shown in Fig. 2. In the core, the charge is negative, i.e., the density of charge carriers is depleted. Outside the core, the charge is positive, i.e. there is charge built up which screens the charge of the core. The radius of the vortex core (radius at which the charge changes sign) found in Fig. 2 is $`1.53\xi `$, comparable to the value $`1.41\xi `$ obtained from the Clem model. The maximum of the screening charge appears at $`2.4\xi `$, while the Clem model yields $`2\xi `$. The maximum of the depleted charge in the core centre is $`510^6`$ electron/cell while the Clem model gives $`6.810^6`$ electron/cell. The core charge per unit lenght is thus reduced to 85% of the Clem model. These numerical results thus show that the approximation based on the Clem model is reasonable inside and in the vicinity of the core.

In conclusion, we have shown that the GL function can be used to evaluate the electric field in superconductors, see Eq. (3), in spite of the fact that the electrostatic potential is absent in the GL theory. For convenient estimates of the core charge we have derived the simple formula (19). Our theory neglects contributions from the condensation energy and the screening by normal particles; more work on this is under way. Due to the neglected screening, the theory applies only for intermediate temperatures where the GL theory still gives reasonable results while the share of the superconducting electrons is already of the order of unity. Numerical results relevant to YBCO at the temperature 77 K and at low magnetic field with well separated vortex cores have been presented.

This work was supported by MŠMT program Kontakt ME160 and GAČR 202000643, GAAV A1010806 and A1010919 grants. The European ESF program VORTEX is gratefully acknowledged.