format: text

text:
# Reversible B↔A Transitions in Single DNA Molecule Immersed in A Water Drop

## Introduction

The transformation from B-DNA Watson:53  to dehydrated A-DNA Franklin:53  in fibers was one of the first reversible structural transitions observed in a biomolecule. These structures represent right-handed helical duplexes with identical topologies and hydrogen bonding, and yet they have strongly different shapes Saenger:84 ; Calladine:92 . The B helix is long and narrow, with stacked base pairs forming the core of the helix. In contrast, the A-DNA structure is short and thick, with strongly inclined base pairs wrapping around a 6 Å solvent accessible cylindrical hole. The B-form is the dominant biological conformation of DNA whereas the A-form is considered as a high energy state that some sequences can adopt temporarily during various biological functions Ivanov:95 . Indeed, local A-DNA motives are rather common in crystal structures of protein-DNA complexes Timsit:99 ; Lu:00 . The A$``$B polymorphism is currently considered as one of the modes for governing protein-DNA interactions Ng:00 .

A complete B$``$A transition in a long mixed sequence DNA can be obtained in vitro in a few special conditions. In crystalline Franklin:53 ; Leslie:83  and amorphous Piskur:95  fibers as well as in films Tunis-Schneider:70 , it is induced by equilibrating the samples under reduced relative humidity or with the concentration of organic solvents increased to $``$80% Zimmerman:83a . Organic solvents can also provoke B$``$A transition in an isolated DNA molecule in solution Ivanov:83 . The transition is reversible and cooperative, Ivanov:73 ; Ivanov:74  which means that it occurs concertedly in an extended DNA fragment rather than randomly in separate nucleotides. The characteristic length of such fragment, however, is is rather small. Originally, it was estimated as 20 base pairs Ivanov:74 , but later this estimate was reduced to 10 base pairs Ivanov:83 , that is approximately one helical turn. DNA fragments shorter than this length generally do not go to the A-form in solution Fairall:89 ; Galat:90 . For some sequences, however, short DNA fragments could be crystallized in the A-form Conner:82 ; NDB: . Under physiological conditions in vitro, the A/B balance depends upon the environment, notably, the types of solvent counterions Ivanov:73  and the temperature Nishimura:86 . It is also sequence dependent, with some sequences exhibiting in solution the features of both forms Fairall:89 ; Wolk:89 ; Lindqvist:01 . The transition is significantly facilitated in groups of consecutive guanines (G-tracts) whereas in A-tracts it is more difficult Leslie:83 ; Peticolas:88 ; Ivanov:96 . The A/B-philicity of different base pair steps and triplets has been parametrized and used successfully for predicting the properties of mixed sequences Ivanov:95 ; Basham:95 ; Tolstorukov:01 .

Several different ideas concerning detailed molecular aspects and the driving forces of A$``$B transitions in DNA were discussed in different years Ivanov:73 ; Alden:79 ; Calladine:84 ; Calladine:92 ; Saenger:86 ; Hunter:93 ; Cheatham:97b ; Cheatham:97c ; Jayaram:98 . Both the sugar phosphate backbone Saenger:86  and base pair stacking Calladine:84  were proposed as one of the possible actors in the transition, with water as an evident second partner. Water bridges between certain DNA atoms may stabilize either B-DNA Calladine:84  or A-DNA Saenger:86 , and the hydrophobic effect possibly favors a B-DNA type of stacking Ivanov:73 ; Hunter:93 . It has been suggested that the reduced water activity is a universal factor that shifts the equilibrium towards A-DNA Malenkov:75 , which was later replaced by a less specific ”economics of hydration” Saenger:86 . Although the B$``$A transition induced by increased salt concentration has been reported for poly-dG, the low water activity in high salt generally does not induce B$``$A transitions Ivanov:73 ; Zimmerman:80 . In any case, it is not clear how the reduced amount of water around DNA affects its structure at the atomic level Fuller:88 . Molecular dynamics simulations suggest, in addition, that accumulation of metal cations in the major groove may specifically stabilize A-DNA Cheatham:97b ; Cheatham:97c .

The cooperativity of the B$``$A transition is poorly understood. Its physical origin was initially attributed to the sugar puckering Ivanov:73 ; Ivanov:74 . The furanose ring conformation is C3’-endo (North) in A-DNA (with the C3’-carbon puckered above the furanose plane towards the nucleobase), and C2’-endo (South) in B-DNA Saenger:84 . For a single nucleoside in vacuum, the C2’-endo conformation is slightly lower in energy, with the barrier of 2-3 kcal/mole between the two states Foloppe:99 , and this is the only significant local energy barrier between the A- and B-DNA conformations. Until recently, it was assumed that the switch in the sugar pucker is strictly linked with the base pair orientation. The latter is restrained due to stacking interactions with the neighbors, therefore, the sugar conformations can change only simultaneously in at least several consecutive base pairs. These views, however, were always questioned Calladine:84  and the recent X-ray and spectroscopic studies really demonstrated that the North sugar puckers are accessible without the B$``$A transition Fairall:89 ; Wolk:89 ; Lindqvist:01 , and even with all sugars in the C3’-endo conformation the overall B-DNA structure can persist Ng:00 . It appears, therefore, that the sugar conformations are only weakly restrained by the stacking and the observed cooperativity can hardly be due to the sugar pseudorotation barrier. At the same time, the ensemble of the currently available X-ray DNA structures agrees with the cooperative character of the B$``$A transition. Although some A-DNA structures have B-like features near the termini Malinina:99 , and there exist a series of structures classified as A/B transition intermediates Vargason:01  the A- and B-forms generally do not mix in the same DNA fragment. The two forms may even coexist in a crystal, but in different samples of the same fragment Doucet:89 .

Many aspect of the A/B polymorphism were earlier studied by theoretical methods. Its thermodynamics is well described by an analytical theory based upon the Ising model Ivanov:74 . Reduced Calladine:84 ; Hunter:93 ; Mazur:89  and all atom Ivanov:79  models in continuous dielectric were used for studying its mechanism and sequence dependence. On the other hand, different analytical theories combined with simplified DNA representation were used for modeling solvent effects Jayaram:96 . Steady progress in MD simulations of DNA since the pioneering calculations by Levitt Levitt:83b  recently made possible detailed modeling of A/B transitions in realistic environment including explicit water and metal counterions Cheatham:97b ; Cheatham:97c ; Jayaram:98 ; Cheatham:96 ; Yang:96 ; Cieplak:97 ; Sprous:98 ; Langley:98 ; Soliva:99 . They shed new light upon the putative molecular mechanisms involved in the DNA A/B polymorphism despite limited time of MD trajectories and some forcefield artifacts Feig:98 .

In the present study I try to get new insight in this problem by using a different MD approach. It has recently been shown that the particle-mesh Ewald (PME) algorithm allows one to follow the dynamics of a DNA double helix put in a small salty water drop Mzjacs:02 . An evident question one asks, what happens if the size of the drop is reduced? By analogy with drying fiber samples it should be expected that DNA would go from B to A-form. A series of simulations described here confirms that this is really the case. Although these conditions have never been used in experiments our simulations give the best currently possible prediction of the properties of such systems. Moreover, the transition is reversible in the time scale accessible in the calculations, that is both B$``$A and A$``$B transitions can be observed in the nanosecond time range. This allows one to test different hypotheses concerning the factors governing the transition. The results suggest that the A$``$B transitions in an isolated water drop are driven by medium range electrostatic interactions between DNA backbone and counterions accumulated in the major groove. Moreover, several evidences suggest that A$``$B transitions in other conditions are actually governed by the same mechanism.

## Methods

In the series MD simulations presented here the Dickerson-Drew double helical dodecamer (CGCGAATTCGCG Wing:80 ) and its derivatives are used as examples. Similar results were also observed for several other test sequences. This particular DNA fragment is preferred because it is neither A-philic or B-philic and because its dynamic structure with Cornell et al. force field Cornell:95  is rather close to experimental data Cieplak:97 ; Young:97b ; Duan:97 . The simulation protocols were similar to the earlier water drop simulations Mzjacs:02 . We use the internal coordinate molecular dynamics (ICMD) method Mzjcc:97 ; Mzbook:01  adapted for DNA Mzjacs:98 ; Mzjchp:99  with the time step of 0.01 ps. The electrostatic interactions are evaluated with the shifted Coulomb law by using the SPME method Essmann:95 ; Mzjacs:02  and the long range cut-off of 50 Å which was always larger than the system size. The Van-der-Waals and SPME direct sum interactions were truncated at 9 Å with the value of Ewald parameter $`\beta 0.35`$.

The initial state for B$``$A transitions was prepared as before Mzjacs:02  with the canonical B-DNA Arnott:72  used a standard conformation. The DNA molecule was first immersed in a large rectangular TIP3P Jorgensen:83  water box of and next external solvent molecules were removed by using a spherical distance cut-off from DNA atoms. The cut-off radius was adjusted to obtain the desired number of water molecules remaining. The drop was neutralized by randomly adding the necessary number of Na+ ions. Simulations of A$``$B transitions started from the final structure of a B$``$A transition obtained with the lowest hydration.

Every system was energy minimized first with the solute held rigid and then with all degrees of freedom. Dynamics were initiated with the Maxwell distribution of generalized momenta at 250K and equilibrated at this temperature during several picoseconds. Production trajectories were computed with the temperature bound to 300 K by the Berendsen algorithm Berendsen:84  with a relaxation time of 10 ps. For better comparison with earlier simulations of A$``$B transitions, the original Cornell et al. force filed Cornell:95  was used. Several test simulations showed, however, that qualitatively similar results are obtained with the subsequent modifications of the AMBER force filed Cheatham:99 ; Wang:00 . Duration of production runs varied from 5 to 25 ns depending upon the character of dynamics observed. The conformations were saved with a 2.5 ps interval. Programs Curves, Curves:  XmMol XmMol:  and Mathematica by Wolfram Research Inc. were employed in the data processing.

