text:
# Ab-initio study of model guanine assemblies: the role of 𝜋-𝜋 coupling and band transport

## I Introduction.

The twofold issue of scaling down the electronic devices and of realizing a high circuit integration density on a single chip has recently seen the rise of the field of molecular electronics, which consists of the use of molecules to realize electrically conducting structures rev-aviramratner ; rev-gimzewski . Because of their sequence-specific recognition properties, DNA molecules are attracting attention for the construction of nanometer scale devices, where they might be used, by virtue of their self-assembling capabilities, to wire the electronic materials in a programmable way braun-nat98 . This research path has led recently to a set of controlled experiments for the direct measurement of the d.c. conductivity.

Using interdigital electrodes, anisotropic conductivity was found in an aligned DNA cast film: At room temperature, a large ohmic current, linearly increasing with the applied voltage, was measured okahata-jacs98 . Ohmic behavior and high conductivity were found, also, for a a 600 nm long rope made f a few $`\lambda `$-DNA molecules fink-nat99 . Instead, nonlinear current/voltage curves, exhibiting a voltage gap at low applied voltage, were measured through a single 10 nm long poly(dG)/poly(dC) DNA molecule trapped between two metal nanoelectrodes porath-nat00 . Large currents were observed, in air and in vacuum, both at ambient temperature and at 4 K. The authors suggested that the observed electron transport is best explained by a semiconductor-like band model where the electronic states are delocalized over the entire length of the base pair stack porath-nat00 .

The guanine (G) base is particularly interesting in view of obtaining conductive molecular aggregates, because of its low ionization potential, which suggests its viability to mediate charge motion along a sequence of bases meggers ; jortner98 . This peculiar property has opened the way to the measurement of the the electrical conductivity of G aggregates. Interestingly, a metallic nanogate filled with a dried solution of a lipophilic derivative of 2’-deoxyguanosine apl ; prl , displayed a current/voltage behavior note-counterions  similar to that of the semiconducting poly(G)/poly(C) sample porath-nat00 . Previous investigations had shown that, in organic solvents, this molecule undergoes extensive self-assembly, mediated by H-bonding among the guanine bases, to give ribbon-like aggregates that, upon drying, form fiber structures within which the guanine cores of the ribbons lie on parallel planes at a distance of about 3.4 Å gottarelli98 ; gottarelli00 . While a direct comparison among the findings for the deoxyguanosine materials and for DNA cannot be made, because of the different overall experimental settings and chemical nature of the molecules, one cannot avoid to notice the similar and peculiar current/voltage characteristics of the nanogates interconnected by molecules featuring self-assembled apl ; prl  or inherent porath-nat00  guanine stacks.

In the following, we have chosen to perform ab-initio calculations of the structural and ground state electronic properties of extended model structures whose building block is the G base alone. Of course, we are well aware that these model structures are only partially related to the structure of a poly(G)/poly(C) duplex porath-nat00 , where each guanine is H-bonded to a cytosine in the opposite strand. However, self-assembled G structures are among the simplest base aggregates characterized experimentally prl ; gottarelli98 , thus allowing accurate theoretical calculations as well as comparison with experimental data. Additionally, our choice is motivated by the role played by guanine, the DNA base with the lowest ionization potential, in the mechanism of charge transport bixon-hhop99 ; schuster-polaron ; GGenhancedtransport .

The main question that we address here is whether the electronic properties of extended G-based structures can account for a band-like mechanism for charge transport. Previous theoretical studies for DNA bases and for their assemblies can be found in the existing literature sponer1 ; sponer2 ; sponer3 ; sugiyama-saito ; machado ; parrinello . Both MP2/Hartree-Fock and density-functional-based calculations were performed in the past with localized basis sets (GAUSSIAN) to treat isolated guanines and small clusters of guanines sponer1 . Here, we focus on density-functional-based calculations using a plane wave basis and ab-initio pseudopotentials fhi-dft : this technique should allow for a correct description of solid aggregates of G’s, resembling those of long-range deoxyguanosine fibers, provided the single molecules and the clusters are well described. Before simulating the structures of our interest, which are planar ribbons and stacked sequences of G’s, we check that our technique is well suitable to describe isolated G molecules as well as pairs of molecules in different (lateral or vertical) configurations. We then consider (GG)n vertical stacks, and (GG)n isolated and stacked periodic ribbons. We study the stability of different vertical configurations as a function of the relative atomic positions between adjacent G molecules, and identify the role of $`\pi `$-$`\pi `$ coupling. We demonstrate that pseudopotential plane-wave Density Functional Theory (DFT) calculations at the current level of accuracy are able to reproduce not only the equilibrium length of isolated hydrogen bonds, but also periodic sequences of such bonds. Concerning the electronic properties, all the model solids that we considered are semiconducting with large energy gaps. The dispersion of the highest valence and lowest conduction bands is always negligible in the $`(x,y)`$ plane containing the molecules (the guanines are connected laterally by H bonds). Conversely, the dispersion along the $`z`$ direction perpendicular to the G planes found to be extremely sensitive to the detailed geometric stacking of the bases, which in turn can be affected by different environments. In vertical stacking geometries that are similar to those present in DNA, the calculated dispersion of both the valence band maximum (VBM) and conduction band minimum (CBM) are extremely small. Consequently, both electrons and holes introduced by “doping” or photoexcitation would have very large effective masses – hence low mobility – in these types of structures. Instead, our results show that a band-like contribution may be responsible for the conduction mechanism in the 2’-deoxyguanosine lipophilic derivative prl .

The paper is organized as follows. Section II describes the computational method. Section III deals with the calculated equilibrium geometries and the electronic structure for the isolated G molecule and for the model guanine assemblies. Section IV presents a discussion of our results with an outlook to the implications for the physics of guanine ribbons and DNA-like stackings. Finally, section V contains a summary of the arguments presented in the paper.

## II Method.

Our calculations are based on the DFT in the Local Density Approximation (LDA) dg-dft . For the hydrogen-bonded pairs and ribbons, we take into account BLYP gradient corrections to the exchange-correlation functional gga-blyp . The electron-ion interaction is described via ab-initio norm conserving pseudopotentials in the factorized form of Kleinman and Bylander pseudo . The search for optimized metastable structures is performed by total energy minimizations with respect to the ionic and electronic degrees of freedom. The former are represented by the ionic coordinates, the latter by the electronic wavefunctions. The ions are treated in a classical formalism, and are displaced according to the forces derived from the potential determined by the full quantum mechanical electronic structure fhi-dft , within a Car-Parrinello-like scheme cp-scheme . For any selected geometry, all the atoms are allowed to relax, until the forces vanish within an accuracy of 0.05 eV/Å. Thus, for each metastable structure, we obtain both the geometry and the consistent single-particle electron energies and wavefunctions. The electron wavefunctions are expanded in a basis of plane waves with kinetic energy up to 50 Ry. This cutoff is very high with respect to standard first-principle calculations for most solids, and is due to the presence of the first-row elements C, N, and O, whose valence electron wavefunctions have strong oscillations in the region around the nucleus, thus needing many plane waves for an accurate treatment. For one of the model systems, we have verified that increasing the precision of the calculations both in the plane-wave kinetic energy cutoff (up to 60 Ry) and in the accuracy within which the atomic forces vanish (0.025 eV/Å), gives changes of the bond lengths smaller than 0.01 Å and of the bond angles smaller than 0.5 degrees, within a guanine molecule, whereas the inter-planar distance in the stacks is not affected at all. Additionally, for the diatomic molecules N2 and O2 the employed pseudopotential was tested up to a cutoff of 80 Ry, with no significant improvement.

