text:
# Dwarfs after Mergers? The case of NGC 520, NGC 772, Arp 141, NGC 3226/7, NGC 3656 and Arp 299

## 1 Introduction

Galaxy interactions provide strong clues about the structure and dynamics of galaxies, and a powerful mechanism for the origin of very massive central star formation events. Besides, the development of long tidal tails of stellar and gaseous material during the encounter, which can be observed and computed in numerical simulations, has posed the question of the fate of this tidally ejected material. Zwicky (1956) was the first to propose that some of this material might survive in the long term in the form of self-gravitating dwarf galaxies. Evidence in favor of this scenario for the origin of dwarfs comes from the identification of blue condensations in the tips of the tails of galaxy mergers with typical Hi masses ranging from $`5\times 10^8`$ to $`6\times 10^9`$ $`M_{}`$. Mirabel et al. (1992) and Duc et al. (1998) have shown that objects similar to dwarf irregulars or blue compacts are found in the tails. These small galaxies of tidal origin are potential progenitors of detached systems, namely, isolated dwarf galaxies (Sanders & Mirabel 1996, Braine et al. 2001).

Numerical simulations of interactions between disk galaxies show that tidal tails are the natural result of such processes (Barnes & Hernquist 1996). Also, if the interaction is strong enough, self-gravitating clumps of gas and stars are likely to form along the tails (Schweizer 1978). Several examples of stellar condensations at the tips of long tidal tails have been reported, like NGC 7252 (Hibbard et al. 1994), Arp 105 (Duc & Mirabel 1994) and Arp 245 (Duc et al. 2000). The possibility that these clumps become independent dwarf companions has been proposed by Barnes & Hernquist (1992) and Elmegreen et al. (1993). In addition, some material – both gaseous and stellar – is ejected from the disks of the parent galaxies, forming lower mass complexes that are probably not self-gravitating. It is expected that only the outermost ejected material may eventually abandon the potential well of the parent galaxy and form dwarf galaxies. The innermost material would slowly fall back to the parent galaxy. Given that the further out the material is, the longer is the time for return, dwarf objects in the surroundings of the interacting pairs would be expected for a long time, even well after the relaxation time. These self-gravitating condensations of stars and Hi, the so-called Tidal Dwarf Galaxies (henceforth TDG), would fade significantly in optical wavelengths in the Gyr that follows their formation, as soon as the gas is appreciably exhausted by the starburst (Weilbacher et al. 2000). In addition, tidal stripping from the parent galaxy may be in the long term a significant source of mass loss, thereby reducing the TDG to a very low surface brightness dwarf (Mayer et al. 2001). All these circumstances conspire to make genuine dwarf galaxies hard to detect, unless Hi cartography as well as spectroscopic surveys around mergers are performed.

In the present study we investigate if strong interactions between galaxies can produce a significant enhancement on the number of dwarf galaxies around the system. Lacking spectroscopic data, which are very difficult to obtain due to the faintness and small size of the large majority of TDG candidates, we choose to address this question statistically, by comparing the photometric sample of extended objects \- or dwarf candidates- catalogued in the previous work of Deeg et al. (1998, hereafter Paper I) with background galaxy number counts. By doing so, we expect to include both young and evolved TDG, arriving at an estimate of the abundance of TDG formed by field interacting galaxies. From the fields around 15 interacting systems catalogued in Paper I, we have selected those fields that were observed in two filters ($`V`$ and $`R`$), had a statistically significant number of detections of extended objects (EO), and in which the central galaxies have similar redshifts (five of the sample galaxies have vr between 2266 and 3132 km s-1, NGC 3226/7 is at 1151 km s-1). Thus, a subsample of fields around six galaxies – namely NGC 520, NGC 772, NGC 3226/7, NGC 3656, Arp 141 and Arp 299 - is analyzed in this paper. Among these, Arp 299 was observed in one filter only, but the high number of detections of EO in this field lead to its inclusion in this study.

In the following, we describe the observations and the detection procedure in Sect. 2. In Sect. 3, we outline the properties of known TDG within the framework of a dynamic evolution scenario and describe the statistical analysis that was performed in order to quantify the possible number of TDG. Conclusions are given in section Sect. 4.

## 2 Observations and Detection of Candidates

Images of the fields were taken in December 1992 with the 2.5m Isaac Newton Telescope at the Observatorio del Roque de los Muchachos, La Palma, with the EEV5 CCD-camera in prime focus, giving an effective field of about 10.5$`\mathrm{}\times 10.5\mathrm{}`$ and 0.55 $`\mathrm{}`$ pix-1 resolution. The observations took place in bright-grey time which affected the limiting magnitude. Atmospheric transparency conditions were very good and seeing ranged from 0.9 to 1.5 $`\mathrm{}`$. Images were taken with the $`V`$ and $`R`$ Cousins filters, with exposure times as given in Table 1. In order to determine the stellar PSF, frames with short 20 s exposures were also taken, in which known bright field stars are not saturated. Landolt (1992) standards were observed for flux calibration.

Table 1 lists relevant parameters of the target list and of the observations.

Raw images were zero-subtracted and flatfielded in the usual manner. The bias frame used was taken from averaging all bias frames in the observations. Also, global flatfields were obtained by averaging the individual sky flatfields taken at the entire observing run. This was possible, since the stability of the flatfields between nights was better than 1%. An absolute flux calibration was performed based on the observations of Landolt standards at different airmasses. The standards were measured using fixed-aperture photometry. For more details about observations, calibration procedure and classification methodology, see Paper I. A procedure was developed to identify objects over the background of the CCD frames, discriminating small galaxy-like objects and star-like objects (see the flow diagram in Paper I). Identified objects were then parameterized by the fitting of elliptical Gaussians, as is described later in the text. The total error in the photometric calibration is about 0.15 mag.

Objects with a size close to the PSF were considered stellar (SO) whereas those with an extent of more than 1.3 times the PSF in either $`x`$ or $`y`$ direction were classified as extended (EO). Objects whose FWHM ranges between 1 and 1.3 times the PSF were classified as ambiguous. All saturated objects, also those having widths larger than the PSF, can safely be classified as stellar, too. Table 2 gives the number of EOs, SOs and ambiguous identifications found around the galaxies of the sample. The selection method used certainly does not exclude that the lists of galaxies still include some semi-resolved double or multiple stars. Also the lists of stars may contain a substantial number of galaxies with a near–stellar PSF. In galaxy–star separations done in a way comparable to ours, Metcalfe et al. (1991) found that later spectroscopy revealed 29% of the star–like objects to be compact galaxies, and 5% of the galaxy identifications turned out to be stars.

It might be argued that some of the detected objects are globular clusters formed during the interaction. In fact, Zepf & Ashman (1993) and Ashman & Zepf (1992) reported two-aged populations of globular clusters around merger remnant candidates. The colors of our objects would be compatible with those of typical globular clusters (Harris 1996), thus opening the possibility that they are globular clusters that had originated during the interaction. However, globular clusters should be detected as point-like sources at the distances of the interacting galaxies of our sample. In addition, only the faintest EOs in our sample would have luminosities comparable to those of globular clusters.

Figure 1 shows calibrated $`R`$ band images of the sample fields, in which the EO, SO and ambiguous objects that have been detected are marked. Below, we give some comments on the galaxy fields:

* NGC 520: This is a gas rich system that is widely studied in the literature. An interesting multifrequency analysis of this interacting system can be found in Hibbard & van Gorkom (1996). Two dwarf objects with radial velocities similar to that of the parent galaxy, detected in both the optical and in Hi, are present in the neighborhood of the merger, and one of them (UGC957) appears permeated by the northern tidal tail of the system. This latter object might be the result of the interaction.
* NGC 772: Clear disruption signatures can be found in the spiral arms of this galaxy, which are likely due to interactions with some of the satellites detected around it (see Zaritsky et al. 1997).
* Arp 141 (UGC 3730): This is quite an irregular system with an optical extension towards the South. Burbidge & Burbidge (1959) attributed the peculiar nature of this interacting pair to the disruption of a spiral galaxy by an elliptical one in a close encounter.
* NGC 3226/7: This peculiar system is composed by NGC 3226 and NGC 3227. The galaxies are embedded in an extended diffuse envelope. Hi was detected along the two plumes emerging from NGC 3227 (see Mundell et al. 1995).
* NGC 3656: This shell galaxy was studied in detail by Balcells (1997). Two low surface brightness tidal tails were reported. At the tip of the western tidal tail, two dwarf galaxies were detected in broad band images. However, it is not certain that these objects are associated with NGC 3656.
* Arp 299: This system consists of two interacting spiral galaxies, IC 649 and NGC 3690. Multi-wavelength HST observations of NGC 3690 highlight two massive star–forming condensations $``$ 3 kpc from the merger nucleus that might become detached dwarf galaxies (Alonso–Herrero et al. 2000). An Hi–rich tidal tail towards the North is apparent in larger optical frames. Surprisingly, no dwarf galaxies were detected along this tail (Hibbard & Yun, 1999).

## 3 Results

### 3.1 TDG properties and the collision process

TDGs, objects which are self-gravitating entities, formed out of the debris of a galactic gravitational interaction (Duc et al. 2000) are characterized by a luminosity comparable to that of typical dwarf galaxies, an extended morphology, high metallicity, and blue colors ($`BV0.5`$ and $`VR0.4`$) as a result of an active starburst. These properties, however, have been outlined for objects found at the tip of the optical tails, some of them kinematically decoupled from the surrounding material but still visually embedded in the collision debris. On the other hand, more numerous ( $`4`$ TDG per collision) TDG candidates have been found, both in broad band optical images (Weilbacher et al. 2000) and in H$`\alpha `$ emission (Iglesias-Páramo & Vílchez 2001), distributed all along the tidal features. However, it remains unclear whether these structures might become truly detached dwarfs in the future. How many TDG does an average merger form? Computer simulations of encounters between equal-mass disk galaxies (Barnes, private communication) show the formation of $`100`$ bound structures along the tidal tails, with a distribution of distances from the merger that is very broad and evolves with time as the tail expands. Only a few bound systems have masses above $`5\times 10^8`$ $`M_{}`$ and they are located randomly throughout the tail, at an approximate mean distance of 100 kpc. Among these most massive self-gravitating structures, those with larger radii are the best candidates for becoming independent dwarf galaxies in the long term (Barnes & Hernquist 1992, Elmegreen et al. 1993). Most of the other bound structures would fall back to the central galaxy within a Gyr. One should however consider that in the case of interactions between unequal mass galaxies, the length of the tidal tail should be much smaller than in the case reported above. Correspondingly, dwarf candidates would be located much closer to the nuclei of the interacting pair.