## Results and Discussion

An example of the molecular transformation we study is shown in Fig. 1. The top row exhibits the canonical B- and A-DNA structures. The process begins with plate 1 where we see a snapshot of the first saved state during a production run started from the canonical B-DNA conformation. This system involves a dodecamer DNA fragment with 22 Na+ and 400 water molecules. It is seen that after the equilibration some ions have already entered the major groove while a thin water shell covers the surface of the DNA molecule. One may notice also a small narrowing in the middle of the minor groove which is a well-known distinctive feature of the crystal structure of this dodecamer. This initial trend is only temporary, however. The number of water molecules in the system is about 17 per nucleotide, which is less than the limiting hydration number for B-DNA Saenger:84 . In dynamics, water tries to form a smooth spherical drop around the DNA fragment and this is achieved owing to a B$``$A transition. Plate 2 shows a snapshot after 2 ns of dynamics. It is seen that the DNA molecule went to the characteristic A-form conformation with a wide minor and narrow major grooves. Note that a substantial part of the Na+ ions is sandwiched in the major groove between the two opposite sugar-phosphate strands. The more compact A-form structure allows water to form a nearly spherical drop that covers the DNA fragment. After the transition the structure remains stable during the time accessible in calculations.

The final state of the previous trajectory was used as a starting point for the second run to obtain an inverse A$``$B transition. For this purpose the system was surrounded by 600 additional water molecules to increase their number to 1000. Plate 3 in Fig. 1 displays a snapshot from the beginning of the second production run. Now the number of water molecules is about 42 per nucleotide, which corresponds to B-DNA hydration numbers. The last plate in Fig. 1 shows a snapshot after 2 ns of dynamics. As we see, the DNA molecule really went back from A- to B-DNA conformation and the narrowing in the middle of the minor groove re-appeared.

The sequence of transitions shown in Fig. 1 suggests that there exists a critical level of hydration characterized by an equal probability of A and B-forms, if the transition is cooperative, or by some structures intermediate between the two. In order to get a rough estimate of this midway point, a series of MD simulations was carried out with the size of the water drop systematically varied. The results obtained are summarized in Table 1. The corresponding trajectories were started from standard A-DNA and B-DNA states and continued as long as necessary to reach either A$``$B or B$``$A transition. For intermediate hydration values, notably, for 21 and 25 water molecules per nucleotide, very long trajectories were necessary. In all cases a transition in one direction only was obtained. For instance, with 25 water molecules per nucleotide an alternative trajectory starting from A-DNA was continued to 10 ns, but no transition occurred. The results summarized in Table 1 suggest that the hydration of 25 water molecules per nucleotide can serve as a rough estimate of the transition point.

As seen in Fig. 1, the A-DNA structure obtained in dynamics is more compact than the canonical conformation. The major groove width in the structure averaged over the last nanosecond is about 2 Å below the canonical value. The final A-DNA structures obtained for higher hydration levels had slightly wider major grooves. Note that, in Table 1, their rmsd from the canonical A-DNA is lower. In all cases, however, structures with the major groove width similar to that in the canonical A-DNA were observed only temporarily in intermediate transition phases, and they seemed to be unstable with respect to the more compact conformation shown in plate 2 of Fig. 1. Visual inspection reveals that in such structures the phosphate groups opposed across the major groove are similarly close or even closer to each other than the neighboring phosphates in the same strand. They are linked by multiple water bridges and at least one such bridge is usually replaced by a Na+ ion slightly shifted either inside or outside the helix. With intermediate hydration numbers, however, this pattern was less clear, with both Na+ and the opposite phosphate groups often hydrated separately.

One of the trajectories mentioned in Table 1 is detailed in Fig. 2. This simulation started from the canonical B-DNA with a water shell of 500 molecules (20.8 per nucleotide). Plate (a) exhibits time plots of several structural parameters. Three of them, namely, the backbone length and the average width of the two grooves deserve preliminary comments. The backbone length is computed as the average distance between the centers of consecutive sugar rings. Note that its value is larger in A-DNA than in B-DNA because, in a B$``$A transition, the sugars are moved to the outer surface of the double helix, with the helical diameter of the sugar trace increased by nearly 6 Å. The widths of both grooves are evaluated by perpendicular cross-sections of spatial traces of C5’ atoms Mzjmb:99 . Near the helical termini, however, the grooves are not measurable and, in fact, the groove lengths vary in dynamics. The last feature is very clear in the dynamics of the major groove profile shown in plate (c). The average widths were computed for the zones where the corresponding groove was defined in each given structure.

The B$``$A transition seen in Fig. 2 is relatively slow, with the final structure established only after 18 ns. An early B$``$A transition occurred close to one of the termini at around 10 ns. However, it was incomplete and inverted two nanoseconds later. According to the two rmsd traces, the final transition started at 13 ns. Until that the dynamics sampled midway conformations between A- and B-forms, which is not that far from normal B-DNA dynamics with the same potentials Cieplak:97 ; Young:97b ; Duan:97 . It may be seen that the minor groove profile periodically switched from narrowing to widening in the middle. At about 13 ns the two rmsd traces drastically moved close to their final levels. Simultaneously the major groove width fell by almost 10 Å and went below the A-DNA level. At that time the molecule already looked very much like A-DNA, with the instantaneous rmsd-A values around 2 Å, but almost all sugars still had their puckers in the South region. The final phase of the transition is slower and during this phase the sugars in the middle AATT tetraplet switch to the North. In this middle zone the transition to the A-form is complete whereas the three terminal base pairs remain in an intermediate conformation with South sugar puckers.

Although the structural parameters shown in plate (a) differ by the amplitudes of fluctuations they all reveal a slow overall drift suggesting a downhill motion of the system on a free energy surface. In contrast, the total energy remains stable throughout the trajectory except the initial phase. This means that, already during heating, the system finds a balance of internal interactions maintained afterwards, and, consequently, the observed steady downhill character of the transition has an entropic origin. The precise nature this entropic force is not clear. The backward A$``$B transitions have qualitatively similar kinetic suggesting that this feature is not due to the properties of the A and B-DNA forms, but results from the specific setup of simulations.

In the backward A$``$B transitions, the sequence of events is inverted, that is, the sugar puckers quickly switch to the North whereas the rmsd’s and the groove widths need sometimes many nanoseconds to reach stationary B-DNA values. Such kinetics does not depend upon the amount of water added and it is slower than in earlier reported AMBER simulations Cheatham:96 ; Sprous:98 . The apparent reason is that during the B$``$A transition a large number of metal ions enter the major groove when it is still wide. A high ion concentration in the major groove is one of the factors that supposedly stabilize the A-DNA conformation Cheatham:97b ; Cheatham:97c . In contrast, when the trajectory is to be started from the canonical A-DNA, the major groove is narrow from the beginning, and very long equilibration period should be necessary for the same number of ions to reach optimal positions.

The possibility of B$``$A transitions in longer and shorter DNA fragments has been checked in a series of simulations summarized in Fig. 3. The foregoing results suggested that the B$``$A transitions in a water drop are mainly driven by electrostatic interactions between counterions and backbone phosphates across the major groove. A characteristic A-form major groove can be formed starting from eight base pairs and, as seen in Fig. 2c, its length in a dodecamer A-DNA is between for and five nucleotides. One could expect, therefore, that B$``$A transitions would progressively become more difficult as the DNA is shortened. All trajectories shown in Fig. 3 started from B-DNA with the lowest degree of hydration used in Table 1. The number of Na+ and the size of the water drop were changed according to the DNA chain length. As expected, relatively rapid B$``$A transitions are observed in the original dodecamer and its longer derivatives. In contrast, the transition did not occur in shorter fragments. The hexamer always remained in a typical B-DNA conformation with high rise and negative inclination. The octamer and decamer fragments exhibited dynamics somewhat similar to that of the dodecamer, with accumulation of Na+ ions in the major groove and noticeable compression of the DNA structure, however, the B$``$A transition did not occur. Moreover, one of the terminal base pairs in the decamer was broken, with water bridges formed between the two nucleotides, which apparently helped to move closer the opposite phosphate strands in the major groove. This effect was reproduced in a repeated run with hydrogen bonds in terminal base pairs enforced by weak additional restraints. The last observation suggested that the three shorter fragments in Fig. 3 would not go to A-form even in much longer dynamics.

Based upon the described qualitative character of DNA dynamics, one can consider two main physical factors as the possible driving forces of the B$``$A transition under low hydration. First, water tends to reduce its contact surface with vacuum and forms a spherical drop forcing the DNA fragment to shrink. This force can be considered as surface tension or as a particular case of the hydrophobic effect, with vacuum in place of a nonpolar solvent. The second factor is the increased concentration of Na+ ions, which reduces the water activity. These two alternatives were checked in additional simulations under modified conditions. To check separately the role of the surface tension, the counterion effects were suppressed by discharging the DNA fragment and next running dynamics without counterions. The partial charges of OP atoms were reduced by by 0.5 eu and two long trajectories were computed starting from the canonical B-DNA structure under the lowest and the highest hydration, respectively. The DNA structures sampled in both these simulations turned out to be similar. These were slightly underwound B-DNA conformations with the average twist around 30. Under the low hydration the duplex looked slightly compressed, but no signs of a B$``$A transition could be noticed.