Our method employs periodic boundary conditions in three dimensions. In order to simulate isolated molecules, we choose a large supercell of size 15.9 Å $`\times `$ 15.9 Å $`\times `$ 10.6 Å. Such a choice ensures that the minimum distance (in any spatial direction) between two molecules is larger than 8.5 Å. In particular, the distance is 10.6 Å in the direction perpendicular to the plane of the molecules, very large with respect to the distance of 3.37 Å between two neighboring bases in B-DNA. The supercells for different assemblies of G molecules are described in section III, where the model structures are discussed in detail. For Brillouin Zone sums, the single high-symmetry $`\mathrm{\Gamma }`$ point has been employed in the case of isolated molecules and dimers, while one or two (depending on the symmetries) special k points monk-pack  have been employed in the case of periodic columns and ribbons, and in the case of stacked ribbons.

The computational technique has been successfully applied in many investigations of the structure and electronic properties of inorganic and organic materials dg-dft ; review-dft-inorganic ; review-dft-molec . In section III, before reporting the results of our calculations for the model guanine assemblies, we address the issue of extending such calculations to biomolecules, by presenting a test on the G base.

## III Results.

In this section we present our results about the isolated G molecules and about their assemblies and discuss them in the frame of the existing theoretical literature, which is limited to isolated G’s and G-pairs sponer2 ; sponer3 ; sugiyama-saito  and does not give any account of periodic G columns and ribbons. We relate the outcome of the calculations to new experimental findings about the structure and electrical behavior of deoxyguanosine-based solids prl .

The section is divided into sub-sections for the different systems: (A) the isolated G molecule; (B) the stacked GG dimers; (C) the stacked (GG)n columns; (D) the hydrogen-bonded GG pairs; (E) the hydrogen-bonded (GG)n ribbons and the stacked ribbons.

### III.1 The Isolated Guanine Molecule.

We start by presenting the calculated structure and selected electron states for the G molecule. This simple system allows us to evaluate the accuracy of our method, by comparison with X-ray data and with the outcome of previous ab-initio calculations. Furthermore, the total energy of the equilibrium structure of the isolated molecules is necessary in order to evaluate formation energies of dimers, ribbons, and columns.

The structure of guanine is well known from X-ray studies saenger ; guaninexp . Figure 1(a,b) shows a plane view of the molecule, and the isosurface plot for the total charge density. The calculated bond lengths and angles are in good agreement with X-ray data (see Table 1). Most of the bond lengths are underestimated within 2% and 3% with respect to the experimental structure, the only exceptions being the underestimate of 4.3% for the C6-O6 bond and the overestimate of 1.6% for the C6-N1 bond. The average C-H and N-H distances, not reported in Table 1, are 1.08 Å and 1.01 Å, respectively, in good agreement with bond lengths in NH3 and CH3 molecules. The bond angles in guanine are also reproduced with a high degree of accuracy (Table 1): the discrepancy with respect to the experimental data is below 2.5%, but the average percentage error is 1%. The formation energy of G, with respect to its elemental components in stable phases (O2, H2, and N2 molecules, crystalline diamond) is 87.5 Kcal/mole. By starting the atomic relaxation with a planar G molecule as the initial condition, we find that the planar configuration is indeed a metastable state. By considering a different initial condition with a non-planar amino group (the NH2 complex bonded to the site C2, see Figure 1a), we find another metastable state. However, the deviation from planarity is small, and the total energy difference between the planar and the puckered geometries is also very small, within the precision of our calculations (estimated to be about 10 meV/G note-eform ). While previous calculations have pointed out a stronger stabilization effect of non-planarity sponer3 , no direct gas phase experimental data exist. Moreover, the hydrogen bonds tend to flatten the structure when forming dimers and ribbons of G’s; thus, we are confident that our results for the periodic structures are not affected by this issue.

The single particle eigenvalues are characterized by a DFT-LDA energy gap of 4.8 eV between the highest occupied (HOMO) and the lowest unoccupied (LUMO) electron states dg-dft . We have identified $`\sigma `$ and $`\pi `$ orbitals. The HOMO (see Figure 1c) has a $`\pi `$ character and is localized on the C8, O6, N2 atoms, and on the C4-C5 and N3-C2 bonds. The LUMO (see Figure 1d) has a $`\pi `$ character and is localized on the N9, C4, C5, C2, and N1 atoms. Because they extend out of the guanine plane, both the HOMO and the LUMO states are well suitable for interactions with adjacent similar states when forming vertical stacks, inducing a splitting of degenerate molecular orbitals by an amount that depends on the strength of such interactions. In the case of an infinite periodic stack, this is the mechanism that may give rise to band dispersion for sufficiently strong coupling, and to mobile carriers if the bands are partially occupied (for instance, as a consequence of doping or photoexcitation).

### III.2 Stacked GG dimers.

In order to select low-energy geometries for the vertical columnar structures, whose electronic properties are the main subject under investigation, we have first considered stacked dimers. We name stacked dimer a pair of G molecules lying in parallel planes whose distance in the perpendicular direction is an output of our calculations. A column is obtained by periodically replicating a dimer along the stacking (perpendicular) direction.

We have analyzed several configurations (Figure 2) characterized by the relative azimuthal rotation angle of the two G’s in the pair, with respect to an axis perpendicular to the G plane, and by an in-plane translation. The supercell used in the calculations was 15.9 Å $`\times `$ 15.9 Å $`\times `$ 19.1 Å. With periodic boundary conditions, two dimers in neighboring supercells are 15.7 Å apart in the stacking direction, sufficient to avoid spurious interactions between them. By allowing all the atoms to relax, the interplanar distance between the two G’s in a pair was a free parameter. In this way we determined the equilibrium interplanar distance to be used in the calculations for the periodic columns. For all the configurations shown in Figure 2 (e.g., independently of the relative rotation angle between the bases), the average interplanar distance was 3.37 Å, typical of base stacks in B-DNA. Each base maintained a planar geometry, with out-of-plane fluctuations smaller than 0.05 Å.

Top views of the computed stacked dimers are illustrated in Figure 2. In Figure 2a, the two G’s of the dimer are perfectly eclipsed, in particular the hexagonal rings are on top of each other, and there is maximum superposition of the $`\pi `$-like HOMO and LUMO of the two guanines (label GGv.A). In Figure 2b, the rotation angle is zero as for GGv.A, but the center of mass of one molecule is shifted with respect to the other, in order to lower the $`\pi `$-$`\pi `$ superposition (label GGv.B). In Figure 2c, the azimuthal rotation angle is zero, as well as the relative translation of the two G’s, but there is a reflection of the upper molecule with respect to its plane: this configuration (label GGv.C) is the only one, among those considered in this work, that exhibits a reflection and allows to discriminate between two molecular faces. In Figure 2d, the rotation angle is 180, and the $`\pi `$-$`\pi `$ superposition is large, though smaller than in GGv.A (label GGv.D). In Figure 2e, the azimuthal rotation angle is 36, and a translation brings the two molecules in a configuration similar to that in B-DNA (label GGv.E).

Although some of these stacked structures are not likely to occur in nature, they allow to understand important microscopic features that may be relevant in real structures. Most notably, the dependence of the stability and of the electronic properties on the azimuthal angle and $`\pi `$ superposition is accessible. The lowest-energy configuration among those shown in Figure 2 is GGv.D, while GGv.A has the highest formation energy, and the other dimers have intermediate formation energies note-faces . However, the energy difference between the two extreme cases GGv.A and GGv.D is small, about 250 meV/G. We attribute the highest formation energy of GGv.A to the electrostatic repulsion due to the complete $`\pi `$-$`\pi `$ superposition. In fact, the $`\pi `$-like HOMO’s of guanines in neighboring planes are mostly responsible for their interaction: the superposition of negative charge in the same region of space (the hexagonal ring) for configuration GGv.A contributes with a Coulomb repulsion. The $`\pi `$-$`\pi `$ superposition is large also in GGv.D, but the repulsion is much smaller, making GGv.D a more viable model for a stacked GG pair sponer2 .