One might imagine a dynamic scenario in which, depending on the nature (morphology, mass) of the galaxies involved in the collision and the stage at which the interaction is being observed, TDG do not necessarily have the properties observed in their young counterparts. TDGs could have formed, developed a massive starburst, and faded away, as well as decoupled from the main stream of merger debris – Weilbacher et al. 2002 reports on 7 knots along the tails of AM 1353 272 which could be moving away from the parent tidal stream with speeds as large as 100 km s-1. In such a picture, it would be possible to find TDGs without blue excess located well away from the tidal features, in numbers that would be hard to predict. Weilbacher et al. (2000) have computed the future luminosity evolution of the tidal knots detected in a sample of 10 interacting galaxies, assuming that a single starburst episode occurs. The colors become significantly bluer (at first mainly in $`VR`$) once the burst commences, and redden again in the next Gyr, as the stars formed in the burst age. Therefore, TDGs around interacting galaxies may comprise a large population, perhaps hidden among background faint galaxies in the neighborhood of the interacting pair. A series of tests to assess this scenario was undertaken and is described in the next section.

### 3.2 Analysis

We have computed the number density distribution $`\mu `$ of EOs and SOs in each field (Fig. 2). This quantity has been used to look for a possible dependence of $`\mu _{\mathrm{EO}}`$ and $`\mu _{\mathrm{SO}}`$ with distance from the parent galaxies, whose center was taken as reference point (RP). Each CCD frame was divided in concentric rings around the RP. Only rings close to the center are completely present on the frames. To include rings further out, an area-normalization algorithm was employed that takes the irregular shapes of incomplete rings into account. Care was taken, through a revision of the flatfield frames, that vignetting is negligible, as it would cause lower detection limits far from the optical axis.

As can readily be seen, there is no clear evidence for a characteristic distance where EOs are preferentially appearing. Within the error bars, the density distributions are essentially flat. Note that the last bin includes the outermost EO detected near the edges of the CCD frames and their small numbers produce unreliable statistical averages. This result does not depend on color. In addition, as is expected for a random spatial distribution, EO and SO distributions do not show any statistically significant difference.

Figure 3 also shows the distribution (small points) of the radial distance of the detected EO as a function of the $`(VR)`$ color index. The thick symbols stand for the average values over a given radial distance interval. Error bars correspond to the standard deviation. Overall, the averaged values are – within the error bars – independent of the distances to the central objects. Therefore, we conclude that the color distribution of EOs is fairly homogeneous, with no privileged position or distance in any of the studied fields.

Figure 4 shows the $`(VR)`$ color index as a function of the $`V`$ magnitude for the EOs detected in all the frames. Crosses with error bars correspond to the averaged values per half unit magnitude and one standard deviation respectively. For comparison, we have over-plotted the mean values (filled circles) and $`1\sigma `$ standard deviation of the colors (dashed lines) for the photometric sample of background galaxies reported by Arnouts et al. (1997). As can be seen, both sets of averaged values agree very well with each other. There is therefore no indication for a foreground population with a distinct color.

Finally, the number of EOs in each field, in intervals of 0.5 $`R`$ magnitudes as well as over the entire magnitude range, was computed and compared with the number of expected background galaxies (Table 3). The purpose of the binning is twofold: first, it is required to compute the errors in the number counts consistently and, second, it allows us to search for possible excesses over background counts within specific magnitude ranges. The abundances of expected background galaxies were calculated from galaxy counts by Metcalfe et al (1991). Possibly, galaxy counts from observations of nearby blank fields might have added support to these results. However, due to the similarity between our star-separation procedure and the one used by Metcalfe et al. (1991 – both are based on a source size criterion) we are confident that this would not have changed the results significantly.

Three factors contribute to uncertainties in the galaxy number counts: Poisson noise, the spatial distribution or clustering properties of galaxies, and the magnitude error. Following Roche et al. (1993), for a square field of side-length $`\theta `$ (degrees), with N galaxies observed in it and a two point angular correlation function $`\omega (\theta )`$ amplitude A, the field-to-field variance in N is computed using:

$$\sigma ^2=N+2.24N^2A\theta ^{0.8}$$
(1)
where a $`\theta ^{0.8}`$ power-law was assumed for $`\omega (\theta )`$ (recent results agree very closely with this value; Connolly et al. 2002), and A has been obtained from Figs. 4 (B band) and 5 (R band) of Roche et al. (1993) (a conversion log(A(R))=log(A(B))$``$0.5 was applied). These errors were calculated for both the detected (EOs) and the expected (Metcalfe’s) counts: $`\sigma `$ and $`\sigma _e`$ respectively, in Table 3. Systematic magnitude errors may also cause changes in the galaxy number counts. Such zero-point errors relate to galaxy count errors through the slope of the log\[N(m)\] distribution (see Huang et al. 1997 for details). In paper I, the uncertainty in the zero point is estimated to be $`\pm `$ 0.1 mag, a value that causes number count errors that are much smaller than the field-to-field variations. Metcalfe et al. also reports on a negligible contribution of the systematic magnitude error to the overall uncertainties.

Table 3 shows the probable number (dN) and maximum number (Nmax) of TDG candidates obtained by following the procedure described previously. We note that in some fields – NGC 520 and Arp 141 – there are a few more objects than expected from background counts, in bins 1-2 magnitudes brighter than the completeness limit (i.e. in the range R=18.5-19.0). NGC 772 shows a $`2\sigma `$ detection in the same range, but one $`2\sigma `$ deviation in a sample of 36 points would be expected just from Poisson statistics and therefore it is not statistically significant. The fact that no discrepancies are found when the total excess of observed counts over background counts is computed (over the entire magnitude range: last row in each field of Table 3.), indicates that the TDG population must be very small in comparison to background objects of comparable size and brightness. In other words, our results rule out that a typical galaxy merger forms 10 or more TDG. Individual cases may produce that quantity, but not more than a few – if any– TDG can be expected from a typical galaxy merger.

## 4 Conclusions

Our results do not support evidence for a significant population of dwarf galaxies around strongly interacting galaxies: the total probable number of TDG is a small negative number in all galaxy fields except NGC 772. In this latter field, only one magnitude bin shows a $`2\sigma `$ detection of TDG, which is statistically insignificant. Tidal condensations may turn into bound systems abundantly, but as self-gravitating, kinematically detached entities, we expect only a few TDG per collision to be formed. The value indicated by our results is $``$ 1 TDG per merger, although as many as $``$ 10 TDG cannot be ruled out in individual cases. For three of the fields – namely NGC 520, NGC 3226/7 and NGC 3656 – TDG candidates have previously been reported (see Balcells 1997; Hibbard & van Gorkom 1996; Mundell et al. 1995), but no more than two per interacting pair and always associated with tidal tails. This result is also consistent with previous Hi and optical observations of six interacting systems (van der Hulst 1979; Hibbard & van Gorkom 1993; Hibbard & van Gorkom 1996; Hibbard et al. 1994), where five dwarf galaxy candidates were reported.

Our results agree with simple theoretical models (Elmegreen et al. 1993), which predict that only the outermost material will gain enough angular momentum and energy to become an independent dwarf galaxy. Most of the bound features along tidal tails found in computer simulations (e.g. Barnes & Hernquist 1992) and a substantial fraction of the luminous condensations reported in the literature (e.g. Weilbacher et al. 2000) are therefore expected to be short-lived structures: truly self-gravitating dwarfs that are the result of galaxy mergers would not be frequent.

###### Acknowledgements.

The Isaac Newton Telescope is operated on the island of La Palma by the Royal Greenwich Observatory in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautic and Space Administration. We thank our referee, J.E. Hibbard, whose valuable comments have certainly contributed to improve and clarify this paper. This work has been partially supported by the Spanish DGC (Grant No. AYA2001-3939). EDD is grateful to the European Union Research Training Network Program for its support.


format: text

text:
# Chromatin dynamics: Nucleosomes go mobile through twist defects

## Abstract

We study the spontaneous ”sliding” of histone spools (nucleosomes) along DNA as a result of thermally activated single base pair twist defects. To this end we map the system onto a suitably extended Frenkel-Kontorova model. Combining results from several recent experiments we are able to estimate the nucleosome mobility without adjustable parameters. Our model shows also how the local mobility is intimately linked to the underlying base pair sequence.

The genetic information of all higher organisms is organized in huge beads-on-a-chain arrays consisting of centimeters to meters of DNA wrapped around globular aggregates of so-called histone proteins. The basic unit of chromatin, the nucleosome, is a tiny $`10\times 5\times 6`$ nm sized spool composed of $`147`$ base pairs (bps) DNA tightly wrapped around an octamer made from 8 histone monomers. Each nucleosome is connected via a stretch of ”linker” DNA to the next such protein spool. The wrapped DNA, being coiled in $`1\frac{3}{4}`$ turns of a left handed helix with radius $`4.2`$ nm, is strongly distorted from its preferred straight ground state due to strong interactions with the histone octamer, namely short range electrostatics (between the negatively charged DNA sugar-phosphate backbone and the positively charged octamer surface) and through extensive hydrogen-bonding – both localized at 14 discrete interaction patches helically arranged along the octamer surface luger97 .

Higher order structures, from the 30nm-chromatin fiber up to the highest level of DNA condensation, the fully folded chromosome, are designed to achieve a huge DNA volume fraction. They all rely on the significant stability of the nucleosome complex. On the other hand, fundamental life processes like transcription (making RNA offprints from the underlying DNA) and DNA replication seem to be in conflict with the picture of a stable nucleosome, as they are all performed by protein machines that track the DNA helix. The latter inevitably implies that every DNA bound obstacle (protein) has to be penetrated or even completely removed from its DNA target. In fact, the numbers are quite dramatic: A typical gene extends over hundreds of nucleosomes, each contributing $`3040`$ $`k_BT`$ net adsorption energy polach95 ; brower-toland02 . Also other mechanisms like the activation of genes rely on regulatory protein binding to specific DNA sequences that are often covered by nucleosomes making them inaccessible.

A key to the understanding of these seemingly contradictory features might be the physical phenomenon of thermally driven nucleosome ”sliding” along DNA (also called nucleosome repositioning) which has repeatedly observed in well-defined in vitro experiments beard78spadafora79 ; pennings91 ; flaus98gottesfeld02 , reviewed in Ref. schiessel03b . Spontaneous repositioning is strongly temperature dependent; at room temperature nucleosomes move a few tens of bps within an hour pennings91 . Despite clear evidence for repositioning the underlying mechanism has been the matter of longstanding controversy, especially due to the lack of any quantitative theoretical treatment of nucleosome statics and dynamics that has to rely on the detailed knowledge of the molecular structure and its underlying parameters.

Only very recently – since the documentation of the high resolution X-ray structure luger97  and the presentation of other new experiments polach95 ; brower-toland02 ; flaus98gottesfeld02  – this has become possible. First theoretical models of nucleosome repositioning schiessel01c ; kulic02  assume that it is based on the formation of DNA loop defects that form on either end of the nucleosomal DNA followed by their thermal diffusion around the octamer, similar to the de Gennes-Edwards reptation mechanism. This model seems to be successful in explaining the apparent 10 bps quantization of the nucleosome ”jump” length pennings91  and it also reproduces the observed diffusion constants. Recent experiments flaus98gottesfeld02 , however, indicate a more local 1 bp-step mechanism that cannot be understood within this model. This lead us here to consider an alternative mechanism: twist diffusion. The carrier of motion in this case is a twist defect that contains one missing or one extra bp.