The effect of the high counterion concentration upon the water activity was separately checked in two simulations with added NaCl. Earlier it was reported that high NaCl does not cause B$``$A transitions in conventional simulations with periodical boundaries Cheatham:97c . To make a more stringent test, here the NaCl concentration was increased under the midway hydration as estimated in Table 1. It was hoped that under such conditions an increased salt concentration can produce a clear qualitative effect upon the dynamics. The initial A and B-DNA states in a water shell of 600 molecules were modified by replacing 66 randomly chosen water molecules with 33 Na+ and 33 Cl- ions. After that, 66 water molecules were added to maintain the earlier hydration level. The effective concentration of Na+ was thus raised from 2M to 5M. The two systems were re-equilibrated with the standard protocol and long trajectories were computed to obtain either A$``$B or B$``$A transition. In contrast to the tests in Table 1, here an A$``$B transition was obtained after about 13 ns dynamics as shown in Fig. 4. The dynamics exhibited in this figure is typical of the A$``$B transitions described above although its kinetics are much slower than those in Table 1. The second trajectory was similarly long, but only B-DNA structures were sampled.

In the animations of these trajectories, progressive crystallization of NaCl was seen, with large blocks of typical cubic NaCl lattice found in the major DNA groove in the end of both trajectories. Addition of NaCl is supposed to stimulate A$``$B transitions in fiber DNA crystalsCooper:66 . Moreover, NaCl crystallization occurs in DNA fibers at 75% relative humidity if the salt content is above 6% of the DNA weightCooper:66 , which gives 4.8M effective concentration of NaCl according to the reported water/DNA adsorption curvesFalk:62 . One may note, therefore, that the overall effect of the added salt in our simulations qualitatively agrees with experiments. The observed shift towards the B-form may be caused by the NaCl crystallization either directly, by forcing the major groove to open, or indirectly, by pushing water out from the major groove together with some free Na+ ions. In any case, however, these results suggest that it is not the water activity reduced by ions that causes the B$``$A transitions in small drops.

The last simulations suggest that B$``$A transitions are induced by a united effect of the limited water shell and counterions, and not by one of these factors considered separately. A simple and attractive mechanism that accounts for all the above results is as follows. The B$``$A transition is caused by direct electrostatic interactions between mobile metal cations and phosphate groups in the opening of the major DNA groove. When the local cation concentration exceeds a certain critical value, the inter-phosphate repulsion is inverted to attraction. A similar idea is long discussed in the context of DNA bendingMirzabekov:79 ; Levene:86 ; Rouzina:98  and it was earlier proposed by Cheatham and Kollman as one of the possible mechanisms responsible for stability of the A-form Cheatham:97b ; Cheatham:97c . The strong inter-phosphate attraction immediately results in local DNA bending and, eventually, in a B$``$A transition at the opposite side of the double helix. When the length of the major groove exceeds one helical turn the apparent bend becomes uniform shrinking. The role of the limited water shell consist in pushing the counterion cloud surrounding DNA inside the diameter of the double helix. The counterions avoid loosing their hydration shells and are forced to come closer to DNA. The surface tension at the drop boundary plays a certain role, but it is the reduced overall space available for counterions that seems more important. This simple model explains well our results, and below we consider why and how it may work in other conditions as well.

### Concluding Discussion

Clarification of the detailed mechanisms involved in the DNA polymorphism is an important challenge for rapidly progressing computational methods in molecular biophysicsCheatham:00 . Here, for the first time, reversible B$``$A transitions were obtained in silico in a simulated titration experiment by smooth variation of water content. The environment conditions used were not earlier encountered in experiments, nevertheless, the qualitative features found seem credible and likely to exist in analogous experimental conditions. First, the computed pattern of DNA transformations is physically sensible and agrees with known experimental trends, notably, as regards the principal effect of the amount of bound water. Second, these calculations give the most accurate currently possible prediction of the true properties of DNA fragments in small water drops. Really, the Cornel at al. force field Cornell:95  was earlier shown to produce DNA dynamics in good agreement with experiments Cheatham:00 . Its parameters were fitted with small molecule data only, and they were never specifically adjusted for periodical boundary conditions. In our case, no artificial interactions are introduced, while those between periodical images are eliminated Mzjacs:02 . There are all reasons to believe, therefore, that our simulations reproduce correctly the dynamics of DNA transformations in such unusual conditions.

Qualitatively, the results obtained here agree with earlier studies. Cheatham and Kollman were the first to obtain an A$``$B transformation in pure water, and stable A-DNA dynamics in 85% EtOH Cheatham:96 ; Cheatham:97b . Visual inspection of the reported A-DNA snapshots Cheatham:97a  shows that our structures are quite similar to those in EtOH as regards the width of the major groove, preferred positions of Na+ ions, and a B-like character of the DNA termini, which gives the molecule a strong apparent bend to the major groove. Moreover, the reported stabilization of A-DNA in EtOHCheatham:97b  required that the molecule was first covered by water and then by EtOH as an outer shell. The necessity of such double phase hydration was confirmed by Sprous et al. Sprous:98 . Although water and EtOH could mix in dynamics, it seems evident that, in order to be stable in simulations, the A-form should be contained in a small water shell. Cheatham and Kollman showed also that a B$``$A transition can be provoked by placing polycationic ligands to the major groove between the opposite phosphate strands Cheatham:97c . Our results confirm that the accumulation of positively charged counterions in the major groove of B-DNA is, perhaps, the major driving force of this transition.

Although one cannot exclude that the A$``$B transitions in different conditions are due to radically different mechanisms, it is attractive to consider the possibility that the electrostatic model outlined above in fact operates in all such cases. This model is obviously transferable from water drop to DNA in fibers under different humidity. One can argue, tough, that the computed A-DNA systematically differs from the canonical and single crystal A-DNA fragments by a very narrow major groove. However, the major groove in the A-DNA conformations available in NDBNDB:  varies so strongly that this parameter can hardly be considered as a well defined attribute of the A-form. It is possible also that the parameterization of Na+ ions in the Cornell et al. force fieldCornell:95  underestimates its binding to water with respect to negatively charged atoms. In our dynamics, correct canonical major groove width was observed during transitions as long as Na+ ions kept their first hydration shells and avoided direct contacts with the phosphate groups.

Another defect of the computed A-DNA consists in the strong B-philicity of the termini accompanied by bending to the major groove. These two features are mutually related because it is long known that the boundaries between A and B-forms are bentSelsing:79 . This defect, however, is only apparent because it qualitatively agrees with experimental data. The B-philicity of DNA ends has to be assumed when A$``$B transitions in solution are consideredMinchenkova:86 . Moreover, many single crystal A-DNA structures appear smoothly bend to the major groove when analyzed with the Curves programCurves: , with this trend being particularly evident for complete helical turns Malinina:99 ; Verdaguer:91 ; Bingman:92 . In our simulations, these features are just accentuated due to the excessively narrow major groove, which requires stronger bending, and also because our sequences are not A-philic as in all single crystal A-DNA fragments.

The solution A$``$B transitions caused by non-polar organic solvents can also be interpreted in terms of the above electrostatic model. For this purpose it is sufficient to imagine what would happen if our water drops were put in EtOH, for instance. The water activity inside the drop is strongly reduced by DNA itself and its counterion shell. In infinite water, the activity coefficient grows with the distance from the double helix and becomes close to one outside the counterion cloud, that is at around 10 Å from DNA. EtOH is easily mixable with normal water, but not with the low activity water near DNA. This argument can be illustrated by the above discussed crystallization of NaCl in the DNA environment at concentrations much below its standard solubility limit. The water molecules would leave the drop and go to EtOH until the outer water activity becomes equal to that inside the drop. In other words, the water activity in the added nonpolar phase determines the size of the drop and, consequently, the state of the A$``$B equilibrium. That is why B$``$A transition curves in different solvents converge when water activity is used as parameterMalenkov:75 . More polar solvents like MetOH and ethylene glycol do not cause B$``$A transitions possibly because they solvate counterions rather than push them inside the double helix.

Both the apparent B-philicity of the DNA ends and the chain length dependence of the B$``$A transitions observed in our simulations agree with experiments in solutionMinchenkova:86 ; Fairall:89 ; Galat:90 . Another salient feature nicely reproduced is the cooperative character of the transition, with the minimal length of the cooperative fragment of around one helical turn. This cooperativity is not caused by the sugar pseudorotation barrier because during B$``$A transitions the sugars switch last. Instead it is conditioned by the A-DNA geometry and the character of the electrostatic interactions in the major groove. In the A-DNA major groove the phosphate of the nth base pair is positioned against that of the (n-8)th. These eight base pairs represent the minimum DNA length that is to be mechanically involved if the transition is caused by a force applied between these groups. The positive energy of this mechanical deformation should be compensated by the negative electrostatic energy of solvent cations sandwiched between the phosphates. One step of such ”electrostatic sandwich” is probably insufficient, and octamer DNA should not generally go to the A-form. However, already in a nanomer the sandwich length is increased by one while the deformation energy is increased only by a fraction of the initial gap. Therefore, the energetic balance is progressively shifted in favor of A-form as the chain length is increased. In long DNA, the boundary between the A and B-forms has higher energy because the energy gap of the initial octamer remains effectively shared between the two opposite boundaries. The foregoing qualitative pattern exactly corresponds to the one dimensional Ising model of cooperative transitionsIvanov:74 .


text:
# References

Renormalization group analysis
of the $`\tau `$-lepton decay within QCD***Electronic version of the published paper, the journal reference: Yad. Fiz. 54, 1114 (1991) and
Z. Phys. C53, 461 (1992).

A.A. Pivovarov
Institute for Nuclear Research of the Russian Academy of Sciences
Moscow 117312, Russia

Abstract