We wish to point out that, although the lower superposition of adjacent $`\pi `$ orbitals of GGv.D with respect to GGv.A decreases the electrostatic repulsion between the two molecules and makes the dimer more viable, such interaction is still rather large and is the origin of energy dispersion in columnar structures based on the GGv.D dimer. This issue is addressed in the next sub-section, by computing the band structures of model periodic stacks built up with the dimers just described. It is also worth stressing that, for the stacked dimers, without inclusion of sugars, phosphates, and water as in real situations, the rotation angle of 36 is not preferred. This result is in line with the fact that the inner core of the base pair stack is not completeley responsible for the stability of double-stranded DNA, but also the backbone and the environment are relevant factors.

### III.3 Periodic stacked poly(G) columns.

The building block of each periodic stack is a GG pair; we have analyzed five different model columns starting from the dimers described above (Figure 2). The supercell for these calculations is 15.9 Å $`\times `$ 15.9 Å $`\times `$ 6.74 Å: no vacuum region separates two adjacent dimers in the stacking direction. The periodic structures are labeled with the same names as the stacked dimers. Direct and reciprocal one-dimensional crystal lattices are associated to the periodic columns: the basis vectors of these lattices are a3=6.74 Å, and b3=$`\frac{2\pi }{a_3}\widehat{\mathrm{a}}_3`$.

The formation energies are reported in Table 2. The energetical order is the same as for the dimers. The most stable configuration is GGv.D. We stress again that configuration A, with very large $`\pi `$-$`\pi `$ superposition, has the highest formation energy. The energy is substantially reduced in the column GGv.D, where the superposition is still large but the repulsion is weaker: the hexagonal rings lie on top of each other, but the O atoms in adjacent planes lie opposite to each other. Thus, we attribute the dominant repulsive contribution of the $`\pi `$-$`\pi `$ interaction to the O atoms.

A detailed analysis of the electronic properties reveals other interesting features. We report the numerical data in Table 2 and we show the bandstructure for columns GGv.A, GGv.D, and GGv.E in Figure 3.

GGv.E, with the two G’s rotated by 36 as in portions of nucleic acids, has flat HOMO- and LUMO-derived bands, very large effective masses, and is therefore incompatible with even a partial band transport mechanism for the conduction in guanine-based devices prl . The periodic columns GGv.A and GGv.D, instead, have electronic properties that may support electronic transport through the base stack. For GGv.A, the HOMO-band disperses downwards by 0.65 eV and the LUMO-band disperses upwards by 0.52 eV, between the center ($`\mathrm{\Gamma }`$) and the edge (A=$`\frac{1}{2}`$b3, b3 being the reciprocal lattice vector along the stacking direction) of the Brillouin Zone. For GGv.D, the dispersions of the HOMO- and LUMO-bands are, respectively, 0.26 eV downwards and 0.13 eV upwards. The effective masses of GGv.A and GGv.D, reported in Table 2, are of the order of 1$`÷`$2 free electron masses. Although these values are much larger than those of conducting organic polymers polymers , they are similar to those of inorganic materials for which band transport is demonstrated. Notable examples of such materials are wide-bandgap semiconductors, such as group-III nitrides nitrides . As a note to support the strength of our results, we wish to point out that a test calculation for the model GGv.A demonstrates that the values of the band dispersions do not depend (within 0.02 eV) on the accuracy required for vanishing atomic forces (0.025 eV/Åinstead of 0.05 eV/Å), and to the approximation employed for the exchange-correlation functional (BLYP instead of LDA).

To illustrate how the $`\pi `$-$`\pi `$ stack may originate channels for charge migration through a band transport mechanism, in Figure 4 we show an isosurface of the HOMO state at the A point for the columnar structure GGv.D: the interaction between the two molecules in the cell is evident in the superposition resulting in a delocalized orbital.

Summarizing our study of the columnar structures, we emphasize that the orbital interaction is strong. At variance, we will soon show that it is practically absent in hydrogen-bonded pairs. This is in agreement with the common knowledge that hydrogen bonds have an electrostatic (rather than covalent) nature. In the following subsections, we compare the behavior of planar G pairs and ribbons with that of stacked G pairs and columns. In particular, we find that hydrogen bonding does not give origin to dispersive electron bands, contrary to what we have seen for the $`\pi `$ stacking.

### III.4 Planar hydrogen-bonded GG dimers.

We have selected one possible arrangement of hydrogen bonds, giving a structure named GG3 sponer1 ; machado  (Figure 5a). We have not considered other documented hydrogen-bonded GG pairs sponer1 , because we focus our interest here in the guanine ribbons present in the fibers of a lipophilic derivative of 2’-deoxyguanosine: Such fibers gottarelli00  have the bond network illustrated in Figure 5(b). The equilibrium structure that we obtain after atomic relaxation is in good agreement with that of the GG3 hydrogen-bonded pair previously described via quantum chemistry and DFT cluster calculations sponer1 . The two individual G molecules in the planar pair remain very similar to their isolated form. The N7(H)$`\mathrm{}`$N1 hydrogen bond has a length of 2.93 Å and forms an angle of 177 (2.96 Å and 173 in the theoretical literature sponer1 ). The N2(H)$`\mathrm{}`$O6 hydrogen bond has a length of 2.85 Å and forms an angle of 169 (3.27 Å and 166 in the theoretical literature sponer1 ). Our calculated value for the DFT-BLYP energy gap dg-dft  of the GG3 dimer is 2.45 eV. The formation energy of the structure, with respect to two isolated guanine molecules, is -310 meV/G. This value, compared to the formation energy of the most stable stacked dimer GGv.D ($``$ -200 meV/G), is in agreement with the previous demonstration that hydrogen bonds are stronger than stacking interactions sponer1 .

The structure of the GG3 dimer was calculated only as a basic building block of the fiber-state ribbons, whose electronic properties we are interested in. Thus, we believe that the underestimate of the N2(H)$`\mathrm{}`$O6 distance with respect to the results of other computations is not a serious issue. A slight underestimate of the hydrogen-bonds calculated in the frame of DFT is documented sponer1 . Furthermore, the contraction of the N2(H)$`\mathrm{}`$O6 hydrogen bridge is due to an imprecise account of secondary interactions such as C8(H)$`\mathrm{}`$O6. It is likely that such secondary interaction is reduced in one-dimensional ribbons, thus minimizing the shortcomings of DFT.

### III.5 Planar hydrogen-bonded ribbons.

In a periodic ribbon obtained by piling up replicas of the GG3 dimer (Figure 5b), the equilibrium structure maintains the bonding characteristics of the dimer. There is a huge energy gain in forming the one-dimensional ribbon, of 820 meV/G with respect to isolated G molecules, and of 510 meV/G with respect to the hydrogen-bonded GG3 pair (Table 3). Note that G has a large dipole moment sponer3 ; the dipole moments of G’s add up to give a non-vanishing dipole also in the GG3 dimer and in the ribbon, parallel to the ribbon axis prl . Such electrostatic interactions account for the high stability of the planar hydrogen-bonded structures, well known both in solution and in the solid state.

The bandstructure of the ribbon was calculated along a symmetry line parallel to its axis. The HOMO- and LUMO-derived bands are separated by a DFT-BLYP energy gap dg-dft  of 3.84 eV, and are both dispersionless. Consequently, these bands have practically infinite effective masses, so that the electrons and holes in these states are not mobile according to this picture. The electronic state analysis shows that the HOMO and the LUMO have a $`\pi `$ character, similar to isolated G. The electronic states are localized around single G molecules, no delocalized intermolecular states, extended through the hydrogen bonds in the ribbon, are present. The dispersion induced by hydrogen-bonding is not compatible with band transport. The only possible conduction mechanism through hydrogen bonds is via electrostatic interactions. Therefore, in a device where dried deoxyguanosine fibers are deposited between two metal electrodes prl , if the ribbons are stretched between the electrodes, band transport cannot contribute to conduction.