The possibility of twist defects was demonstrated as soon as the high resolution crystal structure of the core particle (the octamer plus wrapped DNA) was resolved luger97 . In that study the core particles were reconstituted from palindromic 146 bp DNA and histones assuming that this would result in a complex with perfect two-fold symmetry. However, it turned out that one bp is localized at the dyad axis, the rest being divided into a 73 bp half and a 72 bp half. The missing bp of the shorter half is, however, not localized at its terminus but instead at a 10 bp stretch close to the dyad axis (cf. Fig. 4d in Ref. luger97 ). This is due to the attraction between the DNA termini of adjacent particles in the crystal that come close to mimic a bp step at the cost of forming a twist defect far inside the wrapped chain portion. This allows us to estimate the energy for a single defect to be smaller than the stacking energy that is $`1020k_BT`$ estimate .

In order to model the twist diffusion mechanism we map the nucleosomal DNA on a Frenkel-Kontorova (FK) chain of particles connected by harmonic springs in a spatially periodic potential (cf. Fig. 1). The original FK model was introduced more than sixty years ago to describe the motion of dislocations in crystals frenkel38 . In the meantime variants of this model were applied to many different problems including charge density waves floria96 , sliding friction braun97 ; strunz98 , ionic conductors pietronero81 ; aubry83 , chains of coupled Josephson junctions watanabe96  and adsorbed atomic monolayers uhler88 ; schiessel03 . Here, in the context of DNA adsorbed on the octamer, the beads represent the base pairs. The springs in between have an equilibrium distance $`b=0.34`$ nm and a constant that reflects the coupled DNA twist-stretch elasticity. Specifically

$$E_{elastic}\left(\left\{x_n\right\}\right)=\underset{k}{}C\left(\frac{x_{k+1}x_k}{b}1\right)^2$$
(1)
Here the conformation of the wrapped DNA is given by the set $`\left\{x_n\right\}`$ where $`x_n`$ is the position of the $`n`$th bp measured along the helical backbone; $`C70100k_BT`$ is the combined twist and stretch spring constant including the (here unfavorable) twist-stretch coupling kamien97marko97  and the summation goes over all bp associated with the wrapped DNA. In addition there is the external potential of the 14 contact points to the octamer with neighboring points being 10 bp apart luger97  that we model as follows

$`E_{ads}\left(\left\{x_n\right\}\right)`$ $`=`$ $`U_0{\displaystyle \underset{k}{}}{\displaystyle \underset{l=1}{\overset{14}{}}}\left(\left({\displaystyle \frac{x_k10bl}{a}}\right)^21\right)^2`$ (2)
$`\times \theta \left(a\left|x_k10bl\right|\right)`$
with $`\theta `$ being the Heaviside step function. The two parameters of the external potential, its depth $`U_0`$ and its width $`a`$, can be estimated as follows. $`U_0`$ represents the pure adsorption energy per point contact which follows from competitive protein binding polach95  to be of order $`6k_BT`$. The other parameter, $`a`$, can be estimated from the fluctuations of the DNA in the crystal measured by the B-factor (cf. Fig. 1b in luger97 ) at different nucleosome positions. The ratio of DNA helix fluctuations $`R_{fluct}=x_{middle}^2/x_{bond}^23`$ at positions between the binding sites and at the bound sites is a measure of DNA localization. Using a quadratic expansion of Eq. 2 one finds from a straightforward normal mode analysis that $`a=\left(5U_0/[\left(R_{fluct}1\right)C]\right)^{1/2}bb/2`$, i.e., the adsorption regions lead to a strong localization of the DNA. Knowing all involved parameters the total energy of the DNA chain confined in the nucleosome can be written down

$$E_{tot}=E_{elastic}+E_{ads}+E_{sd}$$
(3)
The last term $`E_{sd}`$ is the sequence dependent part of the total energy which we will neglect first. In the following we study the mechanism for thermal motion of DNA governed by $`E_{tot}`$. Generally two scenarios are possible: (i) The generation of kink-antikink pairs inside the nucleosome and (ii) a kink (or antikink) injection at either nucleosome end. Since the first mechanism is energetically roughly twice as costly than the second one, we will focus here on the (anti)kink injection mechanism only.

How and how fast does the kink step around the nucleosome? Due to the strong DNA localization at the binding sites ($`a/b<1`$) for a realistic range of parameters $`U_0`$ and $`C`$ the kink is localized either between two adsorption positions, i.e., smeared out over $`10`$ bp (denoted by the $`K_{10}`$ state), or between three of them, i.e., smeared out over $`20`$ bp (the $`K_{20}`$ state). It is obvious that the motion of a (anti)kink will consist of an alternation between $`K_{10}`$ and $`K_{20}`$ states similarly to an earthworm creep motion. To model this process we introduce the effective kink coordinate $`x_K`$ describing the coordinate of the DNA bp being pinned/depinned during a single kink step, so that $`x_K0`$ and $`x_Kb/2`$ correspond to $`K_{10}`$ and $`K_{20}`$, respectively, whereas $`x_Kb`$ means that the kink moved by one bp step. The Peierls-Nabarro potential experienced by the kink is then given by $`U_{PN}\left(x_K\right)=C_{\mathrm{𝑒𝑓𝑓}}\left(x_K/b1/2\right)^2U_0\left(x_K^2/a^21\right)^2`$ for $`0<x_K<b/2`$ and $`U_{PN}\left(x_K\right)=U_{PN}\left(bx_K\right)`$ for $`b/2x_K<b`$. Here $`C_{\mathrm{𝑒𝑓𝑓}}=\frac{2}{10\pm 1}C`$ with the ”$``$” sign referring to a kink (1 bp missing) and the ”$`+`$”sign to an antikink (1 additional bp). Depending on the ratio of parameters $`U_0`$ and $`C`$, the state $`K_{20}`$ corresponds to a local minimum or maximum of $`U_{PN}`$ whereas $`K_{10}`$ is always stable for the relevant parameter range. The rate for the kink step process is then given by the expression $`f_{step}=k_BTj_0/b^2\zeta _{\mathrm{𝑒𝑓𝑓}}`$ with $`j_0^1=\left(_0^1e^{U_{PN}\left(sb\right)/k_BT}𝑑s\right)\left(_0^1e^{+U_{PN}\left(sb\right)/k_BT}𝑑s\right)`$ and $`\zeta _{eff}=\frac{4\pi ^2}{10b}\mu _{spin}`$, the effective kink friction constant. Here $`\mu _{spin}=1.3\times 10^{20}Ns`$ is roughly the rotational friction for a single basestep levinthal56nelson99 . To determine the rate at which twist defects are formed at the entry/exit points of the DNA one can now use an argument similar to the one presented in Ref. schiessel01c : The ratio of the life time $`t_{life}`$ of a kink to the time interval $`t_{inj}`$ between two kink injection events at the end of the wrapped DNA portion equals the probability to find a defect on the nucleosome, i.e. $`t_{life}/t_{inj}N_{site}e^{U_{Kink}/k_BT}`$. Here $`N_{site}=13`$ denotes the number of possible positions of the defect between the 14 binding sites and $`U_{Kink}C/10`$ is the energetic cost for a single kink (cf. above).

How is the average life time $`t_{life}`$ of a defect related to $`t_{step}`$, the typical time needed for one step? This can be determined from the mean first passage times $`\tau _{left}`$ and $`\tau _{right}`$ for a defect that forms, say, at the left end to leave the nucleosome at the same or at the other end, respectively. From Ref. vankampen  one finds $`\tau _{left}=\left(25/6\right)t_{step}`$ and $`\tau _{right}=28t_{step}`$. Furthermore, the probability to leave at the left end is $`p_{left}=12/13`$ and at the right end $`p_{right}=1/13`$ vankampen  which gives the life time as the weighted average $`t_{life}=6t_{step}`$. Only a fraction $`p_{right}`$ of the defects reaches the other end and will lead to a repositioning step, i.e., the time of a 1bp diffusion step of the nucleosome along the DNA is given by $`T=t_{inj}/p_{right}`$. Putting all this together we arrive at $`T6b^2\zeta _{\mathrm{𝑒𝑓𝑓}}j_0^1/k_BT\mathrm{exp}\left(C/10k_BT\right)`$. For realistic parameter values $`C=100k_BT`$, $`U_0=6k_BT`$ and $`R_{fluct}=3`$ we find $`T10^3`$ s implying a nucleosome diffusion constant $`D=580\text{ bp}^2/\text{s}=`$ $`6.6\times 10^{17}\text{ m}^2/\text{s}`$.

Hence we find repositioning rates that are orders of magnitude faster than the ones observed in experiments pennings91 . Even worse, the experimental observation of an apparent 10 bp jump length pennings91  seems to be inconsistent with our predictions. We show now how these facts can be explained by the existence of additional barriers with a 10 bp periodicity. To do so we have to extend our simple model to deal with the quenched disorder stored in the DNA bp sequence. The sequence dependent anisotropic bendability, i.e., the propensity of DNA to bend in different directions with different elastic constants turns out to be essential. It has been known for long satchwell ; shrader89shrader90  that (A/T) rich dinucleotide steps (dns) prefer to face the octamer in the minor groove (i.e., at the octamer contact points) whereas (G/C) rich dns prefer to face the octamer in the major groove (i.e., between contact points). This reflects different propensities of the dinucleotides to widen or compress towards the DNA minor groove. To incorporate these anisotropic effects into our model we first note that the bending state of the DNA molecule is fully constrained by its helical path on the octamer surface. Moving a DNA sequence via twist diffusion by a few bp ($`<10`$ bp) along that path changes the relative rotational setting of the bent DNA with respect to its preferred bending direction causing an energetic penalty, whereas a motion by 10 bp restores the initial rotational setting. We address this by introducing a 10bp periodic ”bending field” $`F_{bend}\left(x\right)=\mathrm{cos}\left[2\pi x/\left(10b\right)\right]`$ attached to the octamer surface. We assume the DNA sequence to couple linearly to that field through ”bending charges” $`q_k`$ attached to each of the dns. This gives us finally the third term in Eq. 3:

$$E_{sd}=_kq_kF_{bend}\left(x_k\right)+m_k$$
(4)
In addition to the anisotropic term we also introduced here the isotropic bending parameters $`m_k`$ to include isotropic flexibility effects (which become important when the $`q_k`$’s vanish or average out). The summation involved is again over all base pairs footnote1  incorporated in the nucleosome. $`q_k`$ and $`m_k`$ both have units of energy and can be extracted from competitive protein binding experiments shrader89shrader90  for each of the 10 dns (AA, AT, GC…). To obtain a rough estimate we distribute the dns into three classes: 1) (G/C) containing dns, 2) (A/T) containing dns and 3) mixed dns (like AG, CT etc.) and treat the dns in each class as identical. Using the available experimental data shrader89shrader90 ; Cacchione97  we then arrive at $`q_{G/C}`$ $`95`$, $`q_{A/T}85`$ , $`q_{mixed}0`$ and $`m_{G/C}20,`$ $`m_{A/T}3`$, $`m_{mixed}7`$, where all energies are in cal/mol per dns.