A technique to sum up the regular corrections appearing under the analytic continuation from the spacelike momentum region to the timelike one is proposed. A perturbative part of the inclusive semileptonic decay width of the $`\tau `$-lepton in analyzed in detail.

At present the theory of electroweak interactions is tested with truly great accuracy due to creation of the modern experimental machines  and to the happy possibility to fulfill the reliable theoretical calculations within the perturbative framework in the gauge coupling constants which are small . The standard model as a whole however can not be checked with the same impressive accuracy because of uncertainties introduced by the strong interactions into the theoretical predictions at low energies.

Inclusive process $`\tau \nu +hadrons`$ is one of the most convenient channels to explore strong interactions in the low energy region and, in particular, to determine the numerical value of the strong coupling constant $`\overline{\alpha }_s`$. The reason consists in two crucial observations: first, the problems of describing the final hadronic states disappear while calculating the full width of the decay and, second, the nonperturbative corrections to the decay rate are small due to some specific features of the process kinematics .

In the present note we estimate the accuracy of the calculation for the perturbative part of the full decay width $`\mathrm{\Gamma }(\tau \nu +hadrons)`$ when the analytic continuation effects are summed up in all orders of the perturbative $`\overline{\alpha }_s`$-expansion.

An expression for the full width of the decay $`\tau \nu +hadrons`$ can be written down in the form

$$\mathrm{\Gamma }(\tau ^{})=\mathrm{\Gamma }(\tau ^{}\nu _\tau e^{}\overline{\nu _e})(4.9728+12r_H)$$
(1)
where the perturbative contribution of the strong interactions $`r_H`$ which we are interested in is given by the following spectral integral

$$r_H=2_0^{m_\tau ^2}r(s)\left(1\frac{s}{m_\tau ^2}\right)^2\left(1+2\frac{s}{m_\tau ^2}\right)\frac{ds}{m_\tau ^2}.$$
(2)
The spectral density $`r(s)`$ is determined by the correlator of weak hadronic currents (see refs.  for more details) and is connected with the theoretically calculable quantity $`d(Q^2)`$  by the Källen-Lehmann representation

$$d(Q^2)=Q^2\frac{r(s)ds}{(s+Q^2)^2}.$$
(3)
The perturbative expansion of the function $`d(Q^2)`$ in the running strong interactions coupling constant $`h(Q^2)=\overline{\alpha }_s(Q^2)/4\pi `$ looks like

$$d(Q^2)=h(Q^2)(1+a_1h(Q^2)+a_2h(Q^2)^2+a_3h(Q^2)^3+\mathrm{})$$
(4)
where the coefficients $`a_1,a_2,a_3`$ are determined by the divergent parts of the three, four, five loop Feynman diagrams correspondingly.

To calculate the full decay width through eqs. (1) and (2) it is necessary to find the moments of the spectral density $`r(s)`$

$$r_n=m_{\tau }^{}{}_{}{}^{2(n+1)}_0^{m_\tau ^2}r(s)s^n𝑑s.$$

An accurate definition of such moments can be given with the help of the integration contour which goes along the positive semi axis coinciding with the cut in the complex $`Q^2`$ plane and closes itself far enough from all physical singularities . The result for these moments does not depend on the explicit way of the closing of the contour if it lies in the analyticity region of the function $`d(Q^2)`$. Choosing the circle of the radius $`m_\tau ^2`$ for the closing of the contour we get the representation for the moments in the following form

$$r_n=m_{\tau }^{}{}_{}{}^{2(n+1)}_0^{m_\tau ^2}r(s)s^n𝑑s=\frac{1}{2\pi }_\pi ^\pi p(m_\tau ^2e^{i\varphi })e^{i(n+1)\varphi }𝑑\varphi $$
(5)
where the function $`p(Q^2)`$ is given by the expression

$$p(Q^2)=p(\mu ^2)_{\mu ^2}^{Q^2}\frac{d\xi }{\xi }d(\xi ).$$
(6)
The running strong coupling constant $`h(Q^2)`$ satisfies the standard renormalization group equation

$$Q^2\frac{d}{dQ^2}h(Q^2)=\beta (h(Q^2))=\beta _0h(Q^2)^2\beta _1h(Q^2)^3\beta _2h(Q^2)^4$$
(7)
and the coefficients of the expansion for the $`\beta `$-function are equal to $`\beta _0=9`$, $`\beta _1=64`$, $`\beta _2=3863/6`$ in the $`\overline{\mathrm{MS}}`$-renormalization scheme for $`N_F=3`$ .

At the contour $`Q^2=m_\tau ^2e^{i\varphi }`$, $`\pi <\varphi <\pi `$ the running strong coupling constant $`h(Q^2)`$ becomes the complex function of the angle $`\varphi `$ and the renormalization group equation takes the form

$$\frac{h(m_\tau ^2e^{i\varphi })}{\varphi }=i\beta (h(m_\tau ^2e^{i\varphi })).$$
(8)
Setting $`h(Q^2)=h(m_\tau ^2)(u+iv)^1`$ at the contour one gets for the functions $`u(\varphi )`$ and $`v(\varphi )`$ the system of equations in the given approximation (7) for the $`\beta `$-function

$`{\displaystyle \frac{u}{\varphi }}`$ $`=`$ $`\beta _1h^2{\displaystyle \frac{v}{u^2+v^2}}+\beta _2h^3{\displaystyle \frac{2uv}{(u^2+v^2)^2}}`$
$`{\displaystyle \frac{v}{\varphi }}`$ $`=`$ $`\beta _0h+\beta _1h^2{\displaystyle \frac{u}{u^2+v^2}}+\beta _2h^3{\displaystyle \frac{u^2v^2}{(u^2+v^2)^2}}`$
with the initial conditions $`u(0)=1`$, $`v(0)=0`$.

Formulas (5,6,S0.Ex2) allow to sum up all the regular corrections in $`h(m_\tau ^2)`$ to the moments $`r_n`$ of the spectral density $`r(s)`$ with the renormalization group technique. The naive expression for the moments $`r_n`$ can be obtained by expanding the exact formulas in the series in $`h(m_\tau ^2)`$ and can be also found by the following simple trick. The renormalization group properties of the quantity $`r(s)`$ up to the order $`𝒪(h(\mu ^2)^3)`$ can be explicitly expressed in the form

$$r(s)=h(\mu ^2)+h(\mu ^2)^2\left(a_1+\beta _0L\right)+h(\mu ^2)^3\left(a_2+(\beta _1+2\beta _0a_1)L+\beta _0^2\left(L^2\frac{\pi ^2}{3}\right)\right)$$
(10)
where $`L=\mathrm{ln}(\mu ^2/s)`$ and we have taken into account the representation (3) and expression (4). After integrating eq. (10) over the interval $`(0,m_\tau ^2)`$ (more precisely the corresponding function $`p(Q^2)`$ along the contour in the complex plane of the variable $`Q^2`$) we get formally the following expression for the zero moment, for example, at $`\mu ^2=m_\tau ^2`$

$$r_0=h+h^2(a_1+\beta _0)+h^3\left(a_2+\left(\beta _1+2\beta _0a_1\right)+\beta _0^2\left(2\frac{\pi ^2}{3}\right)\right),$$

where we use the simple notation $`h=h(\mu ^2)`$.

It is clear that such an integration generates a good deal of regular terms connected with the coefficients of the $`\beta `$-function which can and must be summed up with the renormalization group technique .