In order to investigate the competition of $`\pi `$-$`\pi `$ coupling versus hydrogen-bonding interactions, we have simulated two different configurations of stacked ribbons, periodic along the stacking direction. Top views are shown in Figure 5(b,c). Figure 5b shows indeed a single ribbon, but the top view is equivalent for two exactly eclipsed ribbons on top of each other (label SR.A): the stacking between adjacent bases is similar to that of column GGv.A, and all the individual bases are perfectly aligned. Figure 5c shows a stack of two ribbons where only half of the individual G molecules lie on top of each other, with a stacking similar to that of configuration GGv.D (this structure is labeled SR.D). In both the SR.A and SR.D geometries, the projections of the ribbon axes on the $`(x,y)`$ plane coincide, identifying an axis for the stack: this axis defines the $`\mathrm{\Gamma }X`$ direction in the one dimensional BZ.

The supercells are 21.2 Å $`\times `$ 11.3 Å $`\times `$ 6.74 Å for SR.A and 24.3 Å $`\times `$ 11.3 Å $`\times `$ 6.74 Å for SR.D. The choice of different supercells allows to have a vacuum region of the same volume between neighboring ribbons in the two configurations SR.A and SR.D. The relaxed structures maintain the geometry of the single ribbon: no G deformations, no variations of H-bond lengths and angles, and no out-of-plane buckling are observed. This finding indicates that the stacking does not affect the hydrogen-bonding mechanism that determines the structure of the isolated ribbon.

Numerical data for the energetics and the electronic properties are reported in Table 3. Both structures SR.A and SR.D are energetically favorable with respect to isolated G molecules and to stacked columns, but unfavorable with respect to isolated ribbons, as indicated by the smaller energy gain of the stacked ribbons. The configuration SR.D is more stable than SR.A by 280 meV: this trend is consistent with what found previously for the columnar stacks, where a higher stability is achieved by lowering the $`\pi `$-$`\pi `$ superposition in building the stacks. Although the configurations SR.A and SR.D that we have calculated here are not viable models for piling up ribbons into solid-state crystals, ordered phases for deoxyguanosine fibers are known to exist gottarelli98 . We note that other factors, such as presence of the phosphate-sugar backbone, or the possibility of different relative positions of the neighboring ribbons, should be taken into account to achieve a full description, which is beyond the scope of this work.

In Figure 6 we show the calculated bandstructure for the SR.A and SR.D configurations. In both cases, the HOMO- and LUMO-derived bands are dispersionless along the $`\mathrm{\Gamma }X`$ direction, and the charge carriers are not mobile in this direction. This behavior is the same as in isolated ribbons: therefore, it is a further evidence that the base-base interactions due to the stacking do not change the features of the ribbons. In the $`\mathrm{\Gamma }A`$ direction, parallel to the stacking direction, the HOMO- and LUMO-derived bands are dispersive, as in the columnar structures of subsection III.3.

The bands for structure SR.A in the $`\mathrm{\Gamma }A`$ direction are more dispersive than those for structure SR.D, due to the higher $`\pi `$-$`\pi `$ superposition. As a consequence, the corresponding electron and hole effective masses along the $`\mathrm{\Gamma }A`$ direction are smaller. The effective masses of the stacked-ribbon periodic structures are very similar to those of the stacked periodic columns (subsection III.3): this is an indication that the hydrogen-bonding network of the individual ribbons does not affect the electronic properties in the perpendicular direction. The only exception to this trend is me, which becomes infinite for the stacked ribbons SR.D while it is finite for the equivalent stacked column GGv.D: this is due to the fact that in structure SR.D one half of the individual bases are on top of each other (see Figure 5c), while the other half do not participate in the $`\pi `$-$`\pi `$ coupling.

Summarizing the results discussed above, we wish to point out that the $`\pi `$-$`\pi `$ coupling and H-bonding interaction mechanisms are not in competition. H-bonding has an electrostatic nature and accounts for in-plane stabilization: it is not compatible with a band transport mechanism, and it is not modified by base stacking. $`\pi `$-$`\pi `$ coupling is weaker than H-bonding: it is compatible with band transport along the stacking direction and it is not affected by the local details of the planar base sequence.

## IV Discussion.

### IV.1 Guanosine films.

As outlined in the Introduction, the recent lively research activities in molecular electronics have pointed out peculiar guanine assemblies to be exploited as electrical conductors. In fact, the experiments on a lipophilic derivative of 2’-deoxyguanosine prl , demonstrated conduction through the biomolecular material deposited in a nanogate between metal electrodes. The details of the conductivity depend on the on the experimental conditions, in particular on the gate width. Our results indicate that the observed conductivity, resembling that of a semiconductor in the intermediate gate length regime (few hundreds nm) and that of a diode junction in the short gate length regime (less than 100 nm), is compatible with a Bloch contribution to the transport mechanism. We have shown that delocalized orbitals through the base stack may be formed, provided that a relevant superposition of the hexagonal rings of the single guanines is maintained. Therefore, if the guanosine ribbons in the device gate align locally with their plane perpendicular to the direction connecting the electrodes, in a geometry similar to that proposed as SR.A, then extended Bloch-like states may be partially responsible for the observed conduction. This stacking orientation does not need to be complete in order for the proposed mechanism to work: it is sufficient that randomly aligned stacks form locally, and that the total resulting component in the direction connecting the electrodes is non vanishing. In the gate length regime ranging from 100 nm to 300 nm, semiconducting-like conductivity is revealed: we propose that, in such a condition, the ribbons are locally stacked in such a way to form partially delocalized orbitals that give a global band-like contribution. In the gate length regime below a 100 nm, a diode-like behavior is observed: this characteristic is assigned to an interaction of the $`\pi `$ stack with the total dipole moment of the ribbons, as discussed elsewhere prl .

### IV.2 DNA molecules.

While our study is limited to a single type of nucleoside, it also allows a discussion of the electronic properties of DNA molecules, where G sequences in base stacks play an important role grinstaff-rev99 ; barbara-rev99 ; schuster-polaron . Depending on the energetics of the base sequence, and on the overall structural aspects of the system under investigation, the mechanisms proposed for DNA-mediated charge migration include single-step superexchange lewis ; marcus , multistep hole hopping jortner98 ; bixon-hhop99 , phonon-assisted polaron hopping schuster-polaron , and band transport porath-nat00 .

In DNA double strands, the relative arrangement of neighboring bases along the axis of the helix saenger  is characterized by a rotation angle around 36 degrees. Our results show that periodic G columns in such a configuration do not support the formation of extended molecular orbitals. Therefore, band transport would not be effective. For what concerns the mechanism of conductivity in DNA molecules, our results support the conclusion that the contribution of band transport would be very small ye ; bixon-hhop99 ; dePablo , in contrast to recently proposed interpretations of experimental data porath-nat00 , unless structural distortions, possibly activated by temperature effects, may induce rotations that support the formation of partially delocalized electron states. Although in real DNA molecules the backbone phosphate and sugar groups may affect the overall charge mobility, we believe that the presence of the outer mantle would not change our present conclusion that band transport is not supported by the native DNA stack: eventually, it may contribute a hopping or ionic mechanism, but it is unlikely to contribute to the formation of extended electron orbitals. The importance of the base stack for electron transfer through DNA molecules has been highlighted in recent theoretical invesigations ye .

## V Conclusions.

We have reported the results of ab-initio calculations for the structure, energetics, and electronic properties of several guanine assemblies. The geometry of isolated molecules and hydrogen-bonded dimers is well described with our technique, based on plane-wave pseudopotential density functional theory, as demonstrated by comparison to available theoretical and experimental data.

We found that hydrogen-bonding and $`\pi `$-$`\pi `$ coupling are independent mechanisms that control the self-assemblying of guanine bases. Hydrogen-bonding is not responsible for band transport: In fact, we have shown that no band dispersion is present along a planar ribbon of H-bonded guanines. Instead, base stacking is accompanied by $`\pi `$-$`\pi `$ interactions that, for the case of sufficiently large overlap between adjacent $`\pi `$ orbitals (e.g., GGv.D configuration), induce energy dispersion and are consistent with charge mobility. Therefore, band transport may be partially responsible for charge mobility in nucleotide aggregates, in structures characterized by a large base-base superposition. This mechanism is likely complemented by hopping to connect (through space) different regions where such superposition is realized.