It turns out that the nucleosome mobility depends strongly on the underlying bp sequence. When shifting the position of all beads by $`l`$ bp steps, $`x_k`$ $`x_k+lb`$, we find $`E_{sd}\left(l\right)=\left(A/2\right)\mathrm{cos}\left(2\pi l/10\varphi \right)`$ to vary as a cosine function of $`l`$ with phase $`\varphi `$ and amplitude $`A`$ determined by the DNA sequence, which is assumed to be appropriately periodic here. Arranging G/C and A/T tracts properly and taking the sequence dependent $`q`$ and $`m`$ values given above we can easily reach amplitudes $`A`$ (i.e. barriers to repositioning) that exceed $`1012`$ kcal/mol! A very effective sequence arrangement called the ”TG”-sequence which leads to a strong nucleosome stability and localization was experimentally constructed in Ref. shrader89shrader90  by putting G/C tracts around positions $`k=0,10,20\mathrm{}`$ and A/T tracts around $`k=5,15,25\mathrm{}`$. In our picture this means to put the ”bending charges” $`q`$ along the DNA such that they couple favorably to the bending field $`F_{bend}`$ for a distinct rotational setting whereas a 5 bp shift is extremely costly (cf. Fig. 2). The 5S-RNA sequence which was used in most nucleosome mobility experiments shows also the effect of an optimal rotational setting. It is less pronounced than in the ”TG” case, yet it is still detectable. More involved theoretical computations relying on molecular sequence dependent deformability parameters Mattei02  reveal barriers $`A56`$ kcal/mol for this particular sequence. The sequence dependent barrier height $`A`$ exponentially suppresses the bare (sequence independent) diffusion constant $`D`$ obtained above leading to the sequence dependent diffusion constant $`D_{sd}`$:

$$D_{sd}=DI_0^2\left(A/2k_BT\right)\frac{\pi j_0A}{12\zeta _{\mathrm{𝑒𝑓𝑓}}}e^{\left(A+C/10\right)/k_BT}$$
(5)
with $`I_0`$ being the modified Bessel function.

Equation 5 predicts that mobility experiments with highly anisotropic sequences like ”TG” (instead of the standard ”5S-RNA”) would find hardly any appreciable repositioning on the one hour timescale if it would be solely mediated via twist defects ($`D_{sd}=10^610^7\times D=10^410^5\text{ bp}^2/\text{s}`$). The typical path for a nucleosome to escape from such a rotational trap goes very likely via the previously considered loop formation mechanism schiessel01c ; kulic02  that allows ”tunneling” over sequence barriers, thus dominating over twist-diffusion for extremely anisotropic sequences. An experimental test for this prediction would be to increase the free DNA segment length which in this regime should strongly enhance the loop mediated mobility kulic02  whereas it would leave the twist diffusion unaffected. Going to the other extreme, in the most relevant case of random isotropically bendable sequences which make up more than 95% of the eucaryotic genome one should observe that the twist diffusion mechanism is strongly enhanced by 2-3 orders of magnitude as compared to the in vitro measurements on ”5S-RNA”.

In conclusion the following picture is implied: On physiological timescales the majority of genomic nucleosomes seems to be intrinsically highly mobile. However, only a small fraction ($`<5\%`$) of all nucleosomes has strongly reduced mobility due to anisotropic DNA sequences which they populate. We speculate that only the latter require the action of active (ATP consuming) remodelling mechanisms lorch99  making them hotspots and switching elements for global chromatin rearrangements.

We thank R. Bruinsma, K. Kremer, K. Luger, F. Müller-Plathe and J. Widom for helpful discussions.


text:
# Comment on: “Theory of the evaporation/condensation transition of equilibrium droplets in finite volumes”

## Conclusions

The conclusion/moral is self-evident. In general, given a function, we are always entitled to take its Legendre transform and express it in terms of the conjugate variables. In the context of equilibrium statistical mechanics, these transforms are invaluable because the equivalence of ensembles allows for the uninhibited two-way flow of information. For instance, if a particle system is studied at a fixed density then, except at points of phase transitions, we know everything about the fluctuating ensemble with the chemical potential adjusted to produce this density. Even more interesting—and even more useful—are the points of thermodynamic discontinuities. If one ensemble has a forbidden gap (say the particle density in the grand canonical distribution), then forcing the “parameter value” into the gap is essentially guaranteed to have interesting consequences in the other ensemble (e.g., phase separation).

But, the equivalence of ensembles is a mathematical—not to mention physical—fact only in the *thermodynamic limit*. In finite volume, as the droplet formation/dissolution phenomenon dramatically illustrates, the various ensembles are *not* equivalent. In these contexts, the assignment of physical—not to mention mathematical—significance to the conjugate variables is of dubious value. We suspect that this is the generic situation when “phase transitions” on a mesoscopic scale are the object of study. We believe that the dramatic inequivalence of ensembles in finite volume is the *signature* of interesting phenomena taking place below the macroscopic scale.

It is worth pointing out that, in the present context, the natural thermodynamic quantity which exhibits the signature jump is the good old *energy density*. There are several advantages to the use of this quantity as opposed to e.g. the magnetic field considered in . To list a couple, first, there is no numerical difficulty in the dynamical construction of the energy histogram and, second, there is no theoretical dispute in the interpretation of this quantity. Some previous efforts to exhibit the behavior of the energy density can be found in ; but, here we emphasize that the actual energy should be measured directly. Notwithstanding, if the physics of interest concerns *droplets*, it appears most natural to look for the droplet itself. This is evidently numerically feasible  and, presumably, permits the exhibition of *all* the secondary commodities.


text:
# Simulation software for ’Simulating Dynamical Features of Escape Panic’

## I Brief description

The source code  is in plain text format, and the files have been successfully compiled and run on various Linux systems. Documentation is available starting from the file 00\_README.txt .

The program has three different output formats.

* The primary output format is through images on the monitor (using the so-called X server).
* The secondary output format is Encapsulated Postscript images.
* Another possible output format is a compact data file containing all coordinates of all particles during the simulation.

Acknowledgements

The authors thank their host institutions for support, and Z. Csahók and Z. Farkas for permission to include the two library files they had previously developed.


text:
# Fusion cross sections at deep subbarrier energies

## Abstract

A recent publication reports that heavy-ion fusion cross sections at extreme subbarrier energies show a continuous change of their logarithmic slope with decreasing energy, resulting in a much steeper excitation function compared with theoretical predictions. We show that the energy dependence of this slope is partly due to the asymmetric shape of the Coulomb barrier, that is its deviation from a harmonic shape. We also point out that the large low-energy slope is consistent with the surprisingly large surface diffusenesses required to fit recent high-precision fusion data.

The primary ingredient in any nuclear reaction calculation is the nucleus-nucleus potential, consisting of the repulsive Coulomb interaction and an attractive nuclear part. Although the Coulomb term $`V_C(r)`$ is well-known, there are large ambiguities in the nucleus-nucleus potential $`V_n(r)`$, and many attempts have been made to extract information on this quantity from experimental data for heavy-ion reactions. While elastic and inelastic scattering are sensitive mainly to the surface region of the nuclear potential, the fusion reaction is also relatively sensitive to the inner part. They thus provide complementary sources of information.

In heavy-ion reactions, strong channel coupling effects (due to collective inelastic excitations of the colliding nuclei and/or transfer processes) significantly modify the landscape of potential energy surface, replacing the uncoupled single barrier with a distribution of barriers . In order to extract the nucleus-nucleus potential from heavy-ion fusion reactions, it is therefore advisable to use either high-energy fusion data where the barrier penetrability is essentially unity for all the distributed barriers, or very low-energy data where only the lowest barrier contributes to the cross section. Of these, the low-energy data probably provide cleaner information since the high-energy data may be complicated by competing reaction processes such as deep-inelastic scattering.

A recent paper  has reported an attempt to measure the fusion cross section $`\sigma `$ for the 60Ni + 89Y system at deep subbarrier energies, down to the 10-4 mb level. The authors of Ref.  used the Wong fusion formula  to analyse their data and showed that the experimental cross section exhibited an abrupt decrease at extreme subbarrier energies. They also analysed the data in terms of the logarithmic slope, defined by $`L(E)=d(\mathrm{ln}(\sigma E))/dE`$, and showed that this quantity exhibited a continuous increase with decreasing energy, in contrast to the theoretical slope which approached a constant value. They also found similar behavior in a few other systems found in the literature, including the 58Ni + 58Ni and 90Zr + 92Zr reactions.

The main part of the analysis in Ref.  relied on the Wong formula as a reference. A natural question is whether this formula, based on a parabolic approximation to the Coulomb barrier, is adequate at deep subbarrier energies . It was claimed in Ref.  that the Wong formula leads to fusion cross sections similar to those obtained with the coupled-channels approach for the 58Ni + 58Ni system. However, the former was simply a fit to the latter with parameters which had no physical connection to the potential used in the coupled-channels calculations.

The aim of this paper is to reanalyse critically the 58Ni \+ 58Ni reaction with an exact one-dimensional-potential calculation as well as with coupled-channels calculations  and show that the Wong formula is indeed unreliable at very low energies. This is particularly so for a quantity such as the logarithmic slope which accentuates the energy dependence of the cross section. We also discuss the findings of Ref.  in connection with the problem of the large surface diffusenesses of the nuclear potential for subbarrier fusion, discussed for some time in the literature .

Let us first discuss the validity of the parabolic approximation to the potential. Fig. 1 shows the nucleus-nucleus potential for the 58Ni \+ 58Ni system (solid line), along with its parabolic approximation (dashed line)

$$V(r)B\frac{1}{2}\mu \omega ^2(rR)^2,$$
(1)
where $`B`$ and $`R`$ are the barrier height and position, respectively. Here $`\mu `$ is the reduced mass of the system and $`\omega `$ is the barrier “curvature” given by $`\omega ^2=V^{\prime \prime }(R)/\mu `$. We use a Woods-Saxon nuclear potential with $`V_0=`$ 160 MeV, $`r_0`$ = 1.1 fm and $`a`$= 0.65 fm. On the inside, the nuclear potential varies relatively rapidly, while on the outside the Coulomb potential varies slowly, resulting in an asymmetric barrier shape. The deviation from the parabolic approximation (1) becomes larger as the energy goes down and one expects this approximation to break down at energies well below the barrier. It was shown in Ref.  that the parabolic approximation is adequate only for $`|rR|a`$, that is for incident energies within $`\mu \omega ^2a^2/2`$ of the barrier height. In the present example $`\mu \omega ^2a^2/2`$ = 2.62 MeV, and it is evident that the parabolic approximation is valid only in a relatively small range of energies near the barrier top.

An analytic formula for the fusion cross section for the parabolic barrier (1) was derived some time ago by Wong :

$$\sigma (E)=\frac{\mathrm{}\omega }{2E}R^2\mathrm{ln}\left[1+e^{2\pi (EB)/\mathrm{}\omega }\right].$$
(2)
The upper panel of Fig. 2 compares this formula with the fusion cross section obtained by numerically solving the Schrödinger equation with the true potential. No coupling is included in these calculations. As we saw in Fig. 1, the parabolic approximation underestimates the barrier thickness in the tunnelling region, and thus overestimates the penetrability at low energies. The bottom panel of Fig. 2 shows the logarithmic slopes $`L(E)`$. Eq. (2) yields a slope that is constant at low energies and is given by