A following remark is in order here. Within the naive approach, eq. (10), it is not clear how many terms of the expansion in $`h(\mu ^2)`$ must be kept. In ref. , for example, all terms of the same order as in the expansion for the $`d(Q^2)`$-function were kept in the $`h`$-expansion for the spectral density $`r(s)`$ in the $`\overline{\mathrm{MS}}`$-renormalization scheme. This prescription does not work if the so called effective renormalization scheme is used for the function $`d(Q^2)`$ in which the corrections are absent at all due to definition ($`a_1=0`$, $`a_2=0`$, $`a_i=0`$ for $`i>2`$. Within our approach all regular terms determined by the expansion of the $`\beta `$-function are summed up and the action instruction is fully independent of the choice of the renormalization scheme for the quantity $`d(Q^2)`$. After integration by parts in eq. (5) we get the representation for the moments

$$r_n=\frac{1}{n+1}(r_f+\mathrm{\Delta }_n),$$

where

$$r_f=\frac{p(m_\tau ^2+i0)p(m_\tau ^2i0)}{2\pi i}$$

$$=\frac{1}{\beta _0\pi }\left(\mathrm{arctan}(v/u)+k_1h\frac{v}{u^2+v^2}+k_2h^2\frac{uv}{(u^2+v^2)^2}\right),$$
(11)
and $`k_1=a_1c_1`$, $`k_2=a_2c_2+c_1^2a_1c_1`$, $`c_1=\beta _1/\beta _0`$, $`c_2=\beta _2/\beta _0`$, the variables $`u`$ and $`v`$ in formula (11) denote the values of the solution to the system of equations (S0.Ex2) at the point $`\varphi =\pi `$, while the quantities $`\mathrm{\Delta }_n`$ are determined by the integrals

$$\mathrm{\Delta }_n=\frac{1}{2\pi }_\pi ^\pi d(m_\tau ^2e^{i\varphi })e^{i(n+1)\varphi }𝑑\varphi .$$
(12)
Substituting the solution of the system of equations (S0.Ex2) into formulas (11-12) one gets after performing the integration the numerical values for the moments $`r_n`$ of the spectral density $`r(s)`$ with the full summation of the analytic continuation effects by the renormalization group method. The experimental value obtained for the quantity $`r_H`$ from the relation (1) after averaging over several possible ways of extracting the numerical value for $`\mathrm{\Gamma }(\tau \nu +hadrons)`$ and $`\mathrm{\Gamma }(\tau ^{})=\mathrm{\Gamma }(\tau ^{}\nu _\tau e^{}\overline{\nu _e})`$ is equal to

$$12r_H=0.50\pm 0.07.$$
(13)
The theoretical prediction is based on formula (2) and is expressed through the moments of the spectral density $`r(s)`$ as follows

$$r_H=2(r_03r_2+2r_3).$$
(14)
The naive expression for the quantity $`r_H`$ up to the order $`O(h^3)`$ has the form

$$r_H^N=h_N+h_N^2\left(a_1+\frac{19}{12}\beta _0\right)+h_N^3\left(a_2+\left(\frac{19}{12}\beta _1+2\beta _0a_1\right)+\beta _0^2\left(\frac{265}{72}\frac{\pi ^2}{3}\right)\right),$$
(15)
where $`h_N`$ denotes the running coupling constant $`h`$ extracted from the naive formula (15) and the coefficients $`a_1`$, $`a_2`$ are $`a_1=6.564`$ , $`a_2=102.0`$  in the $`\overline{\mathrm{MS}}`$-scheme of renormalization.

The results of determination of the running coupling constant $`h`$ and the QCD-parameter $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ in the different approximations in the coupling constant with the help of formulas (14) and (15) are given in Tables 1-3. The mass scale of QCD $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ is determined through the numerical value of the running coupling constant $`h`$ according to the relation

$$\mathrm{ln}(m_\tau ^2/\mathrm{\Lambda }_{\overline{\mathrm{MS}}}^2)=\frac{1}{\beta _0}\left(\frac{1}{h}c_1\mathrm{ln}\frac{1+c_1h}{h\beta _0}\right)+_0^h𝑑\xi \left(\frac{1}{\beta (\xi )}\frac{1}{\beta _{(2)}(\xi )}\right)$$

where $`\beta _{(2)}(h)=\beta _0h^2\beta _1h^3`$.

It is seen from the results given in the Tables that the difference between the exact value (summed up in all orders in $`h`$) and the naive one becomes smaller with the improvement of the approximation i.e. after taking into account the larger number of corrections in $`h`$. It is explained by the obvious observation that the leading term of the difference between the exact result and the naive one has formally the next order in $`h`$ in comparison with the kept terms (cf. the situation with the usual renormalization group summation of the perturbative corrections). Since the quantity $`h`$ (more precisely $`r_H`$ (13)) is small enough the estimate of the analytic continuation effects by means of the expansion of exact formulas (5,14) into the series in $`h`$ is valid numerically. For the value $`12r_H=1`$, however, the situation is essentially changed even numerically as it can be seen from the last lines of the tables although the estimates remain obviously unchanged parametrically. Even for the moderate numerical values of $`r_H`$, for example, $`12r_H=0.5`$, however, the change of the parameters $`h`$ and $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ under taking into consideration the next approximation in $`h`$ is comparable in the magnitude with the changes influenced by the summation of the analytic continuation effects. The accuracy of experimental determination of the quantity $`r_H`$ is permanently getting improved and after accomplishing the construction of and putting into operation the $`c\tau `$-factories  will be considerably better. So, if $`r_H=\overline{r}_H+\mathrm{\Delta }r_H`$ then at $`\overline{r}_H=0.60`$ and $`\mathrm{\Delta }r_H=0.03`$ one can reliably tell apart the exact and naive results for $`h`$ and $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$.

In conclusion I would like to stress once again the main result of the paper. The knowledge of the analytic and renormalization group properties of the function $`d(Q^2)`$ allows to sum up all the regular contributions into the moments of the spectral density $`r(s)`$ in all orders of perturbation theory in the strong coupling constant $`h=\overline{\alpha }_s(m_\tau )/4\pi `$. The summation of such contributions at the calculation of the inclusive semileptonic width of the $`\tau `$-lepton is necessary for precise comparison of the parameters $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ extracted from the different processes depending on different energy scales.

Table 1. The leading approximation $`a_1=0`$, $`a_2=0`$, $`\beta _1=0`$, $`\beta _2=0`$.

$`\begin{array}{ccccc}& & & & \\ & & & & \\ 12r_H^{\mathrm{exp}}\hfill & 10^2h^{(1)}\hfill & 10^2h_N^{(1)}\hfill & \mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{(1)}\hfill & \mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{N(1)}\hfill \\ & & & & \\ & & & & \\ 0.30\hfill & 1.97\hfill & 1.94\hfill & 107\hfill & 102\hfill \\ 0.40\hfill & 2.51\hfill & 2.44\hfill & 196\hfill & 183\hfill \\ 0.50\hfill & 3.04\hfill & 2.90\hfill & 287\hfill & 262\hfill \\ 0.60\hfill & 3.55\hfill & 3.32\hfill & 374\hfill & 334\hfill \\ 0.70\hfill & 4.08\hfill & 3.71\hfill & 457\hfill & 399\hfill \\ 1.00\hfill & 5.77\hfill & 4.76\hfill & 682\hfill & 555\hfill \end{array}`$

Table 2. The next-to-leading approximation $`a_2=0`$, $`\beta _2=0`$.

$`\begin{array}{ccccc}& & & & \\ & & & & \\ 12r_H^{\mathrm{exp}}\hfill & 10^2h^{(2)}\hfill & 10^2h_N^{(2)}\hfill & \mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{(2)}\hfill & \mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{N(2)}\hfill \\ & & & & \\ & & & & \\ 0.30\hfill & 1.74\hfill & 1.72\hfill & 161\hfill & 155\hfill \\ 0.40\hfill & 2.17\hfill & 2.11\hfill & 277\hfill & 261\hfill \\ 0.50\hfill & 2.57\hfill & 2.45\hfill & 390\hfill & 357\hfill \\ 0.60\hfill & 2.96\hfill & 2.75\hfill & 496\hfill & 441\hfill \\ 0.70\hfill & 3.36\hfill & 3.03\hfill & 596\hfill & 515\hfill \\ 1.00\hfill & 4.79\hfill & 3.74\hfill & 876\hfill & 692\hfill \end{array}`$

Table 3. The next-next-to-leading approximation.

$`\begin{array}{ccccc}& & & & \\ & & & & \\ 12r_H^{\mathrm{exp}}\hfill & 10^2h^{(3)}\hfill & 10^2h_N^{(3)}\hfill & \mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{(3)}\hfill & \mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{N(3)}\hfill \\ & & & & \\ & & & & \\ 0.30\hfill & 1.70\hfill & 1.70\hfill & 140\hfill & 140\hfill \\ 0.40\hfill & 2.09\hfill & 2.07\hfill & 237\hfill & 232\hfill \\ 0.50\hfill & 2.45\hfill & 2.39\hfill & 330\hfill & 316\hfill \\ 0.60\hfill & 2.81\hfill & 2.68\hfill & 416\hfill & 387\hfill \\ 0.70\hfill & 3.17\hfill & 2.95\hfill & 498\hfill & 448\hfill \\ 1.00\hfill & 4.69\hfill & 3.62\hfill & 751\hfill & 585\hfill \end{array}`$


text:
# Proximity effect in ultrathin Pb/Ag multilayers within the Cooper limit

## Abstract

We report on transport and tunneling measurements performed on ultra-thin Pb/Ag (strong coupled superconductor/normal metal) multilayers evaporated by quench condensation. The critical temperature and energy gap of the heterostructures oscillate with addition of each layer, demonstrating the validity of the Cooper limit model in the case of multilayers. We observe excellent agreement with a simple theory for samples with layer thickness larger than 30Å. Samples with single layers thinner than 30Å deviate from the Cooper limit theory. We suggest that this is due to the ”inverse proximity effect” where the normal metal electrons improve screening in the superconducting ultrathin layer and thus enhance the critical temperature.

The superconducting proximity effect is a well known phenomenon and has drawn a lot of interest both from the fundamental and the practical points of view. In a high transmission contact between a superconductor and a normal metal the superconducting wave-function varies smoothly across the interface causing a suppression of the pair amplitude in the superconductor and an enhancement of pairing on the normal side. The characteristic distance in which superconductivity “leaks” into the normal region is the normal state coherence length $`\xi _N=\sqrt{\text{h}\text{ }\text{ }D/2\mathrm{\Pi }kT}`$, where D is the diffusivity, and the length-scale at which superconductivity is suppressed in the superconducting side is the superconducting coherence length, $`\xi _S`$.

In the Cooper limit cooper  (in which both the superconductor and the normal metal layers are smaller than the characteristic coherence lengths $`\xi _S`$ and $`\xi _N`$ respectively) the behavior of an S-N bilayer is well described within the framework of the de Gennes model de Gennes ; deutscher . The electrons experience an average paring interaction (average between the two materials) and $`T_C`$ can be described within the BCS weak coupling form adapted to the Cooper limit proximity effect in the following way :

$$kT_C=1.13\text{h}\text{ }\text{ }\omega _De^{\frac{1}{N(0)V}}$$
(1)
$$[N(0)V]_{S+N}=\frac{d_S[N(0)V]_S+dN[N(0)V]_N}{d_S+d_N}$$

here $`\omega _D`$ is the debye frequency, N(0)V is the pairing interaction (N(0) being the density of states at the Fermi level and V the pairing interaction), and dS,N is the thickness of the superconducting and normal metal films respectively. This model assumes that the debye frequencies in the two metals are similar and it is based on the understanding that the electrons spend $`\frac{N_Nd_N}{(N_Nd_N+N_Sd_S)}`$ of their time in the normal metal and $`\frac{N_Sd_S}{(N_Nd_N+N_Sd_S)}`$ in the superconductor ($`N_N`$ and $`N_S`$ are the density of states in the normal metal and in the superconductor respectively).

The above picture is modified in the case of strong coupled superconductors where $`T_C`$ is given by the McMillan expression mcmillen  :

$$T_C=\frac{\omega _D}{1.45}\mathrm{exp}[\frac{1.04(1+\lambda )}{\lambda \mu ^{}(1+0.62\lambda )}]$$
(2)
here $`\lambda `$ is the dimensionless electron-phonon coupling coefficient and $`\mu ^{}`$ is the effective coulomb repulsion. In a normal-metal/strong coupled superconductor within the Cooper limit one can expect $`\lambda `$ in equation 2 to be replaced by an average electron-phonon coupling silbert :

$$\lambda _{S+N}=\frac{\lambda _Sd_S+\lambda _Nd_N}{d_S+d_N}$$
(3)
The proximity effect in the Cooper limit has been observed experimentally in numerous systems of superconductor/normal metal bilayers. In this paper we study systems of multilayers of superconductor/normal metal (Pb/Ag) in which the total thicknesses of the metals are comparable to the relevant coherence lengths. The transport and tunneling experimental results confirm that the basic idea of the Cooper limit notion, i.e. that the electrons experience an average pairing, is valid also in a multilayer sample. We observe good *quantitative* agreement with theory as long as the thickness of the layers is larger than $`30`$Å while thinner layers deviate from the theoretical predictions. We suggest that the simple picture is modified in ultra-thin layers due to effects of screening of e-e interactions (the inverse proximity effect as previously reported inverse ).

The Pb and Ag layers were prepared by quench condensation, i.e. evaporation on a cryogenically cooled substrate within the measurement apparatus strongin . This allows in-situ sequential depositions under UHV conditions and simultaneous transport and tunneling measurements on a single sample. This method has several essential advantages for the study of the superconducting proximity effect prox . Since Pb and Ag are immiscible, the alloying of these two materials is very improbable even at room temperatures. Alloying becomes practicably negligible when the samples are quench condensed and are kept below T=10K throughout the entire experiment (sample growth and measurement). In addition, the extremely clean environment of the experiment leads to very high quality interfaces between the superconducting layers and the normal metal layers. Such clean interfaces provide ideal conditions for proximity effects.

The geometry of the samples used in this study is shown in the insert of figure 1. We begin by evaporating two strips of Al onto a room temperature polished glass substrate. These strips are used both as voltage leads and as the base electrode for tunneling measurements. Then we allow a native oxide barrier to grow on the Al for 20 minutes. We connect metallic leads to the substrate and place it in an evaporation chamber which is then pumped out and immersed in liquid helium. A 10Å Ge adhesion layer is then quench condensed across the Al strips, permitting the evaporation of a continuous layer of Pb as thin as 5Å strongin . The thicknesses are monitored by a calibrated quartz crystal situated in the chamber. Multilayers, made of base units of Pb/Ag bilayers having different thicknesses, are quench condensed across the Al strips. This configuration allows the simultaneous measurements of transport properties of the multilayers (thus determining the critical temperature), and tunneling into the multilayers (for determination of the energy gap). We performed these measurements at incremental stages of the evaporation of each of the Pb or Ag layers. All measurement were performed in a screened room using standard AC methods.

In this paper we present data on three Pb/Ag systems having different base unit thickness: Sample #1 where $`d_{pb}=d_{Ag}`$32Å, sample #2 where $`d_{pb}=d_{Ag}`$15Å, and sample #3 where $`d_{pb}=d_{Ag}`$9Å. Each of these unit bilayers is within the cooper limit; $`\xi _N`$of the Ag at low temperatures is of the order of 1 $`\mu m`$ and $`\xi _S`$ of the amorphous Pb (which is in the dirty limit) at low temperature is given by $`\xi ^{}=\sqrt{\xi _0l}`$ where $`\xi _0`$ is the clean limit coherence length (for Pb this is 800Å) and *l* is the mean free path which is approximately an interatomic distance, about 2.7Å in our samples remark2 . These values yield a zero temperature coherence length, $`\xi ^{}46`$Å, which is larger than each one of our Pb layers. Furthermore, as the temperature increases the coherence length increases until it diverges at $`T=T_C.`$

We begin by presenting the results for the initial bilayer. Figure 1 shows the resistance versus temperature of sample #2 as a layer of 14.5 Å thick Ag is grown in steps of $``$1.55Å on top of a 14.4Å thick initial Pb layer. Such measurements are used to extract $`T_C`$ by taking the midpoint of the transition (RN/2) remark1 . It is seen that the initial Pb layer has a $`T_C`$ which is smaller than the bulk value of 7.2K. This is typical for very thin superconducting layers which are known to exhibit a superconductor - insulator transition as a function of film thickness (and Rsq) graybeal ; white ; goldman . It has been argued uniform ; lee  that as Rsq (the disorder) increases, the electronic screening is reduced. Coulomb interactions are then enhanced leading to a weakening of the superconducting coupling. This explains the decrease of $`T_C`$ in the Pb films when the thickness is reduced from 30Å to 9Å as is seen in figure 2. The enhancement of Coulomb interactions for high Rsq is also evident in figure 3 which depicts the dI/dV-V measurements (proportional to the density of states) of a 9Å thick Pb film with increasing layers of Ag, taken at temperatures above $`T_C`$. It is seen that in the ultrathin Pb layer the density of states at the Fermi level is relatively low. The tunneling conductance exhibits a strong voltage (energy) dependence demonstrating that this is a strongly coulomb correlated system. The density of states dramatically increases when layers of Ag are evaporated on top of the Pb, presumably due to enhancement of screening and, hence, suppression of the Coulomb interactions uniform .

As subsequent layers of Ag are added (and Rsq is reduced as $`1/d_{Ag}`$ ) $`T_C`$ is decreased even further due to the proximity effect as evident in figure 1. Since the bilayer is in the cooper limit, $`T_C`$ is indeed expected to decrease as a function of the normal metal thickness. We note, however, that this is not always the case. We have shown inverse  that when the initial Pb layer is extremely thin (thinner than that of sample 3) the first added Ag sub-layers cause $`T_C`$ to increase in what seems to be an inverse proximity effect . This behavior was interpreted as being due to the fact that the Ag layers enhance screening in the system thus suppressing the e-e interactions which are responsible for the $`T_C`$ suppression. This trend opposes the usual proximity effect and becomes more important in thinner superconducting films.

The evolution of $`T_C`$ as more Pb and Ag layers are added is shown in figure 2 for the three samples. It is seen that $`T_C`$ oscillates when we change between Pb and Ag depositions reaching a local minimum or maximum at the completion of each layer. The oscillation amplitude decreases with increasing number of bilayers and TC appears to approach a value of $`3K`$ in all three samples. The oscillation behavior is also observed in the energy gap measurements which are illustrated in figure 4 for sample #1. The dI/dV curves were measured at T=1.65K well below the $`T_C`$ of the multilayer. Note that despite that fact that tunneling is performed into the bottom Pb layer, the results are still affected by adding the 6th layer on top. This clearly demonstrates that within the cooper limit the tunneling electrons probe the *entire* sample and the superconducting parameters are determined by the mean properties of all the layers in the multilayer. Figure 4 also shows the ratio $`\frac{2e\mathrm{\Delta }}{kT_C}`$ as a function of the number of layers. The superconducting gap was evaluated from a classical temperature dependent BCS fit. It is seen that this ratio crosses over from a value of 4.7 (which is close to the strong coupling limit of $``$4.8) to 3.6 (which is close to the week coupling limit of $``$3.5). The fact that this ratio is relatively constant as we add the last few layers gives us further confidence that all our multilayers are within the Cooper limit.

We have attempted to fit the curves in figure 2 to the proximity effect modified McMillan expression of equation 2 using an average $`\lambda `$ for the Cooper limit (equation 3) in the spirit of the de-Gennes model. These are shown in figure 2 were we used the known values for Pb: $`\lambda _S=1.55,\omega _D=105K,\mu ^{}=0.11`$ and taking a low electron-phonon coupling for Ag, $`\lambda _N=0.2,`$ without any fitting parameters. It is seen that the experimental results are in excellent agreement with this model for sample #1 (where $`d_{pb}=d_{Ag}`$32Å). We note that the values for $`\lambda ,\omega _D`$ and $`\mu ^{}`$ are not adjusted, but accepted values and so the agreement is quite impressive. However, when we used the same parameters for the other samples with thinner base layers, we get substantial deviations from the model for the first few layers (figure 2). In fact, in sample #1 (and in the theory) the oscillations are superimposed on a global decrease of $`T_C,`$ in sample #2 the background is roughly constant and in sample #3 TC globally increases with the number of bilayers. We argue that this behavior is an extension of the inverse proximity effect observed in the ultrathin bilayers inverse . The large Coulomb interactions characteristics of these thin samples are continuously weakened with increasing numbers of Pb/Ag bilayers, thus, as more layers are added, we approach the Cooper limit proximity effect theory.

In summary, we have shown that the de Gennes considerations regarding the Cooper limit are valid for the case of multilayers until the layers become very thin. Using a simple model we were able to describe very accurately the variation of the critical temperature and the energy gap as a function of the number of layers. We showed that for very thin layers the model breaks down indicating that other physics is involved. We suggest that this breakdown is due to the effect of strong Coulomb interactions in the superconducting ultrathin films. These interactions are screened out as normal metal is added to the superconductor resulting in an inverse proximity effect which must also be taken into account in the thinnest layers.

We gratefully thank O. Naaman and W. Teizer for fruitful discussions. The work was supported by the Israel-US Binational Foundation #1999332 and by the NSF grant #DMR0097242.


text:
# Self-interaction effects in (Ga,Mn)As and (Ga,Mn)N

## I Introduction

The discovery of ferromagnetism in Mn-doped GaAs and InAs ohno1  has changed conventional thinking about electronics. In fact these diluted magnetic semiconductors (DMSs) combine the functionalities of semiconductors with those of magnetic materials, paving the way for a complete integration of data storage and logic on the same electronic component wolf . Unfortunately the most widely studied DMS, (Ga,Mn)As, shows a Curie temperature ($`T_C`$) far below room temperature and therefore appears to be unsuitable for commercial applications. This has led to a considerable world-wide effort in searching for new DMSs with higher $`T_C`$s.

Part of this research is guided by a pioneering theoretical work by Dietl and collaborators, who predicted remarkably high $`T_C`$s for wide-gap magnetic semiconductors at moderate Mn doping dietl . Their model is based on the Zener mechanism for ferromagnetism, according to which a conduction electron can hop across two magnetic atoms through a double exchange of place with the valence electrons of the bridging anions. Spin-conservation implies that this conduction can only occur between spin-aligned magnetic centers, thus a ferromagnetic ordering of the Mn ions is expected. In the model, $`T_C`$ is crucially dependent on hole and Mn concentration, and on the strength of their mutual coupling. High $`T_C`$ ferromagnetism in wide-gap DMSs has recently been confirmed by several groups who reported the existence of above-room-temperature ferromagnetism in Ga1-xMnxN gamnn1 ; gamnn2 ; gamnn3 ; gamnn4 . However, in contrast to “conventional” As-based DMSs, where the existence of double-exchange induced ferromagnetism is clear, the origin and the nature of the observed magnetism in Ga1-xMnxN are controversial. On one hand, both $`p`$\- and $`n`$-doped samples are, apparently, ferromagnetic, thus indicating that hole-mediated double-exchange may not be, in fact, the driving mechanism. More importantly, recent magnetic circular dichroism (MCD) measurements ando  on high $`T_C`$ samples show that the MCD signal in the vicinity of the GaN band-gap comes from a paramagnetic phase. This suggests the possibility of two phases in (Ga,Mn)N: a paramagnetic phase, attributed to substitutional Mn at the Ga site, and dominated by $`s,p`$-$`d`$ interactions between Mn and the host GaN, and an unidentified ferromagnetic phase, not detected by MCD, which might be due to small particles that are invisible to diffraction methodsando .

In summary, there are sound indications that (Ga,Mn)As and (Ga,Mn)N may not be, in fact, very similar as far as magnetic ordering is concerned. Thus a careful examination of their microscopic properties is necessary in order to give a correct interpretation of the experimental results.

With this goal in mind, in this paper we investigate the electronic structures of zincblende (Ga,Mn)As and wurtzite (Ga,Mn)N using first-principles density functional theory (DFT) calculations. Band energy calculations for these compounds obtained within the familiar local-spin density functional theory (LSDA) are already present in the literature.mer ; me ; kalotov ; fong  Here we present the results obtained with our recently developed pseudopotential self-interaction free approach (pseudo-SIC).fh\_sic

It is known that the band energies calculated within LSDA have systematic errors caused by the presence in the one-electron Kohn-Sham (KS) potential of its self-interaction (SI)si  (that is the interaction of an electron with the potential generated by the electron itself). The SI is especially relevant in ionic and mixed ionic-covalent materials. Indeed the poor LSDA description of the electronic properties of II-VI semiconductors vogel , transition metal oxides tmo , and magnetic perovskites pero  is a consequence of the SI fh\_sic . In general, any compound characterized by localized electron charges is, in principle, affected by a non-vanishing, unphysical, SI contribution.

Regarding the compounds investigated in this paper, it is known that the LSDA band structure of bulk GaN shows severe discrepancies with photoemission results (for example the fundamental energy gap is severely underestimated and the position of the $`d`$ bands is incorrect). Furthermore, the LSDA is generally inaccurate in describing the electronic properties of oxides of manganese, and in particular fails to give the correct position of the Mn d states with respect to the oxygen-derived valence manifold.

Since the band gap of the host material and the position of the transition metal d states are crucial in the description of these diluted magnetic semiconductors, there is the necessity to confront the LSDA description with alternative approaches. Our pseudo-SIC is ideally suited for this task since it repairs to a large extent the LSDA failures and provides a reliable description of the electronic properties for a large range of materials, including GaN and Mn-based perovskites. In addition, and in contrast with other beyond-LSDA approaches, our method preserves the conceptual simplicity and the computational feasibility of the LSDA, and thus the possibility to approach large-sized systems such as the diluted magnetic semiconductors.

In this paper we determine the extent to which correction of the SI changes the electronic properties and the chemical pictures of wurtzite (Ga,Mn)N and zincblende (Ga,Mn)As. We find that in (Ga,Mn)As the localization of the Mn 3d electrons is weak because both on-site and off-site hybridizations occur, and, as a consequence, the SI is not strong. Therefore the LSDA already provides a good description of (Ga,Mn)As and the pseudo-SIC does not make large qualitative changes. In (Ga,Mn)N, however, the SI is much stronger, due to the large N electronegativity, and consequent higher ionicity of GaN. Therefore the LSDA description is inadequate and the pseudo-SIC gives qualitatively different results. The difference in localization between the two compounds in turn has a profound effect on their electronic properties.

This paper is organized as follows: in section II we describe the methodology and the technicalities used for the calculations. Sections III and IV present the results for (Ga,Mn)As, and (Ga,Mn)N, respectively. Finally, in sections V and VI we present our conclusions.

## II Methodology

In this paper we use our recently developed self-interaction free local-density scheme fh\_sic  based on ultrasoft pseudopotentialsuspp  and a plane wave basis. The scheme is described in detail in Ref. fh\_sic ; here we summarize the key points. The Kohn-Sham equations are modified in order to effectively extract from the one-electron potential the spurious self-interaction contribution:

$`\left[^2+\widehat{V}_{\text{pp}}+\widehat{V}_{\text{hxc}}^\sigma \widehat{V}_{\text{sic}}^\sigma \right]|\psi _{n𝐤}^\sigma =ϵ_{n𝐤}^\sigma |\psi _{n𝐤}^\sigma ,`$ (1)
where $`\widehat{V}_{\text{pp}}`$ is the usual ion-core pseudopotential projector, $`V_{\text{hxc}}^\sigma `$ the LSDA-KS potential, and $`\widehat{V}_{\text{sic}}^\sigma `$ the SIC term. The latter is cast as a nonlocal Kleinman-Bylander-type pseudopotential projector:

$`\widehat{V}_{\text{sic}}^\sigma ={\displaystyle \underset{i}{}}{\displaystyle \frac{|V_{\text{hxc}}^{\sigma ,i}\varphi _ip_i^\sigma \varphi _iV_{\text{hxc}}^{\sigma ,i}|}{\varphi _i|V_{\text{hxc}}^{\sigma ,i}|\varphi _i}}.`$ (2)
Here $`\varphi _i(𝐫)`$ is the pseudowavefunction, $`V_{\text{hxc}}^{\sigma ,i}[\varphi _i^2(𝐫)]`$ the SI potential of the $`i^{th}`$ atomic orbital ($`i`$ runs over angular quantum numbers and atomic positions, $`\sigma `$ over the spin), and $`p_i^\sigma `$ is the corresponding occupation number calculated self-consistently by projecting the Bloch state manifold onto the basis of atomic pseudowavefunctions. This procedure is no more computationally demanding than the usual LSDA, and, in contrast with other SIC schemes vogel , it can be applied to both metals and insulators, a feature which is crucial for the present application.

In earlier work our method was applied to a number of different systems including III-V and II-VI ionic semiconductors (e.g. GaN and ZnO),fh\_sic  transition metal monoxides (MnO, NiO)fh\_sic  and perovskites (CaMnO3, BaTiO3, YMnO3)bec\_sic . In all cases the pseudo-SIC obtained both the correct physical structures, and band structures in excellent agreement with photoemission experiments. This series of successes lends strong support to the credibility of the present results.

We use ultrasoft pseudopotentials to allow moderate cut-off energies (35 Ryd and 40 Ryd for (Ga,Mn)As and (Ga,Mn)N, respectively). For the pseudo-SIC calculations the Ga 3d electrons are treated as core states (this choice is fully justified in Ref. fh\_sic .) For zincblende (Ga,Mn)As and wurtzite (Ga,Mn)N 16-atom and 32-atom supercells are used respectively, in each case with one Mn ion substituting for one Ga ion, corresponding to doping levels of $``$ 12% and $``$ 6%. The lattice parameters are fixed to the values calculated for the host materials (a = 10.68 bohr for GaAs, a = 6.03 bohr, c= 9.802 bohr, u= 0.377 for GaN), which are in good agreement with the experimental data.

## III Zincblende (Ga,Mn)As

### III.1 LSDA

The calculated LSDA band structure of (Ga,Mn)As (in figure 1) is that of an half-metal with metallicity in the majority spin band. Here the most important feature for the electronic and magnetic behavior is given by the three bands crossing the Fermi energy in the majority manifold. They are rather dispersed and span $``$2 eV between $`\mathrm{\Gamma }`$ and M. An orbital decomposition of these bands shows that they consist of Mn d $`t_2`$ character, as well as Mn and As $`p`$, with significant contributions from both nearest and next-nearest neighbor As ions. This indicates that part of the charge is extended beyond the Mn-centered tetrahedron. The next set of bands lower in energy corresponds to the valence band top of bulk GaAs at the $`\mathrm{\Gamma }`$ point ($`E=0.73`$eV) but mixes in Mn character elsewhere in the Brillouin Zone, as does the doublet below it. An orbital decomposition shows that the Mn d character is not restricted to a single set of bands close to the Fermi level, but instead is widely distributed through the GaAs sp valence manifold.

In the upper part of Table 1 we report the occupation numbers of individual orbitals. The total magnetization per cell is $``$ 4 $`\mu _B`$, almost entirely due to the Mn (3.72 $`\mu _B`$). The Mn d electron charge is 5.48, of which 4.61 and 0.87 comes from spin-up and spin-down electrons, respectively.

In conclusion, according to the LSDA description (and in agreement with previous calculations mer ; me ) (Ga,Mn)As is dominated by Mn d - As p hybridization rather than on-site electron correlation. However, since the overestimation of p-d hybridization and the suppression of on-site Coulomb energy are general tendencies of the LSDA, a comparison with the pseudo-SIC is necessary in order to assess the reliability of this picture.

### III.2 pseudo-SIC

In Figure 2 we show the band structure calculated within pseudo-SIC. It is immediately apparent that the main features are not changed significantly with respect to the LSDA, apart for the expected increase of the fundamental energy gap of the bulk GaAs due to the downward shift of the occupied valence band manifold with respect to the empty bands. Also the band manifold occupies an energy region $``$ 1 eV wider than that in the LSDA result. Interestingly, the three majority, half-occupied bands have an even larger dispersion than in LSDA, but the shape of the bands is very similar.

In the lower part of Table 1 we report the orbital occupations calculated within pseudo-SIC. The total magnetization is 4 $`\mu _B`$, as in LSDA. However, as a consequence of a larger Mn d spin splitting, the Mn magnetic moment (4.58 $`\mu _B`$) is larger than in LSDA, but it is compensated by a non discardable spin-polarization on the As (0.17 $`\mu _B`$) which is antiparallel to the Mn moment. (Table 1 reports the values for a single orbital (px) which has a magnetization of -0.06 $`\mu _B`$ in pseudo-SIC. The total magnetization for px \+ py \+ pz is -0.17 $`\mu _B`$.)

In conclusion, the pseudo-SIC result is a sound confirmation that the p-d hybridization and the weakly-correlated character of the Mn d electrons is in fact a reliable prediction of the LSDA, and validates the extensive earlier computational literature on (Ga,Mn)As. (for a review see Ref.mer ).

## IV Wurtzite (Ga,Mn)N

### IV.1 LSDA

In the hexagonal crystal symmetry of the wurtzite structure, the d orbitals are split into two doublets ((dxy, d$`_{x^2y^2}`$), and (dxz, dyz)) and one singlet (d$`{}_{}{}^{2}{}_{z}{}^{}`$), where $`x`$ and $`y`$ lie in the hexagonal plane and $`z`$ is parallel to the c axis. Similarly, the p orbitals are split into a doublet (px, py) and a singlet (pz). In Figure 3 we report the band structure of wurtzite (Ga,Mn)N calculated within the LSDA. Like (Ga,Mn)As, this material is half-metallic and has three partially filled bands close to the Fermi energy. However this three-band manifold (consisting of a doublet composed of Mn dxy-d$`_{x^2y^2}`$ and the px-py orbitals of the three planar N, and a singlet derived from Mn d$`_{z^2}`$-pz and on-top N pz) is much narrower in (Ga,Mn)N than in (Ga,Mn)As. The doublet below these states ($`E`$ = –1.39 eV) at $`\mathrm{\Gamma }`$ is also a p-d combination (of Mn dxy-d$`_{x^2y^2}`$, and N px) with a strongly localized electron charge confined almost entirely to orbitals internal to the Mn-centered tetrahedron. The small band dispersion (less than 1 eV) is evidence of the strongly-correlated character. The same features can be appreciated in the orbital-resolved density of states, reported in Figure 4 (only majority spin components are shown). The hybridization between Mn d and the p orbitals of the four surrounding N is clearly visible, while almost no contribution comes from second neighbor N.

The orbital occupation numbers for (Ga,Mn)N are given in the upper part of Table 2. The Mn magnetic moment is 3.63 $`\mu _B`$ with a total Mn d electron charge of 5.3. Notice how the much stronger electronegativity of N with respect to As is reflected in the comparison between the orbital occupations in (Ga,Mn)As (Table 1) where roughly 60% and 30% of the s, p electron charge is distributed to As and Ga, and in (Ga,Mn)N, where the proportion is 70% to N and 20% to Ga.

These band structures agree qualitatively with earlier LSDA calculationsleeor ; mark , although the exact position of the Mn d levels differs according to the details of the pseudopotentials or the lattice parameters used in the calculations.

### IV.2 pseudo-SIC

Strong electron localization makes the LSDA results for (Ga,Mn)N potentially affected by a large self-interaction. Our calculations indeed confirm these expectations. The pseudo-SIC band structure and the orbital-resolved density of states are shown in Figures 5 and 6), respectively.