## VI Acknowledgements.

We acknowledge the allocation of computer resources from INFM Progetto Calcolo Parallelo. Fruitful discussions with M. Buongiorno Nardelli, R. Cingolani, G. Gottarelli, and R. Rinaldi, are also sincerely acknowledged. We are grateful to F. Grepioni for communicating results prior to publication.


format: text

text:
# MS-TP-01-12 Bicocca-FT-01-21 Quenched twisted mass QCD at small quark masses and in large volumeBased on a poster presented by R. Frezzotti at the XIX International Symposium on Lattice Field Theory ”Lattice 2001”, August 19-24, 2001, Berlin, Germany.

## 1 INTRODUCTION

The lattice twisted mass QCD (tmQCD) action  in infinite space-time volume reads

$`S_\mathrm{W}[U,\overline{\psi },\psi ]=S_\mathrm{g}+a^4{\displaystyle \underset{x}{}}\overline{\psi }_x[D_{\mathrm{tm}}\psi ]_x,`$
$`D_{\mathrm{tm}}=\gamma \stackrel{~}{}\frac{a}{2}^{}+m_0+i\mu _\mathrm{q}\gamma _5\tau ^3,`$ (1)
where $`\stackrel{~}{}=(+^{})/2`$, $`\psi `$ is a flavour quark doublet and the matrix $`\tau ^3`$ acts in the flavour space. If $`m_0m_\mathrm{c}`$, the quark mass is entirely given by the multiplicatively renormalized parameter $`\mu _\mathrm{q}`$.

Lattice tmQCD is a regularization of QCD with $`N_\mathrm{f}=2`$ mass degenerate quarks  where an action term that is proportional to $`\mu _\mathrm{q}`$ and is not aligned in flavour chiral space with the Wilson term contributes to the soft breaking of the chiral symmetry. As a consequence, the lowest eigenvalue of the Hermitean square of the Dirac matrix on any gauge background is bounded from below – in lattice units – by $`(a\mu _\mathrm{q})^2`$. The problem of the exceptional configurations in quenched and partially quenched QCD is thus solved as long as the quark mass is not vanishingly small.

One can think of the two mass-degenerate quark flavours as the $`u`$ and $`d`$ quarks in a chirally twisted basis, with the obvious consequences on the physical meaning of the correlation functions . Heavier flavours of Wilson quarks can still be added e.g. in the usual way, while the neglected tiny mass difference between $`u`$ and $`d`$ quarks can safely be taken into account by means of chiral perturbation theory. Moreover, lattice tmQCD can also be formulated for a doublet of non-degenerate quarks, whilst retaining the protection against exceptionals on all gauge backgrounds.

## 2 OBSERVABLES AND SIMULATIONS

In this study we restrict attention to the pseudoscalar and vector meson masses, $`M_{\mathrm{PS}}`$ and $`M_\mathrm{V}`$, the pseudoscalar decay constant, $`F_{\mathrm{PS}}`$, and the renormalized quark masses, $`\mu _\mathrm{R}`$ and $`m_\mathrm{R}`$, which are defined according to eq. (3.3) and eq. (3.20) of Ref. . It is convenient for our discussion to define polar quark mass ”coordinates” via

$$\mathrm{tan}\alpha =\mu _\mathrm{R}/m_\mathrm{R},M_\mathrm{R}=(\mu _\mathrm{R}^2+m_\mathrm{R}^2)^{1/2}.$$
(2)
We work in the quenched approximation considering systems of physical size $`L^3T`$ at $`\beta =6`$ and $`\beta =6.2`$, with $`L`$ such that $`M_{\mathrm{PS}}L4.5`$ to suppress finite volume effects. Following closely Ref. , we adopt Schrödinger functional (SF) boundary conditions for tmQCD and compute the SF correlators

$`f_\mathrm{A}^{11}(x_0),f_\mathrm{P}^{11}(x_0),f_\mathrm{V}^{12}(x_0),f_1^{11},`$
$`k_\mathrm{V}^{11}(x_0),k_\mathrm{T}^{11}(x_0),k_\mathrm{A}^{12}(x_0)`$ (3)
as a function of the Euclidean time $`x_0`$, for $`L`$ ranging from 1.5 to 2.2 fm, depending on the quark mass values, and $`T/L=`$ 2 or 3. We hence construct the renormalized and O($`a`$) improved SF correlators $`[f_\mathrm{A}^{}^{11}]_\text{R}(x_0)`$, $`[f_\mathrm{P}^{}^{11}]_\text{R}(x_0)`$, $`[k_\mathrm{V}^{}^{11}]_\text{R}(x_0)`$ and $`[k_\mathrm{T}^{}^{11}]_\text{R}(x_0)`$, which are specified in eq. (3.11) of Ref.  with $`\alpha `$ given by eq. (2).

In the limit of large $`x_0`$ and large $`Tx_0`$ and up to cutoff effects, the correlators $`[f_\mathrm{A}^{}^{11}(x_0)]_\text{R}`$ and $`[f_\mathrm{P}^{}^{11}(x_0)]_\text{R}`$, which correspond to the insertion of the operators $`(A_\mathrm{R}^{})_0^1`$ and $`(P_\mathrm{R}^{})^1`$, see eq. (3.10) of Ref. , are expected to be dominated by the lowest isotriplet pseudoscalar state. In the same limit the correlators $`[k_\mathrm{V}^{}^{11}(x_0)]_\text{R}`$ and $`[k_\mathrm{T}^{}^{11}(x_0)]_\text{R}`$, which correspond to the insertion of the operators $`(V_\mathrm{R}^{})_k^1`$ and $`(T_\mathrm{R}^{})_{k0}^1`$, see eq. (3.10) of Ref. , are dominated by the lowest isotriplet vector state.

In our simulations we employ the non-perturbatively improved version of the action (1), with $`c_{\mathrm{sw}}`$ taken from Ref. . An overview of our simulation parameters, statistics and preliminary results is given in Table 1, while $`\alpha =\pi /2+\mathrm{O}(a)`$.

Our most critical simulation (set C) required $`230`$ GFlops $`\times `$ day by employing a CGNE solver for the SSOR-preconditioned Dirac matrix $`D_{\mathrm{tm}}`$. For this (and similar) choice(s) of parameters the BiCG solver often does not converge, while a CGNE solver works fine with small fluctuations in the number of iterations.

## 3 RESULTS

We find that lattice tmQCD allows, as expected, to safely work in a region of parameters which would be inaccessible with ordinary Wilson quarks: see e.g. Fig. 1. For a given number of independent measurements, the statistical errors on $`M_{\mathrm{PS}}`$ and $`M_\mathrm{V}`$ are comparable, up to a factor of one to three, to those found e.g. with domain wall quarks . The CPU time effort, e.g. for the data sets A1, A1’ and A2, is in line with the computational cost for ordinary Wilson quarks.

Concerning the extraction of hadron masses and matrix elements, we recall that the lattice action (1) enjoys reduced parity and isospin symmetries : it just preserves the third isospin generator and parity combined with a flavour exchange. At finite $`a`$ the correlators may hence receive contributions of order $`a\mu _\mathrm{q}`$ from the allowed states with different parity and isospin. While deferring details to Ref. , we show in Fig. 2 an example of effective masses extracted from SF correlators: the correlator $`[f_\mathrm{A}^{}^{11}]_\mathrm{R}(x_0)`$ receives contributions of order $`a\mu _\mathrm{q}`$ that are peculiar to our SF setup and would be absent at $`\alpha =0`$. Analogous effects in $`[f_\mathrm{P}^{}^{11}]_\mathrm{R}(x_0)`$ are expected to be very small .