$$L(E)\frac{d}{dE}\mathrm{ln}\left[\frac{\mathrm{}\omega }{2}R^2e^{2\pi (EB)/\mathrm{}\omega }\right]=\frac{2\pi }{\mathrm{}\omega }.$$
(3)
On the other hand, the slope computed from the exact results shows a continuous increase with decreasing incident energy (solid line). This is reminiscent of the experimental findings of Ref. .

At low energies, the logarithmic slope is related to the s-wave barrier penetrability $`P_0`$ by $`L(E)=d\mathrm{ln}(P_0(E))/dE`$. In the WKB approximation, the penetrability is given by

$$P_0(E)=e^{2S(E)/\mathrm{}}=\mathrm{exp}\left[2_{r_1}^{r_2}𝑑r\sqrt{2\mu (V(r)E)/\mathrm{}^2}\right]$$
(4)
at energies well below the barrier. Here, $`r_1`$ and $`r_2`$ are the inner and the outer turning points, respectively. Defining $`\mathrm{\Delta }(E)`$ as the difference between the true action integral $`S(E)`$ and its value in the quadratic approximation we have (ingoring an unimportant constant factor):

$$S(E)=_{r_1}^{r_2}𝑑r\sqrt{2\mu (V(r)E)}=\frac{\pi }{\omega }(BE)+\frac{\mathrm{}}{2}\mathrm{\Delta }(E).$$
(5)
Since the Coulomb barrier $`V(r)`$ has a non-symmetric shape, $`\mathrm{\Delta }(E)`$ increases as the energy decreases and the logarithmic slope $`L(E)=2\pi /\mathrm{}\omega d\mathrm{\Delta }(E)/dE`$ is always larger than $`2\pi /\mathrm{}\omega `$. Furthermore, one can show that the second derivative of this action integral is a positive quantity and thus $`L(E)`$ is a decreasing function of $`E`$. For example, this is the case for the sharp-cut potential, $`V_n(r)=\left[V_C(r)V_0\right]\theta (R_0r)`$, for which the action integral can be evaluated analytically . These facts are consistent with the numerical result shown in the lower panel of Fig. 2 as well as with the experimental findings discussed in Ref. . We thus conclude that the continuous increase of the logarithmic slope with decreasing energy is not in itself evidence of anomalous behavior of the fusion cross section at very low energies, as claimed in Ref..

We now discuss the relation between the logarithmic slope $`L(E)`$ and the surface property of the nuclear potential. For scattering processes, it seems well accepted that the surface diffuseness parameter $`a`$ should be around 0.63 fm if $`V_n`$ is parametrised by a Woods-Saxon form . In marked contrast, recent high-precision fusion data suggest that a much larger diffuseness, between 0.8 and 1.4 fm, is needed to fit the data . This is not just for particular systems but seems to be a rather general result. Note that fusion depends strongly on the potential on both sides of the barrier, in contrast to the elastic scattering which depends mainly on the potential on the outside. At high energies, the fusion cross section changes with the diffuseness due to the way the position and height of the $`l`$-dependent barrier change with increasing $`l`$. At lower energies, the main effect comes from the overall width of the barrier. A large diffuseness seems to be derirable in both these respects .

For a fixed value of the barrier height $`B`$, the barrier curvature $`\mathrm{}\omega `$ is approximately proportional to $`a^{1/2}`$ . Eq. (3) then indicates that the logarithmic slope $`L(E)`$ is roughly proportional to $`a^{1/2}`$. The large experimental slope found in Ref.  may therefore be another indication of the large surface diffusenesses already noted in heavy-ion fusion. In order to assess this, we perform the exact coupled-channels calculations for the 58Ni + 58Ni reaction using the computer code CCFULL  with different values of the surface diffuseness. This code uses the isocentrifugal approximation to reduce the dimensionality of the coupled-channels equations (see Ref.  for details) but we have checked that this is still valid at energies well below the Coulomb barrier. In the calculations, we include the double quadrupole-phonon excitations in both the projectile and target nuclei. A similar coupling scheme successfully explained the experimental fusion cross section and barrier distribution for the very similar 58Ni + 60Ni system . The dynamical quadrupole deformation parameter $`\beta _2`$ for the Coulomb coupling is estimated to be 0.177 from the experimental $`B(E2)`$  with the radius parameter $`r_{\mathrm{coup}}=1.2`$ fm. We require a somewhat larger value of $`\beta _2`$ = 0.261 (with $`r_{\mathrm{coup}}=1.06`$ fm) for the nuclear coupling in order to fit the data. The fusion reaction often requires a radius parameter of around 1.06 fm, smaller than the usual value of around 1.2 fm, used to extract a deformation parameter from the electromagnetic transition probability. This results in a larger deformation parameter as well as in a larger deformation length $`\beta r_{\mathrm{coup}}`$. Although the Coulomb-coupling hamiltonian is independent of the value of the radius parameter to be used, the nuclear coupling term depends on it through the combination $`\beta r_{\mathrm{coup}}`$. Therefore, this problem may also be related to the parametrization of the nucleus-nucleus potential and thus to the large surface diffuseness problem, though the value of $`r_{\mathrm{coup}}`$ = 1.06 fm should be reasonable for finite nuclei with a diffuse surface .

In Fig. 3, we show the dependence of the fusion cross section (upper panel) and of the logarithmic slope (lower panel) on the surface diffuseness parameter $`a`$ for the 58Ni + 58Ni reaction. The figure also includes the experimental data  for comparison. The experimental slope was computed using point-difference formulae with both two and three successive data points. The dotted line is the result with the nuclear potential shown in Fig. 1, that is with $`a`$=0.65 fm, while the dashed line is obtained with the potential parameters $`V_0=`$ 195 MeV, $`r_0`$ = 0.94 fm and $`a`$= 1.0 fm. The former leads to a cross section whose slope is not steep enough to account for the experimental data at energies below the barrier. As a consequence, the logarithmic slope $`L(E)`$ is underestimated at these energies, as in Ref. . On the other hand, the potential with $`a`$=1.0 fm improves the agreement considerably both for the cross section and the logarithmic slope. We also include in the figure a calculation with $`a=1.3`$ fm (solid line). This further improves the fit to the logarithmic slope, although it somewhat worsens the fit to the cross section itself at incident energies around 97 MeV. (We have confirmed that none of these results depends on the value of $`r_0`$ as long as $`V_0`$ is adujusted so that the barrier height remains unchanged.) Clearly, the experimental data favor a large value of the surface diffuseness, as in many other systems in the literature.

In summary, the “unexpected” behavior of heavy-ion fusion cross sections at extreme subbarrier energies claimed in Ref.  has two causes. One is the use of the Wong formula which is inadequate at energies far below the barrier. The exact numerical calculation is vital in discussing the fusion cross section and especially the logarithmic slope $`L(E)`$ at low energies. We pointed out that the exact calculation shows a similar energy dependence of the logarithmic slope as in the experimental data even without coupling. The other reason for this apparent anomaly is the use of a diffuseness parameter which is widely used in calculations for scattering processes, that is $`a`$ 0.63 fm. This potential leads to fusion cross sections whose logarithmic slope is much smaller than for the experimental data at deep subbarrier energies. If such a calculation is used as a reference, the experimental data may appear to fall much more steeply than expected . However, if one uses a larger value of the diffuseness parameter in the phenomenological potential, the data can be reproduced within the present coupled-channels framework. The need for a large diffuseness to describe the fusion process has also been found consistently in other systems. However, the reason for the large differences in diffuseness parameters extracted from scattering and from fusion analyses remains an open problem. In particular, it is still not clear whether a large surface diffuseness reflects the true nature of the potential or simply mocks up other effects which cause a rapid decrease of fusion at low energies. In this context, we mention that neither the double-folding potential  (which is usually much deeper and narrower than the Woods-Saxon) nor the geometrical corrections to the coupling potential  seem to resolve this problem.

More experimental and theoretical studies of fusion at deep subbarrier energies are needed to improve our understanding of this process which may be especially important in astrophysical fusion reactions. Isotopic dependences may also be of interest, particularly for exotic nuclei whose surface properties may be modified by the presence of weakly bound nucleons.

K.H. thanks the Nuclear Theory Group at the IPN Orsay for its hospitality, and the Kyoto University Foundation for financial support. M.D. acknowledges financial support from the Australian Research Council.


text:
# The spectral properties of the Falicov-Kimball model in the weak-coupling limit

## 1 Introduction

The Falicov-Kimball model (FKM) has become, since its introduction in 1969, one of the most popular examples of a system of interacting electrons with short-range interactions. The model was originally proposed to describe metal-insulator transitions  and has since been investigated in connection with a variety of problems such as binary alloys , the formation of ionic crystals , and ordering in mixed-valence systems . It is the latter language we shall use here, considering a system of localized $`f`$ electrons and itinerant $`d`$ electrons coupled via the local interaction $`U`$ and hybridization $`V`$. The Hamiltonian of the spinless FKM is

$$H=\underset{ij}{}t_{ij}d_i^+d_j+U\underset{i}{}f_i^+f_id_i^+d_i+E_f\underset{i}{}f_i^+f_i+V\underset{i}{}d_i^+f_i+h.c.,$$
(1)
where $`f_i^+`$, $`f_i`$ are the creation and annihilation operators for an electron in the localized state at lattice site $`i`$ with binding energy $`E_f`$ and $`d_i^+`$, $`d_i`$ are the creation and annihilation operators of the itinerant spinless electrons in the $`d`$-band Wannier state at site $`i`$.

The first term of (1) is the kinetic energy corresponding to quantum-mechanical hopping of the itinerant $`d`$ electrons between sites $`i`$ and $`j`$. These intersite hopping transitions are described by the matrix elements $`t_{ij}`$, which are $`t`$ if $`i`$ and $`j`$ are the nearest neighbors and zero otherwise (in the following all parameters are measured in units of $`t`$). The second term represents the on-site Coulomb interaction between the $`d`$-band electrons with density $`n_d=\frac{1}{L}_id_i^+d_i`$ and the localized $`f`$ electrons with density $`n_f=\frac{1}{L}_if_i^+f_i`$, where $`L`$ is the number of lattice sites. The third term stands for the localized $`f`$ electrons whose sharp energy level is $`E_f`$. The last term represents the hybridization between the itinerant and localized states.