Due to the self-interaction, the fully-occupied dxz, and dyz bands are pushed down in energy, well below the p-like valence band top of the host material. Around the Fermi energy two electrons occupy the other three majority Mn d bands, similarly to what we have seen within LSDA. However now these three bands are closer to the top of the p band manifold, and, due to the large increase of the GaN energy gap, much farther from the bottom of the s-like conduction bands.

In the lower part of Table 2 the occupation numbers are shown. The local magnetic moments are 4.18 $`\mu _B`$ for Mn, 0.06 $`\mu _B`$ for the on-top N, and –0.10 $`\mu _B`$ for each of the three planar N. This means that the coupling between the on-top N pz orbital and the majority d$`_{z^2}`$ orbital is ferromagnetic, whereas the coupling between planar N $`p_x`$ and py orbitals with dxy and d$`_{x^2y^2}`$ orbitals is antiparralel, thus contributing negatively to the total magnetic moment M=4 $`\mu _B`$.

This remarkable increase of Mn magnetic moment with respect to the LSDA values is a typical self-interaction correction effect: the spatial localization is increased and the hybridization reduced.

Finally, notice that these results are not very sensitive to the Mn concentration: band energy calculations of wurtzite (Ga,Mn)N with 12% Mn doping within both LSDA and pseudo-SIC shows no relevant changes with respect to the corresponding calculations at 6% concentration.