As detailed in Ref. , we expect the relations among our observables and the renormalized parameters $`r_0`$ and $`M_\mathrm{R}`$ to be O($`a`$) improved. In particular, when working at $`\alpha =\pi /2+\mathrm{O}(a)`$, an O($`a`$) improved estimate of $`F_{\mathrm{PS}}`$ is obtained in terms of $`M_{\mathrm{PS}}`$ and the matrix element of e.g. the operator $`2M_\mathrm{R}(P_\mathrm{R}^{})^1`$ between the charged pion and the vacuum. This approach to $`F_{\mathrm{PS}}`$, which is also valid in unquenched lattice tmQCD, requires to know no other counterterms than $`Z_\mathrm{g}`$, $`m_\mathrm{c}`$ and $`c_{\mathrm{sw}}`$.

In order to check for the residual scaling violations, we produced data at $`\beta =6.2`$ (sets A2 and B2), while keeping $`\alpha `$, $`M_\mathrm{R}`$ and $`r_0`$ fixed. More precisely, the ratio $`m_\mathrm{R}/M_\mathrm{R}`$ was tuned111The small mismatch in $`M_\mathrm{R}r_0`$ for the set A2 was corrected by employing estimates of the dependence of our observables on $`M_\mathrm{R}r_0`$. to the value $`0.000(1)`$ for $`M_\mathrm{R}r_0=0.195(1)`$ and $`0.015(1)`$ for $`M_\mathrm{R}r_0=0.273(2)`$. For the simulation C we found $`m_\mathrm{R}/M_\mathrm{R}=0.083(5)`$. We also reanalysed the data of Ref. , which were produced at $`\beta =6,6.1,6.2,6.45`$, by imposing precisely the same renormalization conditions as in this study of tmQCD. We then performed a continuum extrapolation, assuming a purely quadratic dependence on $`(a/r_0)^2`$ and discarding the data at $`\beta =6`$. However, owing to the use of the one-loop estimate of $`b_\mathrm{A}`$, the resulting estimate of $`F_{\mathrm{PS}}`$ is not fully O($`a`$) improved. The outcome of this exercise is compared with the results from tmQCD in Fig. 3. We omit the case of $`M_\mathrm{V}r_0`$, as cutoff effects are hardly visible within statistical errors.

## 4 CONCLUSIONS

Lattice tmQCD solves the problem of exceptional configurations and allows for rather precise results with a moderate computational effort and tiny scaling violations, at least for the observables of this study. The framework is suitable for extension beyond the quenched approximation.

This work is part of the ALPHA Collaboration research programme.


text:
# Braking and bouncing

## 1. Introduction

The commonly accepted model for the formation of cataclysmic variables (CVs) with hydrogen–rich main-sequence or brown dwarf donor stars involves four main assumptions.

1. The progenitor system was a wide enough binary to allow the progenitor star of the white dwarf, an intermediate–mass star, to evolve to giant dimensions prior to Roche–lobe contact.
2. The system evolved through a common envelope (CE) phase which exposes the giant’s degenerate core, the future white dwarf, and tightens the orbit.
3. Orbital angular momentum (AM) losses drive the post-common envelope system into contact, and drive mass transfer in the semi–detached system.

There are various clues to the magnitude of the AM loss rate.

The overwhelming majority of CVs below 2 hrs orbital period, $`P`$, are dwarf novae, while at longer periods novalikes and dwarf novae are observed with roughly equal frequency. Dwarf novae are unstable disc accretors where the mass transfer rate $`\dot{M}`$ is below a critical limit $`\dot{M}_{\mathrm{crit}}`$ so that the hydrogen–ionisation instability is not suppressed. The relative distribution of novalike CVs and dwarf novae with orbital period thereferore suggests that the transfer rate $`\dot{M}`$ is significantly smaller than $`\dot{M}_{\mathrm{crit}}`$ at short $`P`$, while $`\dot{M}`$ must be close to $`\dot{M}_{\mathrm{crit}}`$ at longer $`P`$. As the critical rate typically scales as $`\dot{M}_{\mathrm{crit}}P^2`$ (e.g. Shafter 1992, and references therein), the secular mean transfer rate, and hence the AM loss rate, must be significantly smaller at short orbital periods ($`P<2`$ hr) than at longer periods.

The well known disrupted orbital braking model for the orbital period gap (Spruit & Ritter 1983, Rappaport, Verbunt & Joss 1993) is an extreme manifestation of this dependence: the AM loss rate drops essentially discontinuously by more than an order of magnitude, causing the systems to detach at the upper edge of the gap. The width of the gap requires that the mass transfer rate at this point is $`\dot{M}12\times 10^9`$ M$`_{}`$ yr-1 for unevolved Pop I donor stars (e.g. Stehle, Ritter & Kolb 1996). This is almost two orders of magnitude larger than the rate driven by gravitational radiation (GR). Conversely, the short–period cut–off of the observed period distribution at 78 min, when interpreted as the minimum period that coincides with the donor’s transition from a main–sequence star to a brown dwarf (Paczyński 1981), seems to require a braking rate which is at most a factor of 3 above the value set by GR.

Constraints on $`\dot{M}`$ from estimated absolute magnitudes of CVs suffer from distance uncertainties, difficulties in determining the duty cycle of dwarf novae, and in mapping visual magnitudes to disc accretion rates. The deduced transfer rates are, on average, larger above the gap than below the gap, and there seems to be a considerable intrinsic spread of the $`\dot{M}`$ values at least for systems above the gap (e.g. Figure 9.8 in Warner 1995).

The spectral type of CV donor stars above the gap is, on average, later than that of an unevolved main–sequence star in thermal equilibrium that would fill its Roche lobe at the same orbital period (e.g. Baraffe & Kolb 2000). In the period range 3–6 hrs this is a natural consequence if the stars experience mass loss at a rate consistent with the period gap width.

Townsley and Bildsten (2002) considered the accretional heating of white dwarfs. They assert that the observationally deduced white dwarf temperatures are consistent with a mean accretion rate of order a few $`\times 10^9M_{}`$ yr-1 in CVs above the gap, and $`<10^{10}M_{}`$ yr-1 below the gap.

The last three clues — gap width, late spectral types, accretional heating — constrain the long–term mean of the transfer rate, averaged over typically a thermal time of the donor star.

In the following three sections I briefly discuss the status of our understanding of magnetic braking, the effect of nova outbursts on the observed mass transfer rates, and the shape of the period distribution below 2 hrs. A recent review that covers also other applications of the standard disrupted braking model of CV evolution is given by Kolb (2001).

## 2. Magnetic braking

It has been suggested that magnetic stellar wind braking, usually thought to be the dominant AM loss mechanism in CVs, is much weaker than previously thought (Andronov, Pinsonneault & Sills 2001; Pinsonneault, these proceedings). This work represents a long overdue generalisation of the so–called Skumanich law for solar–type stars to rapidly rotating stars (in some cases almost as rapid as CV secondaries) and spectral types other than G. The AM loss rate $`\dot{J}`$ has a Skumanich form $`\dot{J}\omega ^3`$ at small values of the rotational angular speed $`\omega `$, but appears to saturate for rapid rotators to give $`\dot{J}\omega `$.

A note of caution applies to this “new” braking law. Building on the work of Krishnamurti et al. (1997), Sills et al. (2000), and others, Andronov et al. deduced this AM loss rate from the observed rotation rates in young open clusters. The procedure involves the deconvolution of various physical mechanisms that all affect the surface rotation of cluster stars: the pre–main sequence contraction of stars, the redistribution of angular momentum inside the star, interaction with protostellar discs, and the rotational braking by magnetic braking. It is clear, though, that the fastest rotators in young open clusters do not brake (by far) as fast as the original Skumanich law would have suggested.