In this paper we perform exhaustive numerical studies of the spectral properties of the FKM with and without hybridization. While the static properties of the FKM are well understood at present  (including the picture of valence transitions), the dynamical properties of the model are still unclear. Even, the spectral properties of the $`f`$ electrons are not understood satisfactorily nor for $`V=0`$, where only a few exact results are known for the infinite-dimensional systems . No exact results are known for nonzero hybridization and $`T=0`$, with the exception of numerical results obtained on very small clusters . The first systematic study of dynamical properties of the FKM with hybridization has been performed recently by Craco  within so called static approximation. He studied the model in the strong-coupling limit (large $`U`$) and found that the system is insulating for small values of hybridization (the size of the $`f`$ and $`d`$-electron gap coincides, including the case $`V=0`$) and with increasing $`V`$ an insulator-metal transition takes place at some critical value of $`V=V_c`$. In our preceding paper we reexamined these strong-coupling results by exact-diagonalization calculations and obtained fully different conclusions . In particular, we have found that for zero hybridization the gaps in the $`f`$ and $`d`$ electron density of states do not coincide, and almost all $`f`$ electron spectral weight is located outside the $`d`$ electron subbands. For nonzero hybridization the $`f`$ and $`d`$ electron gaps coincide, for both the symmetric and unsymmetric case, but no insulator-metal transition driven by hybridization is observed in one as well as in two dimensions. In this contribution we extend numerical calculations to the opposite limit (the weak-coupling limit). From this point of view the paper represents the first systematic (exact) study of dynamical properties of the FKM in the weak-coupling limit. Here, the special attention is devoted to the behavior of the $`f`$ and $`d`$-electron density of states with increasing $`E_f`$. Such an analysis is very important since a parametrization of $`E_f`$ with applied pressure $`p`$  can, in principle, provide an interpretation of some experimental data, e.g., the behavior of the activation gap (the gap between the occupied and unoccupied states) with increasing $`p`$. In our previous paper  the problem of pressure dependence of the activation gap was analysed through the behavior of the energy gap in the $`d`$-electron spectrum of the FKM and many similarities with the experimental data on the valence fluctuating compound SmB6 were found. These similarities are however only qualitative, since in the correct analysis one should take into account also the behavior of the $`f`$-electron spectral functions. This was done by Park and Hong  and a nice correspondence of theoretical and experimental results has been obtained. Unfortunately, these results have been obtained on a very small cluster (consisting of only eight sites) and thus cannot be considered as definite. Here we reexamine the behavior of the activation gap in the FKM for both $`V=0`$ and $`V>0`$. At zero hybridization we are able to present results on relatively large clusters for both the $`d`$ and $`f`$ electron density of states, and thus our results can be extrapolated satisfactorily to the thermodynamic limit. For nonzero hybridization we were able to perform small-cluster exact-diagonalization calculations on lattices only slightly larger than used by Park and Hong  ($`L12`$), however a fundamental different behavior of the model is observed already on such small clusters.

## 2 Results and discussion

Let us start the discussion of our results with the case $`V=0`$. In Fig. 1 we present numerical results for the $`d`$ and $`f`$ electron density of states obtained for $`U=0.6`$ and several values of the $`f`$-level energy $`E_f`$. Since the ground state configurations of the FKM in the weak coupling limit are well-known  (the most homogeneous (insulating) configurations for $`|E_f|<E_{fc}1.34`$, and the phase separated (metallic) configurations for $`|E_f|E_{fc}`$), the $`d`$-electron density of states can be calculated directly from the single particle spectrum of the FKM model for $`V=0`$. For the periodic configurations this can be done analytically in the thermodynamic limit , and for an arbitrary $`f`$-electron concentration a numerical diagonalization is possible on very large clusters ($`L=64000`$). Unlike this case, the $`f`$-electron density of states has to be calculated by exact diagonalization calculations (Lanczos method ), thereby the cluster sizes are strongly limited ($`L24`$). To minimize the finite-size effects on the $`f`$-electron spectra we have performed numerical calculations for several cluster sizes at each value of $`E_f`$. The typical behaviors are shown in Fig. 2 for two selected values of $`E_f`$. It is seen that finite-size effects are small and thus already results obtained for $`L=24`$ can be used satisfactorily to represent the behavior of large systems. Although we have used in the numerical calculations the relatively large value of $`ϵ=0.01`$ for the resolution , the formation of the gap in the $`f`$-electron density of states at the Fermi level (in all examined cases the Fermi level is located between the first and the second peak) is apparent. It is also apparent that the gaps in the $`d`$ and $`f`$ electron spectra do not coincide. Moreover, one can see (in accordance with results obtained in the strong-coupling limit ) that practically all $`f`$-electron spectral weight is located inside the principal $`d`$-electron gap (the gap at the Fermi energy). With increasing $`E_f`$ both the principal gap as well as the $`f`$-electron spectrum shift to higher energies. Since the finite size effects on the $`f`$-electron spectra are negligible for the cluster sizes used in our numerical calculations ($`L=24`$), and even the $`d`$-electron spectra can be obtained exactly on much larger clusters, one can try to construct the $`E_f`$ dependence of all relevant energy gaps. In particular, we have calculated the $`d`$-electron gap $`\mathrm{\Delta }_d`$ (the gap in the $`d`$-electron spectrum at the Fermi energy), the $`f`$-electron gap $`\mathrm{\Delta }_f`$ (the gap in the $`f`$-electron spectrum at the Fermi energy), and the $`fd`$-electron gap $`\mathrm{\Delta }_{fd}`$ (the gap between the occupied $`f`$ subband and the empty $`d`$ subband (the activation gap). To obtain the $`d`$-electron gap $`\mathrm{\Delta }_d`$ it is sufficient to know the $`f`$-electron distribution that minimizes the ground-state energy of the FKM for given $`E_f`$ and $`U`$. As mentioned above, the ground-state configurations of the FKM in the weak-coupling limit are the most homogeneous configurations for $`|E_f|<E_{fc}`$, and the phase separated configurations for $`|E_f|E_{fc}`$. While the most homogeneous distribution of $`f`$ electrons can be easily generated for arbitrary $`E_f(n_f)`$, the $`f`$-electron distribution in the phase separated configuration has to be determined numerically. This distribution can be found in principle exactly, as shown in  or approximately as done in . Here we adopt the latter method, since it allows to treat several times larger clusters ($`L300`$) and still to keep the high accuracy of computations. The results of numerical calculations for $`\mathrm{\Delta }_d`$ obtained using this method on the cluster consisting of 240 sites are presented in Fig. 3. It is seen that $`E_f`$ dependence of the $`d`$ electron gap qualitatively mimics the pressure dependence of the activation gap in SmB6  (we note a parametrization of $`E_f`$ with pressure). With increasing pressure (the $`f`$-level position $`E_f`$) the gap decreases and vanishes discontinuously at some critical pressure $`p_c`$ ($`E_{fc}`$). The finite-size scaling analysis that we have performed for a wide range of $`L`$ values ($`L=40,60\mathrm{}240`$) showed that $`\mathrm{\Delta }_d`$ is practically independent of $`L`$ for all $`|E_f|<E_{fc}`$. On the other hand the same analysis performed for $`|E_f|E_{fc}`$ revealed that the finite-size effects are still present in this region, but a disappearance of $`\mathrm{\Delta }_d`$ can be satisfactorily confirmed by extrapolation of results obtained for different $`L`$.

To determine the $`E_f`$ dependence of the $`f`$-electron gap $`\mathrm{\Delta }_f`$ , and the $`fd`$-electron gap $`\mathrm{\Delta }_{fd}`$ we have performed exhaustive numerical studies of the $`f`$-electron density of states for a wide range of $`E_f`$ values. In particular, we have calculated the $`f`$-electron density of states for each $`E_f`$ from the interval \[0,1.4\] with the step $`\mathrm{\Delta }E_f=0.05`$, and with the step $`\mathrm{\Delta }E_f=0.01`$ for $`E_f`$ near the point of the insulator-metal transition. To reveal the finite-size effects on $`\mathrm{\Delta }_{fd}`$ and $`\mathrm{\Delta }_f`$ the calculations have been done for several cluster-sizes $`(L=12,16,20,24)`$ at each $`E_f`$. It was found that finite-size effects on $`\mathrm{\Delta }_f`$ and $`\mathrm{\Delta }_{fd}`$ are negligible for $`|E_f|<E_{fc}`$, while small finite-size effects have been observed for $`E_f>E_{fc}`$. The resultant $`E_f`$ dependence of $`\mathrm{\Delta }_f`$ and $`\mathrm{\Delta }_{fd}`$ is presented in Fig. 3. It is seen that $`\mathrm{\Delta }_{fd}`$ exhibits the same behavior as $`\mathrm{\Delta }_d`$, of course with one exception and namely that $`\mathrm{\Delta }_{df}\mathrm{\Delta }_d/2`$. Since $`\mathrm{\Delta }_{fd}`$ is the gap between the occupied $`(f)`$ and unoccupied $`(d)`$ states (the activation gap) its behavior can be directly compared with the behavior of the activation gap in SmB6 . Although data for $`\mathrm{\Delta }_{fd}`$ are more scattered, both activation gaps exhibit qualitatively the same behavior. This indicates that the spinless FKM, probably the simplest model of correlated electrons can, in principle, provide the correct physics for describing properties of real materials.

Let us now examine what happens if the hybridization is switched on. One could expect that just the hybridization will improve the accordance between the theoretical and experimental results for the activation gaps, since the hybridization generally smears behaviors. This conjecture support also results by Park and Hong  obtained for the FKM with hybridization on a very small cluster consisting of only eight sites. Their results obtained for the $`E_f`$ dependence of the activation gap and the $`d`$-electron concentration $`n_d`$ are summarized in Fig. 4. It is seen a nice correspondence of theoretical and experimental data for both the activation gap as well as the $`d`$-electron concentration. Unfortunately, these results have been obtained on a very small cluster and so cannot be considered as definite. Even, the actual behavior of the activation gap and the $`d`$-electron concentration $`n_d`$ on large lattices can be fully different from one presented in Fig. 4. Indeed, a comparison of their results obtained for $`n_d`$ at $`V=0`$ with our results presented in Fig. 3b reveals fully different behavior of $`n_d`$. While $`n_d`$ calculated by Park and Hong  for $`L=8`$ is constant for all values of $`E_f`$ from $`E_f=0.8`$ to the insulator-metal transition point $`E_{fc}=1.34`$, our results calculated on much larger clusters ($`L=240`$) exhibit strong dependence on $`E_f`$. In the next we reexamine exactly the behavior of the FKM with hybridization on the cluster consisting of eight sites, as well as on clusters slightly larger ($`L=10,12`$). We present new results that significantly improve results obtained by Park and Hong

To show the hybridization effects on $`\mathrm{\Delta }_d`$,$`\mathrm{\Delta }_f`$ and $`\mathrm{\Delta }_{fd}`$ we have calculated exactly the $`f`$ and $`d`$ electron density of states for $`V=0`$ and $`V=0.02`$ on a small finite cluster consisting of $`L=8`$ sites. We chose this cluster size to be compatible with results of Park and Hong . The resultant behaviors are presented in Fig. 5. The most prominent difference between the results obtained for $`V=0`$ and $`V=0.02`$ is that a nonzero spectral weight appeared in the $`d`$-electron density of states at the Fermi level for finite hybridization. Of course, this fact will change dramatically the picture discussed above for $`V=0`$. This is clearly demonstrates in Fig. 6a, where the $`E_f`$ dependence of $`\mathrm{\Delta }_d`$, $`\mathrm{\Delta }_f`$ and $`\mathrm{\Delta }_{fd}`$ is displayed for $`V=0.02`$ and $`L=8`$. It is seen that all gaps coincide for nonzero $`V`$, what strongly contradicts to the case of $`V=0`$. For a comparison we have displayed in Fig. 6a also the behavior of the single particle excitation energy defined as $`\mathrm{\Delta }_s=E_G(L+1)+E_G(L1)2E_G(L)`$, where $`E_G(N)`$ is the ground state energy for $`N`$ electrons. As one could expect (on the base of the strong-coupling results ) the single particle excitation energy $`\mathrm{\Delta }_s`$ coincides with $`\mathrm{\Delta }_d`$, $`\mathrm{\Delta }_f`$ and $`\mathrm{\Delta }_{fd}`$ for nonzero hybridization. This result is very important from the numerical point of view since the single particle excitation energy can be calculated easily by other methods (that allow to treat much larger clusters, e.g., the density matrix renormalization group method (DMRG)), and so it can be used satisfactorily for describing conducting properties of the model. As one can see from Fig. 6a the conducting properties of the FKM with hybridization are described very inaccurate within the small-cluster exact-diagonalization calculations. In the region where the metallic phase has been identified for $`V=0`$, the single particle excitation energy increases and it is not clear if this is a consequence of finite hybridization or finite-size effects. Since calculations performed for $`L=10`$ and $`L=12`$ revealed relatively large finite-size effects we have decided to use DMRG method to verify the actual behavior of $`\mathrm{\Delta }_s`$. In comparison to exact diagonalization calculations the DMRG method allows to treat several times larger clusters ($`L50`$) and still to keep the high accuracy of computations. We typically keep up to 128 states per block, although in the numerically more difficult cases, where the DMRG results converge slower, we keep up to 400 states. Truncation errors , given by the sum of the density matrix eigenvalues of the discarded states, vary from $`10^8`$ in the worse cases to zero in the best cases. The typical behavior of $`\mathrm{\Delta }_s`$ is displayed in Fig. 6b for two finite clusters of $`L=32`$ and $`L=36`$ sites. Although the finite-size effects are still present for cluster sizes treated by DMRG method the metallic region seems to be satisfactorily verified also for nonzero hybridization. However, there is one important difference between the case $`V=0`$ and the case $`V=0.02`$. While for $`V=0`$ both $`\mathrm{\Delta }_d`$ and $`\mathrm{\Delta }_{fd}`$ vanish discontinuously, they seem to vanish continuously for nonzero hybridization. This indicates that the spinless FKM behaves fully differently for zero and non-zero hybridization and fully different are also the corresponding pictures of insulator-metal transitions. Since rare-earth compounds exhibit both types of insulator-metal transitions (discontinuous as well as continuous) their different behavior under the external pressure can be interpreted directly as a consequence of absence (presence) of hybridization in a given material. According to this conjecture, e.g., SmB6 that exhibits a discontinuous pressure-induced insulator-metal transition should be a material with zero hybridization between $`d`$ and $`f`$ states. Indeed, SmB6 is the highest symmetry of the $`O_h`$ point group and the on-site hybridization between $`d`$ and $`f`$ states is forbidden by inversion symmetry . This confirms that the spinless FKM, in spite of its relatively simplicity, can yield the correct physics for a description of real materials, e.g. SmB6. Besides of a qualitative correspondence between the theoretical and experimental data for the activation gap, we can state another example that supports this conclusion. In Fig. 7 we have displayed in detail the $`d`$ and $`f`$ electron density of states near the Fermi level for several values of $`E_f`$ ($`E_f=0,0.1,0.2`$). For this set of $`E_f`$ values the ground state is characterized by $`n_f=0.5`$ what models very well the situation in SmB6 at low temperatures . This electronic structure consisting of two wide and two narrow subbands strongly mimics the electronic structure that we have proposed several years ago for a description of low temperature resistivity data of SmB6 . The analysis of these data showed that they can be explained by introducing a fine structure (consisting of two narrow bands) into the principal gap, in accordance with our results presented in Fig. 7. Of course, all these comparisons are only qualitative since our results have been obtained for the one-dimensional system at $`T=0`$, while the real experimental systems are three dimensional and measurements are done at finite temperatures. In future work, we plan to perform a similar analysis in higher dimensions and $`T0`$. Moreover, we also want to examine the influence of other factors (the electron-phonon interaction, the orbital dynamics, etc.) that have been neglected in this version of the model.

In summary, the $`f`$ and $`d`$ electron density of states of the spinless Falicov-Kimball model were studied in the weak-coupling limit by exact diagonalization calculations. The resultant behaviors were used to examine the $`d`$-electron gap $`\mathrm{\Delta }_d`$, the $`f`$-electron gap $`\mathrm{\Delta }_f`$, and the $`fd`$-electron gap $`\mathrm{\Delta }_{fd}`$ as functions of the $`f`$-level energy $`E_f`$ and hybridization $`V`$. It was shown that the spinless Falicov-Kimball model behaves fully differently for zero and finite hybridization between $`f`$ and $`d`$ states. At zero hybridization the energy gaps do not coincide ($`\mathrm{\Delta }_d\mathrm{\Delta }_f\mathrm{\Delta }_{fd}`$), and the activation gap $`\mathrm{\Delta }_{fd}`$ vanishes discontinuously at some critical value of the $`f`$-level energy $`E_{fc}`$. In all examined cases practically all $`f`$-electron spectral weight is located outside the $`d`$-electron subbands. On the other hand, at finite hybridization all energy gaps coincide and vanish continuously at the insulator-metal transition point $`E_f=E_{fc}`$. The results obtained are compared with experimental behaviors observed in SmB6.

This work was supported by the Slovak Grant Agency VEGA under Grant No. 2/7021/20 and the Science and Technology Assistance Agency under Grant APVT-51-021602. Numerical results were obtained using the PC-Farm of the Slovak Academy of Sciences.

Figure Captions

Fig. 1. The $`f`$ and $`d`$ electron density of states of the FKM obtained for $`V=0,U=0.6`$ and several different values of $`E_f`$. The results for the $`f`$ electron density of states are plotted for $`L=24`$, while the results for the $`d`$ electron density of states are plotted for $`L=64000`$. The corresponding $`f`$-electron concentrations are: $`n_f=0.5`$ for $`E_f=0`$, $`n_f=4/10`$ for $`E_f=0.5`$, $`n_f=1/3`$ for $`E_f=0.8`$, and $`n_f=1/4`$ for $`E_f=1.2`$.

Fig. 2. The $`f`$ electron density of states of the FKM obtained for $`V=0,U=0.6`$ and two different values of $`L`$ and $`E_f`$.

Fig. 3. The $`d`$-electron gap $`\mathrm{\Delta }_d`$, the $`f`$-electron gap $`\mathrm{\Delta }_f`$, the $`fd`$-electron gap $`\mathrm{\Delta }_{fd}`$, and the $`d`$-electron concentration $`n_d`$ as functions of the $`f`$-level energy $`E_f`$ calculated for $`V=0`$ and $`U=0.6`$. The solid lines are guides to the eye.

Fig. 4. The activation gap $`\mathrm{\Delta }_{fd}`$ and the $`d`$-electron concentration $`n_d`$ as functions of the $`f`$-level energy $`E_f`$ calculated for $`U=0.6`$ and $`L=8`$ (Ref. 7).

Fig. 5. The $`f`$ and $`d`$ electron density of states of the FKM obtained for $`U=0.6,E_f=1.2,L=8`$ and two different values of $`V`$. The arrow denotes the most prominent difference between the results obtained for $`V=0`$ and $`V=0.02`$.

Fig. 6. a) The $`d`$-electron gap $`\mathrm{\Delta }_d`$, the $`f`$-electron gap $`\mathrm{\Delta }_f`$, the $`fd`$-electron gap $`\mathrm{\Delta }_{fd}`$, and the single particle excitation energy $`\mathrm{\Delta }_s`$ as functions of the $`f`$-level energy $`E_f`$ calculated for $`V=0.02,L=8`$ and $`U=0.6`$. b) The single particle excitation energy $`\mathrm{\Delta }_s`$ as a function of the $`f`$-level energy $`E_f`$ calculated by DMRG method for the same values of $`V`$ and $`U`$. The solid and dashed lines are guides to the eye.