## V Discussion

The chemical picture of (Ga,Mn)N obtained within the pseudo-SIC (which we believe to be the most accurate representation) is markedly different both from the LSDA description of (Ga,Mn)N, and also the behavior of (Ga,Mn)As. The difference between (Ga,Mn)N and (Ga,Mn)As can be traced back to the large electronegativity of N, which in turn favors ionicity and charge transfer over hybridization. This aspect, already evident within LSDA, is further emphasized within pseudo-SIC, due to its tendency to favor electronic configurations with completely filled or empty orbitals.note

The differences are so large that they will certainly lead to different charge mobilities in the two materials, and affect the nature of the magnetic ordering as well. Indeed the band alignment obtained within pseudo-SIC for (Ga,Mn)N may not be appropriate for validating the usual Zener picture for ferromagnetism. In fact, according to our calculation, there are no free holes in the GaN valence band. Instead, the holes are localized near the Mn impurity, and occupy a triplet of 1 eV wide bands around the Fermi energy. This could lead to two possible scenarios. On the one hand substitutional Mn in GaN could be responsible only for the paramagnetic component of the magnetization as indicated by MCD measurementsando . On the other hand this could be indicative of the formation of a Zhang-Rice polaron in (Ga,Mn)N as has been recently suggested ZR . According to this model, Mn doping in GaN behaves as an effective mass acceptor, with a configuration of d5 plus a hole localized on the magnetic ion. This local singlet could move through the sublattice of Mn2+ ions and cause ferromagnetic ordering (still mediated by a double-exchange mechanism) at much higher Tc then that of the usual Zener “free hole” double exchange. Our calculations show a magnetic moment for Mn of 4.18 $`\mu _B`$, which is suggestive of a d4 configuration. However the presence of very flat, partially occupied bands with d character allow us to re-interpret our results in terms of a Mn d5 configuration plus a tightly bound hole of d-character. These are the conditions for the formation of a Zhang-Rice polaron. To further assess the validity of these hypotheses, larger-sized calculations, capable of describing the properties of these compounds in different magnetic orderings (e.g. paramagnetic and antiferromagnetic) will be necessary in the future.