It is equally clear that the AM loss rate deduced by Andronov et al. cannot account for the large accretion rates observed above the period gap, nor for the degree of thermal disequilibrium the spectral type evidence seems to suggest. Likewise, the period gap model would have to be modified because the braking rate is continuous over the fully convective boundary. Until there is a viable alternative to the current thinking that the gap arises because systems are detached while they cross the gap region, and unless all novalike systems represent, in evolutionary terms, short–lived states with a high accretion rate, the AM loss rate in CVs above 3 hrs orbital period must be much larger than the AM loss rate Andronov et al. derive for single stars.

Potential mechanisms for additional orbital braking include circumbinary discs (Spruit & Taam 2001) and accretion disc winds (Livio & Pringle 1994), the latter only in systems close to mass transfer instability (King & Kolb 1995). However, such alternative scenarios usually fail to reproduce sufficiently sharp period gap boundaries. It will be interesting to watch the evolution of the observed period distribution; the most recent list of systems (Downes et al. 2001, period as of 1 July 2001) displays gap boundaries which are less sharply defined than in earlier samples (e.g., compare Figure 1, left panel, with Fig. 1 in Kolb 1996).

## 3. The effect of nova outbursts

The effect of nova outbursts on the distribution of mass transfer rates in CVs with similar orbital period (henceforth transfer rate spectrum for short) has been revisited by Kolb et al. 2001. Prior work by Shara et al. (1986) pointed out that as a result of the ejection of a nova shell the transfer rate can drop by one or more orders of magnitude. They hypothesized that CVs “hibernate”, i.e. most of the time between nova outbursts they are effectively detached with no mass transfer, and accrete only in short intervals that precede a classical nova outburst. However, Shara et al. (1986) did not consider how the transfer rate $`\dot{M}`$ varies with stellar radius $`R`$ and Roche lobe radius $`R_\mathrm{L}`$, and how both of these radii change with time.

Adopting a simple exponential law,

$$\dot{M}\mathrm{exp}\left(\frac{RR_\mathrm{L}}{H}\right),$$

where $`H`$ is the photospheric scaleheight (e.g. Ritter 1988), and assuming that the star’s thermal relaxation is the same as in the mean evolution (averaged over individual nova outbursts), reveals an area of parameter space where significant widths of the transfer rate spectra can occur. Generally, this occurs when the ejected envelope mass $`M_{\mathrm{ej}}`$ is larger than a critical limit, and the effect is strongest if the ejection is ballistic, i.e. if the affected shell does not gain additional AM while ejected.

Conditions for wide spectra are most favourable if the WD mass ($`M_1`$) is small and the system is dynamically unstable against conservative mass transfer. For a typical WD mass, $`M_10.60.7M_{}`$, the transfer rate spectrum is wider than 1 dex if

$$M_{\mathrm{ej}}>ϵ\times 3\times 10^4M_{}.$$
(1)
The scale factor $`ϵ`$, the ratio of photospheric scaleheight $`H`$ to stellar radius $`R`$ in units of $`10^4`$, is about unity for CV donor stars. Nova ejecta masses of this order have been suggested both by theoretical model calculations (e.g. Prialnik & Kovetz 1995) and by observational estimates (e.g. Starrfield 1999).

Systems that are dynamically unstable in the absence of nova outbursts do not develop the instability as the transfer rate is kept in check through nova outbursts. The phase space of systems that are dynamically unstable in the absence of novae but dynamically stable in the presence of novae is quite large, and may be populated by a non–negligible fraction of the observed CVs. This is illustrated in Figure 1 (right panel) where we plot the critical lower limit on the WD mass for dynamically stable mass transfer, both for the case of conservative mass transfer (upper curve) and for mass transfer with ballistic nova ejections (lower curve). The dashed lines (from top to bottom) indicate the lower mass limit on C/O WDs, the upper mass limit on He WDs and the lower mass limit on He WDs as obtained by standard population synthesis calculations (e.g. Willems, Mundin & Kolb, this volume).

Any CVs residing in between the two stability curves contribute to a significant widening of the mass transfer spectra if condition (1) is met. The overall impact of nova outbursts on the transfer rate spectra must be assessed with population synthesis calculations. As an illustraton, Figure 2 shows results from a standard CV model population, i.e. with standard assumptions about the common envelope evolution and magnetic braking strength, following the method used by Howell, Rappaport & Nelson 2001. The solid curves refer to a population with constant ejecta mass $`M_{\mathrm{ej}}=3\times 10^4M_{}`$, while the dashed curves represent a population where the transfer rate was averaged over individual nova outbursts.

The figure shows that the transfer rate spectra with novae are wider than in the averaged model for periods in the range $`3<P/hr<5`$, while the differential effect is negligible everywhere else. The mean value of the two different distributions is necessarily the same and set by the AM loss rate, while the width of the nova–affected distributions is insensitive to the adopted systemic AM loss rate. In other words, the effective widening of the spectra works also for a much less efficient magnetic braking rate.

A higher weighting of brighter (higher $`\dot{M}`$) systems would skew the distributions towards larger values of $`\dot{M}`$. In that case the observed average transfer rate of the weighted distribution immediately above the gap is larger than the secular mean rate set by the AM loss rate. This could alleviate the problem that the transfer rate at the upper edge of the gap required to explain the width of the period gap is smaller than the critical rate for disc instability, in apparent conflict with the observed underabundance of dwarf novae in the $`34`$ hr period interval (e.g. Shafter 1992).

At longer orbital periods the effect of novae is negligible. Podsiadlowski, Han & Rappaport (2001) point out that the relatively large fraction of CVs with a somewhat evolved donor star could give rise to a greater width of the mass transfer spectra above 5 hrs than the ones shown in Figure 2.

## 4. CVs on the edge

The shape of the CV period distribution below 2 hrs continues to challenge the standard model of CV evolution. The well known short–period cut–off at $`78`$ min is thought to represent the period minimum that systems with hydrogen–rich donor stars reach when their evolution is driven by gravitational radiation (e.g. Paczynski 1981, Kolb & Baraffe 1999, and references therein). The problems with this interpretation are summarised in Figure 3. The upper panel shows the lower end of the CV period distribution ($`P<2`$ hrs; data from an updated version of Ritter & Kolb 1998; see also Downes et al. 2001). Even with the most up–to–date input physics which is so successful in describing very low–mass stars and brown dwarfs (e.g. Chabrier & Baraffe 2000) the calculated value for the minimum period $`P_{\mathrm{min}}`$ is more than $`10`$ min smaller than the observed cut–off (middle panel), consistent with earlier findings (e.g. Kolb & Ritter 1992). Perhaps even more seriously, the very fact that systems evolve through a state with $`\dot{P}=0`$ implies that there will be an accumulation of systems at the period minimum (a “period spike”). The lower panel of Figure 3 shows the discovery probability $`1/\dot{P}`$ for an ensemble of CVs that form continuously with $`0.21M_{}`$ secondary mass and $`0.6M_{}`$ WD mass.

Various mechanisms have been proposed that would move the calculated period minimum $`P_{\mathrm{min}}`$ to longer periods. Most assume an additional form of AM losses, e.g. magnetic braking. Numerical experiments show that a rate of roughly three times the rate given by gravitational radiation (GR) alone is required to raise the calculated $`P_{\mathrm{min}}`$ to the observed value (Kolb & Baraffe 1999). The magnetic braking strength claimed by Andronov et al. 2001 is about two times the GR rate. Barker & Kolb (2002; see also this volume) show that consequential AM losses such as accretion disc winds or a magnetic propeller would raise the period minimum by $`<5`$ min. Increased AM loss is preferred by Patterson (2001), as it both avoids an overpopulation of the “CV graveyard” — systems with a brown dwarf donor star that have evolved past the period minimum and continue to transfer mass at a very low rate — and removes a claimed systematic discrepancy in the mass–radius relation of low–mass single stars and donors in short–period CVs.