Fig. 7. The $`f`$ and $`d`$ electron density of states of the FKM obtained for $`V=0,U=0.6`$ and several different values of $`E_f`$. The results for the $`f`$ electron density of states are plotted for $`L=24`$, while ones for the $`d`$ electron density of states are plotted for $`L=64000`$. In all cases the $`f`$-electron concentration is equal to $`n_f=0.5`$. To see clearly the $`f`$-electron gap, the $`f`$-electron density of states has been plotted with high resolution ($`ϵ=0.001`$).


text:
# Non-cubic layered structure of Ba1-xKxBiO3 superconductor

## Abstract

Bismuthate superconductor Ba1-xKxBiO3 (x=0.27-0.49, Tc=25-32K) grown by an electrolysis technique was studied by electron diffraction and high-resolution electron microscopy. The crystalline structure thereof has been found to be non-cubic, non-centrosymmetric and of the layered nature, with the lattice parameters a$``$ap, c$``$2ap (ap is a simple cubic perovskite cell parameter) containing an ordered arrangement of barium and potassium. The evidence for the layered nature of the bismuthate superconductor removes the principal crystallographic contradiction between bismuthate and cuprate high-Tc superconductors.

Superconducting bismuthates including the firstly discovered BaPb1-xBixO3  and Ba1-xKxBiO3 (BKBO)  with the highest transition temperature (Tc=32-35K for x$``$0.4) have the basic characteristics similar to those of high-Tc copper-oxide superconductors . This might favor for a common pairing mechanism. On the other hand, bismuthates have been generally considered to be fundamentally different from cuprates due to two reasons. Bismuthates are non-magnetic, and they have three-dimensional structures rather than layered two-dimensional ones characteristic of cuprates. Of these two contradictions, the most important one seems to be the latter. The absence of magnetic fluctuations in bismuthates may point to a non-magnetic nature of the pairing mechanism, whereas it is unlikely that a common superconducting scenario in cuprates and bismuthates does not depend on lattice dimensionality.

The study of BKBO is of a particular interest also for understanding the relation between charge-density-waves (CDWs) and superconductivity in high-Tc oxides. The parent compound for BKBO is the perovskite BaBiO3 containing a CDW formed of an ordered arrangement of non-equivalent bismuth ions referred to as Bi3+ and Bi5+ . This CDW is assumed to be responsible for the semiconducting behavior of BaBiO3 and Ba1-xKxBiO3 materials with low potassium content (x $`<`$ 0.25). The widely known notion that BKBO superconductors have a simple cubic ABO3 solid-solution structure of a non-layered nature, with barium and potassium randomly occupying the A-position, was inferred from long-range structural studies of ceramic samples by X-ray  and neutron diffraction . A simple cubic structure excludes the existence of a CDW that leads to the conclusion of the total incompatibility of CDWs and superconductivity. However, studies of BKBO by methods sensitive to short-range symmetry, in particular, by Raman scattering spectroscopy , a paired-distribution function analysis of neutron diffraction data , extended X-ray fine structure analysis (EXAFS)  evidence that the local structure of BKBO superconductors is not cubic. An X-ray diffraction study  of Ba0.6K0.4BiO3 single crystals grown by an electrolysis technique  has revealed supercell reflections assigned to a CDW.