## VI Conclusions

In summary, we have compared the electronic structure of ferromagnetic (Ga,Mn)As and (Ga,Mn)N within standard LSDA and within our pseudo-SIC approach which sets the KS potential free of the spurious self-interaction part. Our results indicate two main conclusions. First, that the LSDA picture of (Ga,Mn)As as a half-metal, with the Fermi energy crossing three weakly correlated bands which derive from Mn d - As p hybridized states, is indeed a good description, largely unchanged when the self-interaction is corrected. However for (Ga,Mn)N the LSDA description is inadequate because the d electrons are strongly correlated, and thus the pseudo-SIC description is markedly different from that of the LSDA. Second, that the behavior of (Ga,Mn)N is qualitatively different from that of (Ga,Mn)As. Three flat bands of dxy, d$`_{x^2y^2}`$, and d$`_{z^2}`$ character are partially occupied by two spin-polarized electrons, confined within the Mn-centered tetrahedron, thus effectively describing a highly confined hole around the Fermi energy which may be more consistent with the formation of a Zhang-Rice polaron, rather than with the free-hole mediated double-exchange mechanism of the Zener model.

S.S. thanks Enterprise Ireland (grant SC-2002/10) for financial support. A.F. and N.A.S.’s development of the pseudo-SIC method has been supported by the NSF under grant number DMR-0312407, and DMSs work by the ONR under grant number N00014-00-1-0557.