Other attempts to reconcile the calculated and observed minimum period focussed on the Roche model and the approximation of the lobe–filling component in a semi–deatched binary by a one–dimensional stellar model. Nelson et al. (1985) used corrections for rotational and tidal deformations of the donor star (Chan & Chau 1979) and found a differential increase in $`P_{\mathrm{min}}`$ by 10%. However, more recently Kolb & Baraffe (1999) obtained a much smaller differential effect ($`1`$ min) using the same prescription but up–to–date stellar models.

Using a polytropic equation of state, Rezzolla et al. 2001 constructed sequences of self–consistent three–dimensional hydrostatic models of binaries to investigate the quality of the Roche approximation where the component stars are treated as point masses rather than extended mass distributions. The volume–equivalent radius of the critical lobe, the orbital AM, and the GR AM loss rate found from the self–consistent models agree with those from the Roche model to within 1-2%. They also confirmed the validity of the commonly used approximation for the Roche lobe radius in units of the orbital separation given by Eggleton (1983).

As an independent check, Renvoizé et al. (2001) calculated three–dimensional polytropic stars using SPH models with 20,000 particles. Significantly, the volume–equivalent radius of the star at the start of mass loss through the $`L_1`$ point (as depicted in Figure 4) is larger than the radius of the same star when it is non–rotating and in isolation. The radius increase depends on the polytropic index $`n`$ and is 6% for $`n=1.5`$ (as appropriate for fully convective stars near $`P_{\mathrm{min}}`$) and $`>10\%`$ for $`n=3`$, the case shown in the Figure.

A simple rescaling of the radius of one–dimensional models by a factor 1.06 increases the calculated minimum period of CVs that are driven by GR by about 5 mins ($`P_{\mathrm{min}}71`$ min if $`M_1=0.6M_{}`$). This is still well short of the cut–off period; a radius increase of almost $`20\%`$ is required to match up $`P_{\mathrm{min}}`$ with the cut–off (e.g. Barker & Kolb 2002). However, it is not straightforard to include the findings from the SPH star into one–dimensional stellar models. The expanded star is likely to have a lower luminosity and hence a longer thermal time. This in turn could lead to period bounce at a longer period. Such a second–order effect cannot be modelled by a simple rescaling.

None of the considerations presented so far address the problem of the missing period spike at $`P_{\mathrm{min}}`$. Barker & Kolb (2002) calculated models that effectively “smear out” the period spike, e.g. by allowing a mixed population with consequential AM losses of varying efficiency, or with donor stars that are subject to different degrees of bloating, e.g. due to the effect of magnetic pressure inside the donor star (as suggested by D’Antona 2000). They find that none of the theoretically calculated period distributions fits the observed one as well as a distribution that is simply flat in $`P`$. It is significant that the missing period spike affects both magnetic and non-magnetic CVs alike, so it seems unlikely that the cause should be related to effects in the accretion disc.

The radical hypothesis put forward by King & Schenker (2002) to abandon the identification of the cut–off with the minimum period altogether may appear to solve both problems, the mismatch and the missing spike, in one stroke. The suggestion that the cut–off is due to an age limit of the population implies that observed CVs close to 80 min orbital period are still evolving to smaller periods. The age limit hypothesis appears to reconcile current stellar models for short–period CVs that are driven by GR with the observed period distribution. The same models then tell us that systems that are currently at 78 min orbital period will continue to decrease $`P`$ for another 0.8 Gyr or so, before they reach the minimum period at about 70 min. This is a measure of how finely tuned the age limit has to be. More importantly, the mass of all hydrogen–rich donor stars in CVs that adhere to the age limit must be larger than $`0.09M_{}`$. Although there is no CV known where we have unambiguous proof that the donor mass is smaller than this limit, there are several excellent candidates where this is likely, most notably WZ Sge (e.g. Patterson 1998, but see Steeghs et al. 2001). The quest to identify brown–dwarf donor CVs should continue; the first few found would make the age limit hypothesis history.

## 5. Conclusions

This is an exciting era, challenging our understanding of the evolutionary state of CVs and X–ray binaries. Old paradigms will have to be re–assessed. The magnetic braking of rapidly rotating young, single stars seems to be different from the braking that old, rapidly rotating CV secondaries experience. The best stellar models to date continue to ignore the observed period minimum, and despite an ever increasing number of CVs with determined orbital period they still refuse to pile up at 80 minutes. Do we really understand the long–term evolution of CVs? The answer must be a clear no.

## References

Andronov N., Pinsonneault M., Sills A. 2001, ApJ, submitted (astro-ph/0104265)

Baraffe I., Kolb U. 2000, MNRAS, 318, 354

Barker J., Kolb U. 2002, in preparation

Chabrier G., Baraffe I. 2000, ARA&A, 38, 337

Chan K.L., Chau W.Y. 1979, ApJ, 233, 950

D’Antona F. 2000, in Evolution of Binary and Multiple Star Systems, eds. Ph. Podsiadlowski et al. ASP Conf Series 299, 21

Downes R.A., Webbink R.F., Shara M.M., Ritter H., Kolb U., Duerbeck H. 2001, PASP, 113, 764

Eggleton P.P. 1983, ApJ, 268, 368

Howell S.B., Nelson L.A., Rappaport S. 2001, ApJ, 550, 897

King A.R., Schenker K. 2002, this volume

King A.R., Kolb U. 1995, ApJ, 439, 330

Kolb U. 2001, in The influence of binaries on stellar population studies, ed. D. Vanbeveren, Kluwer, Dordrecht, ASSL, 264, 321

Kolb U. 1996, in Cataclysmic Variables and Related Objects, eds. A. Evans, J.H. Wood, Kluwer, Dordrecht, IAU Coll. 158, 433

Kolb U., Baraffe I. 1999, MNRAS, 309, 1034

Kolb U., Ritter H. 1992, A&A, 254, 213

Kolb U., Rappaport S., Schenker K., Howell S. 2001, ApJ, in press (astro-ph/0108322)

Krishnamurthi, Pinsonneault M.H., Barnes S., Sofia S. 1997, ApJ, 480, 303

Livio M. Pringle J. 1994, ApJ, 427, 956

Nelson L.A., Chau W.Y., Rosenblum A. 1985, ApJ, 299, 658

Paczyński B. 1981, Acta Astr., 31, 1

Patterson J. 1998, PASP, 110, 1132

Patterson J. 2001, PASP, 113, 736

Podsiadlowski Ph., Han Z., Rappaport S.A. 2001, MNRAS, submitted (astro-ph/0109171)

Podsiadlowski Ph., Saul Rappaport S.A., Pfahl E. 2001, MNRAS (astro-ph/0107261)

Prialnik D., Kovetz A. 1995, ApJ, 445, 789

Rappaport S.A., Verbunt F., Joss P.C. 1993, ApJ, 275, 713

Renvoizé V., Baraffe I., Kolb U., Ritter H. 2001, in preparation

Rezzolla L., Uryu K., Yoshida S. 2001, MNRAS, in press (gr-qc/0107019 )

Ritter H., Kolb U. 1998, A&AS 129, 83

Shafter A.W. 1992, ApJ, 394, 268

Shara M., Livio M., Moffat A.F.J., Orio M. 1986, ApJ, 311, 163

Sills A., Pinsonneault M.H., Terndrup D.M. 2001, ApJ, 540, 489

Skumanich A. 1972, ApJ, 171, 565

Spruit H.C., Ritter H. 1983, A&A, 124, 267

Spruit H.C., Taam R.E. 2001, ApJ, 548, 900

Starrfield S. 1999, Physics Reports, 311, 371

Stehle R., Ritter H., Kolb U. 1996, MNRAS, 279, 581

Steeghs D., Marsh T., Knigge C., Maxted P., Kuulkers E., Skidmore W. 2001, ApJL, submitted (astro-ph/0109311)

Townsley D.M., Bildsten L. 2002, this volume

Warner B. 1995, Cataclysmic Variables, CUP

### Acknowledgments.

I thank Isabelle Baraffe, John Barker, Steve Howell, Andy Norton and Saul Rappaport for help in compiling this review.