It has been hypothesized  that BKBO superconductors may have layered structures with an anisotropic ordered arrangement of potassium and barium. Detailed studies  of phase relations in the Ba-Bi-O system have revealed a series of BanBin+mOy oxides of ordered layered structures . BanBin+mOy are assumed to be matrices transformed, when intercalated with potassium, into superconducting oxides retaining the layered nature of the matrices structures. Besides, Ba-rich oxides with a solid-solution structure that attributed to BKBO superconductors were discovered . Such oxides are formed in two-phase regions of the Ba-Bi-O system. Assuming that the same phenomenon is inherent also in the Ba-K-Bi-O system, we looked for conditions of single-phase growth of BKBO crystals. By varying parameters of the electrolysis technique, we examined the growth process and obtained the data suggesting the existence of individual superconducting phases BanKmBin+mOy with different Tc (8-35K) .

This paper reports electron diffraction (ED) and high-resolution electron microscope (HREM) studies of superconducting crystals Ba1-xKxBiO3 (x=0.27-0.49, Tc=25-32K) grown by a modified electrolysis technique . We present the evidence obtained by a long-range structural method that the BKBO superconductor has a non-cubic structure with layered ordering of barium and potassium.

Superconducting crystals were produced by electrolysis of KOH-Ba(OH)2-Bi2O3 melt (K:Ba:Bi=72:1.33:2) with the current of 5.3 mA at 300oC for 5 h. The anode deposit was a polycrystalline boule consisting of intergrown single crystals. According to an X-ray powder diffraction analysis, the deposit contained several phases of pseudocubic perovskite structure (ap=0.4277-0.4310 nm). The temperature dependence of the magnetic susceptibility $`\chi `$ was measured for four crystals of the cubic shape of 0.5-1 mm3 in volume chosen in the deposit. The crystals displayed similar curves, $`\chi `$(T), showing bends at 25, 27, 30 and 32K that indicates the presence of different superconducting phases.

One of the four crystals was selected for electron microscope studies. It was ground to prepare a suspension with particles of a few $`\mu `$m in size which was deposited on holey carbon films. Electron diffraction studies were performed in an electron microscope JEOL JEM-2000FX equipped with a system for an energy dispersive X-ray (EDX) elemental analysis. High-resolution studies were performed in a microscope Hitachi HF-3000. In order to avoid beam-induced modulations , the experiments were run at minimal electron beam intensity. We examined 40 particles of the ground crystal, EDX spectra and ED patterns were taken from each particle. Although the full composition range of Ba1-xKxBiO3 particles was found to be x=0.27-0.49, majority of them had K:Ba:Bi ratio close to average value of 0.38:0.62:1 obtained by summing over all measurements.

Half of the particles exhibited diffraction patterns (fig.1) containing supercell reflections with the vector q=$`\frac{1}{2}`$ (indexing in terms of a simple cubic perovskite cell) in addition to the basic perovskite reflections. As only two out of three $``$100$``$ axial directions are observable at electron microscope studies (due to a limited range of specimen tilt), some of particles not displaying supercell reflections were, in fact, also of the ordered nature, but those particles were observed in the  zone axis when the supercell reflections could not be excited. Assuming that the predominant cleavage at crystal grounding is weak for BKBO superconductors, a lower limit of the relation of all ordered particles to non-ordered can be estimated from the above data as 3:1. It indicates that the bulk of the studied crystal was of the ordered state. From the temperature dependence of the magnetic susceptibility, which changed much and sharp at 30K and smoothly at lower temperatures, it follows that the ordered state is related to Tc=30K, whereas the non-ordered part of the crystal became superconducting at Tc=25-27K.

The anisotropic orientation of the supercell reflections along only one of the axial directions indicates that the ordered lattice is of non-cubic symmetry. As no clear splitting of perovskite spots was observed, which would indicate the presence of twins of a phase with lower crystallographic symmetry, it follows that tetragonal symmetry can be assigned to the ordered lattice. From the magnitude and type of the supercell reflections one may state that the supercell is primitive with the the lattice parameters a$``$ap and c$``$2ap, it consists of two perovskite blocks.

A  HREM image taken from the particle, which displayed supercell spots in ED patterns, is shown in fig.2a. The ordering manifests itself as intensity modulations along the  direction with the period of 2ap. A detailed picture of the image contrast is shown in fig.2b presenting the enlarged image of the area near the edge of the particle. It can be visualized that the supercell consists of two perovskite blocks. Matching the corresponding spots in the image and measurements of the distances between them disclose that the image contrast of the blocks is different, and the blocks are slightly different in size (4-5%). This indicates the the blocks in the supercell are not equivalent. An important feature of the HREM image in fig.2b can be noticed, namely, the absence of a symmetry center, which suggests a non-centrosymmetric nature of the supercell structure.

Two crystallographic models of a different nature may be responsible for the appearance of supercell reflections with the vector q=$`\frac{1}{2}`$: (1) the model of solid solution with the occurrence of common barium and potassium positions in the perovskite blocks, (2) the model of ordering of barium and potassium.

In the model of solid solution the appearance of the supercell reflections is related with distortions (tilting and/or deformation) of oxygen octahedra surrounding bismuth ions. This approach was applied to explain the origin of supercell reflections with the vectors q=$`\frac{1}{2}`$$``$111$``$ observed in diffraction patterns of the perovskite BaBiO3 . For our case, a scheme of one-dimensional deformation of oxygen octahedra along the  direction without tilting, when octahedra are alternatively different in size with the period of 2ap, can only give rise to reflections with q=$`\frac{1}{2}`$, because the apxapx2ap supercell is not compatible with the presence of tilting . In this scheme, there is a non-equivalency of bismuth ions in neighbouring planes rectangular to the  axis, which is interrelated with the existence of a CDW.

To establish if the solid solution model is credible, we performed HREM image simulations using the NCEMSS program  with variation of possible parameters affecting the image contrast: ion coordinates, specimen thickness, objective lens defocus value, specimen deviation from the exact zone orientation. The simulation has revealed that no asymmetry in HREM images of centrosymmetric cells is observed, whereas it does occur for non-centrosymmetric ones. However, a trustworthy agreement between simulated images of non-centrosymmetric cells and the experimental images cannot be achieved by varying the parameters. So, we conclude from structural studies that the model of solid solution should be rejected. The model of only oxygen octahedra deformations is not compatible also with superconducting properties of BKBO crystals of the compositions studied here. Quite intensive reflections with q=$`\frac{1}{2}`$ (more intensive than reflections with q=$`\frac{1}{2}`$$``$111$``$ characteristic of semiconducting BaBiO3) would imply large octahedra deformations and hence the existence of a very strong CDW, which would result in dielectric properties.

In the model of ordering one of the two cation positions between the oxygen octahedra (A-positions) is fully occupied by barium. The second A-position is a common site for barium and potassium with the highly prevailing amount of potassium (for Ba0.6K0.4BiO3 composition being close to the average one of investigated particles, the ratio of the occupation factors Ba:K is 1:4 in the second A-position). The image simulations (Fig.2c) have shown that the contrast of experimental images can be described in the model of ordering with non-symmetric ion coordinates. Intensity of the supercell reflections is dictated by a difference in scattering factors of Ba and K.

Fig.3 illustrates the model of the structure consistent with the HREM images. The tetragonal space group P4mm can be taken as a possible non-centrosymmetric group to describe the lattice symmetry. Barium and potassium are displaced against the centers of the perovskite blocks in parallel manner. Apical oxygens O(2) and O(4) are synchronously displaced from the centers of the edges, that results in distinct positions thereof relative to bismuth ions. Bismuth ions in neighbouring BiO2 planes with planar oxygens O(1), O(3) are non-equivalent in this model as well as in the solid solution one.

The authors of initial studies  of BKBO superconductors by EXAFS assumed the existence of only one Bi-O bond length in the structure. The recent EXAFS study  has found two distances of oxygen relative to bismuth. It can be assumed that the bond lengths of bismuth with planar oxygens O(1), O(3) (Fig.3) are close or the same, whereas the bond lengths of bismuth with apical oxygens are clearly different. As the number of apical bonds is smaller than that of planar ones, the contribution of the formers into the general set of Bi-O bonds is insignificant and gives rise to the second-order effect . However contrary to the assumption put forward in Ref.10 about the dynamic nature of the occurrence of different Bi-O bonds in BKBO superconductors, our study shows that it is of a static nature.

BKBO superconductors have the formal oxidation state of bismuth ions in the range of 4.35-4.55 . A part of bismuth ions may be present in the stable oxidation state +3. It is known that Bi3+ ions in the structures of bismuth-containing oxides are coordinated with surrounding oxygen ions asymmetrically due to the presence of stereoscopically active lone electron pairs of Bi3+ ions. It can be assumed that a collinear orientation of these electron pairs gives rise to the asymmetric structure in fig.3.

It is significant that the existence of non-equivalent metal-oxygen distances in the ordered structure of BKBO points to the likely presence of a CDW. This CDW is obviously other than the semiconducting CDW existing in the parent perovskite BaBiO3, because the non-equivalent Bi ions in BaBiO3 have all (six) different Bi-O bonds, whereas such ions in the ordered BKBO differ with respect to one or two apical oxygens. This is consistent with the X-ray study  giving an evidence that splitting of the supecell reflections resulting from the domain structure was sensitive to electric and magnetic fields. Further studies with the use of a quantitative diffraction technique (e.g. neutron diffraction) are needed to reveal the nature of the CDW and details thereof which may be related to the pairing mechanism.

An important point is that the ordered arrangement of barium and potassium ions along the  axis results in an anisotropic structure, where stacking of BiO2 planes with planar oxygens O(1),O(3) along the  axis (fig.3) endows the structure with the layered nature. This makes the crystalline structure of the ordered BKBO with the separate BiO2 planes qualitatively similar to crystalline structures of cuprate superconductors with CuO2 planes that cancels their principal crystallographic contrasting. The difference between them is that cuprates have layered structures with distinct anisotropy, whereas the layered BKBO has a weakly anisotropic structure due to close ion radia of Ba2+ (0.138 nm) and K+ (0.133 nm).

In summary, we have observed a non-cubic structure of the superconducting bismuthate Ba1-xKxBiO3 (x=0.27-0.49, Tc=25-32K). The structure is ordered with the lattice parameters a$``$ap, c$``$2ap in terms of a simple cubic perovskite cell parameter ap. It has been revealed that the ordering nature is related to the ordered arrangement of barium and potassium ions, which endows the structure of Ba1-xKxBiO3 with a layered character being typical for structures of superconducting cuprates. The crystalline cell of the superconductor has been found to be non-centrosymmetric with non-equivalent metal-oxygen distances that may point to the existence of a charge-density wave.

The work at Chernogolovka is supported in part by the Superconductivity section of the Russian State Research and Engineering program ”Topical Problems in the Physics of Condensed Matter” and the project no. 02-02-1678 of the RFBR.