text:
# Magnetically Stabilized Luminescent Excitations in Hexagonal Boron Nitride

## 1 Introduction

The study of luminescence in hexagonal boron nitride (h-BN) informs us about properties intrinsic to the material as well as its defects and impurities. Techniques used to study h-BN include electron paramagnetic resonance, thermally stimulated luminescence, thermally stimulated conduction, and luminescence induced by light, cathodic, electric, ionic, and radioactive sources. Identified defects include F-centers, nitrogen vacancies and elemental impurities . Information about the intrinsic nature of h-BN has been obtained through studies of the near-band-gap luminescence.

A continuum photon emission spectrum in h-BN has been observed as a result of thermal neutron exposure which produces energetic He+ (1.47 MeV) and 7Li+ (0.84 MeV) ions in the material via the 10B(n,$`\alpha `$)7Li reaction. Other irradiation studies used 1.0 MeV H+ and 1.8 MeV He+ ions and resulted in the characterization of a continuum photon emission spectrum of 280 nm to 600 nm. Here we present observations of magnetic suppression of cold neutron induced luminescence in h-BN at low temperature. The luminescence is stimulated by the by-products of neutron capture and can be suppressed by application of a magnetic field. The studies we present here were prompted by the use of h-BN as a neutron shielding material in an experiment to magnetically trap ultracold neutrons (UCN).

## 2 Experimental Apparatus

Hexagonal boron nitride is used in our neutron trapping apparatus to collimate an incident neutron beam and to absorb those neutrons scattered out of the collimated beam. Boron compounds are commonly used as shielding materials for absorbing low energy neutrons. Natural boron has a 19.9 % abundance of 10B, which has a large thermal neutron absorption cross section, $`\sigma _{\mathrm{thermal}}=3840`$ b. Hexagonal boron nitride is particularly suitable because of its high purity and mechanical properties. Hexagonal BN can be obtained with $`<0.4\%`$ oxygen and $`<0.1\%`$ all other total elemental impurities. It is easily machinable and has a low thermal expansion coefficient. The h-BN used in our experiment is obtained commercially as a solid cylinder and machined into the form of interlocking tubular pieces with inner diameters of 3.86 cm and outer diameters of 4.16 cm. The assembled h-BN tube fits inside a cupronickel (CuNi) cell body and extends its entire length, 110 cm (see figure 1). Two cylindrical rings are used as neutron collimators and are positioned at distances of 36 cm and 74 cm relative to the front (beam entrance) of the h-BN tube. These rings each have an inner diameter of 2.0 cm, length of 1 cm and fit snugly within the h-BN tube. The cell is filled with liquid 4He and is cooled by a dilution refrigerator to 0.1 K to 1.0 K.

A cold neutron beam (moderated to 40 K) traverses the cell parallel to the long axis with a fluence of approximately $`5\times 10^8\mathrm{s}^1`$. About 60 % of the neutrons entering the cell are absorbed by the h-BN. A set of superconducting magnets resides just outside the cell body and when at maximum current, the h-BN is in a magnetic field ranging from 1.0 T to 1.8 T.

Under normal operating conditions, a light collection system is situated just inside the tube of h-BN between the two collimators (see figure 1). This system consists of a diffuse reflector made of expanded polytetrafluoroethylene (PTFE) onto which is evaporated an ultraviolet-to-blue downconverter (tetraphenyl butadiene, TPB) and two helically coiled blue-to-green wavelength shifting fibers. The fibers exit the cryogenic apparatus and are each coupled to a photomultiplier tube. This system is sensitive to photons over the entire spectrum from the extreme ultraviolet (EUV; $`<100`$ nm) to the visible blue ($`500`$ nm) and is not capable of further wavelength discrimination111The fiber light collection system was not actually used in the measurement made in Ref. ..

## 3 Results

In a typical run, the detection region is exposed to neutrons for 900 s. During this irradiation period, the photomultiplier tubes (PMTs) are turned off. After 900 s, the neutron beam is blocked, the PMTs are brought up to voltage and luminescence is detected.

The time dependence of the luminescence signal is shown in figure 2 (solid line), starting at the point when the beam is turned off ($`t=0`$ s). This decay does not fit to a single exponential and varies roughly as inverse time. When the magnet is energized during the irradiation and observation, a significant fraction of this luminescence signal is suppressed (see figure 2, dashed line). When the field is lowered (starting at 1275 s), luminescence from the “stored” component is emitted.

Numerous tests were performed to identify the source of the luminescence signal and to characterize its behavior. By the process of elimination using multiple experimental setups, we determined that the luminescence signal was originating from the boron nitride and simplified the detection system to verify this conclusion. In the experimental setup shown in figure 3, only h-BN and an acrylic lightguide were present (the TPB downconverter and wavelength-shifting fibers were removed), and the geometry was arranged such that neutrons could not directly impact the acrylic. The luminescence signal was observed. Both components of this decay signal can be seen in figure 4. In separate tests performed at 4 K, we did not observe luminescence in an acrylic sample under similar neutron irradiation conditions.

Studies to characterize the behavior of the magnetically stabilized component of the luminescence signal were performed using the fiber cell apparatus (figure 1). Data resulting from lowering the magnetic field in various ways is shown in figure 5. The initial irradiation takes place with the magnets at maximum current. In figure 5A the magnet current is held constant for 800 s after the neutron beam is turned off. Then the magnetic field is reduced in four steps: from $`100\%75\%`$, $`75\%50\%`$, $`50\%25\%`$ and $`25\%0\%`$ of its maximum value. The largest fraction of the luminescence is emitted in the last two reduction steps. Figure 5B further demonstrates the suppression of luminescent emission. In this measurement, the magnetic field is slowly ramped down from 100 % to zero over 2700 s. Immediately after this, the field is rapidly ($`100`$ s) turned back to 100 % for 570 s and then back to zero. As can be seen, the luminescence is suppressed while the field is on. We have verified, in a separate experiment, that the magnetic field has no significant effect on the efficiency of the detection system.

To characterize the time dependence of the decay of the stabilized component, the superconducting magnet was quenched resulting in a reduction in field from its maximum value to zero in less than 1 s. The time dependence of this decay signal is roughly the same as the initial decay (see figure 6). In addition, the total integrated number of counts using the magnet quench signal is roughly equal to that obtained using a slow field ramp. Presumably the luminescence is the result of relaxation of some type of solid state excitation in the h-BN.

The persistence of these excitations under conditions of high magnetic field is most dramatically demonstrated by leaving the magnetic field on for 28,800 s after irradiation. After this time the magnetic field is lowered to zero and luminescence is observed. The number of photons emitted after 28,800 s was within 50 % of the number emitted after only 1000 s of storage (under the same irradiation conditions).

We also observed that raising the temperature of the cell has a similar effect to lowering the magnetic field. Figure 7 shows data for different field and temperature conditions. In both sets of data the cell was irradiated by the neutron beam for 900 s. In figure 7A the magnetic field is lowered from its initial value to zero in approximately 100 s beginning 1275 s after the beam is turned off. Traces are shown for initial magnet currents of 100 %, 75 %, 50 % and 25 % of its maximum current. The data shown in figure 7B is taken maintaining the magnet at its maximum current while raising the temperature of the cell from its initial value to 1 K in approximately 100 s beginning 1275 s after the beam is turned off. The traces are for initial temperatures of 170 mK, 400 mK, 525 mK and 650 mK. It is important to note that temperature changes brought about by eddy current heating while ramping the magnetic field ($`<20`$ mK) are much smaller than those required to initiate the release of luminescence ($`>200`$ mK).

## 4 Discussion

Given the myriad of mechanisms for producing luminescence, we can only offer general hypotheses about the magnetically stabilized component of the luminescence signal and offer some general points. The energy of interaction of the magnetic moment of an electron or hole in h-BN with a magnetic field of 1 T in temperature (Kelvin) units is 1.3 K. This is the same energy scale as that observed in figure 7. Our highest applied field of 1.8 T results in an energy of interaction $`2.3`$ K. The applied magnetic field should lead to polarization of paramagnetic species as the thermal fluctuations present in the solid are not energetic enough for depolarization. In addition to electrons and holes, other possibilities for paramagnetic species present in the h-BN include triplet state excitons. These typically have a magnetic moment of one Bohr magneton ($`1\mu _B`$) or greater. Spin polarization of trapped electrons, holes or excitons may take place and then affect the possible relaxation paths to luminescence.

One possibility is that the magnetic field stabilizes a triplet state of a defect center such as an F-center. Triplet states can have very long lifetimes and the stabilization would prevent conversion to the much shorter lived singlet state. Observations of magnetic suppression of singlet state formation, resulting in enhanced triplet state formation has been observed in alkali halides. Unlike our case, these triplet states exhibit pure exponential decay which is uneffected by magnetic field changes. Another possibility is that the interaction of an electron or exciton with a paramagnetic center causes trapping, similar to that seen in indium antimonide. The applied magnetic field could keep the two species aligned. As the temperature is raised (or the field lowered) thermal fluctuations would release the excitations which could then undergo radiative recombination. The apparent distribution of trapping energies may indicate a combination of thermoluminescent activity with a magnetic interaction. That is, the usual thermoluminescent trapping is modified by some type of magnetic interaction, perhaps one of those discussed above. In this case, the spatially varying magnetic field would provide a continuum of trapping energies.

## 5 Conclusion

We have observed a magnetically stabilized component in the luminescence of hexagonal boron nitride under neutron irradiation. The excitations that give rise to the luminescence can be stabilized for long periods by high magnetic fields ($`1`$ T) at low temperatures ($`600`$ mK). Further experiments are necessary to elucidate the nature of the excitations, including spectroscopic analysis of the luminescence, quantitative analysis of the magnetic field dependence and optically detected electron paramagnetic resonance.

## 6 Acknowledgments

This work was supported in part by the National Science Foundation under grant No. PHY-9424278.


format: text

text:
# Experimental Evidence for a Spin-Polarized Ground State in the 𝜈=5/2 Fractional Quantum Hall Effect

\[

## Abstract

We study the $`\nu =5/2`$ even-denominator fractional quantum Hall effect (FQHE) over a wide range of magnetic field in a heterojunction insulated gate field-effect transistor (HIGFET). The electron density can be tuned from $`n=0`$ to $`7.6\times 10^{11}`$ cm-2 with a peak mobility $`\mu =5.5\times 10^6`$ cm2/Vs. The $`\nu =5/2`$ state shows a strong minimum in diagonal resistance and a developing Hall plateau at magnetic fields ($`B`$) as high as 12.6T. The strength of the energy gap varies smoothly with $`B`$-field. We interpret these observations as strong evidence for a spin-polarized ground state at $`\nu =5/2`$.

\]

Ever since its discovery in 1987  the origin of the even-denominator fractional quantum Hall effect (FQHE) state at Landau level filling factor $`\nu =5/2`$ remains mysterious. Results of recent ultra-low temperature experiments  leave no doubt that this even-denominator state is a true FQHE with vanishing resistivity and Hall plateau formation. Such an even-denominator state does not fit the normal odd-denorminator rule of the FQHE and requires other or additional electron correlations. An early, so-called Haldane and Rezayi (HR) model  arrives at a spin-unpolarized trial wave function.

Eisenstein et al.  tested the spin-polarization of the $`\nu =5/2`$ state by tilted magnetic field experiments  in a traditional, fixed-density sample. While the orbital motion of the electrons and hence their correlation energy ($`E_c`$) is subject only to the perpendicular component of the magnetic ($`B`$) field, the Zeeman energy ($`E_z`$) depends on the total $`B`$-field. Varying angle and $`B`$-field, the specimen can be kept in the $`\nu =5/2`$ state while the Zeeman energy is raised. Such a procedure should leave a spin-polarized state intact, but should be detrimental to a spin-unpolarized state, once the Zeeman energy cost surmounts the gain in correlation energy. In the experiment the strength of the $`\nu =5/2`$ state decreased quickly upon tilting and the state disappeared totally at $`\theta 50^{}`$. This was taken as evidence of a spin-singlet state at $`\nu =5/2`$.

In recent years, with the advent of the composite fermion (CF) model , there has been a renewed interest in the $`\nu =5/2`$ state . Moore and Read (MR) proposed a ground state of $`p`$-wave paired CF’s . Unlike the HR state, the MR state is spin-polarized. It is now being argued that the earlier disappearance of the $`\nu =5/2`$ state under tilt may be the result of the compression of the wave function due to the in-plane component of the $`B`$-field. This reduction of the $`z`$-extend of the wave function affects electron correlation and is the cause for the gap collapse, rather than the suspected increase in Zeeman energy. Indeed, recently two tilted field experiments  showed that the added in-plane magnetic field not only destroys the FQHE at $`\nu =5/2`$ but also induces an electronic transport anisotropy. Theoretical modeling  suggests that this is due to a phase transition from the MR pairing state to an anisotropic state and unrelated to any spin-effect. Hence, the spin-polarization of the $`\nu =5/2`$ FQHE remains unresolved and there is presently little experimental input into the debate over the nature of the even-denominator FQHE state.

Here we pursue the spin-polarization of the $`\nu =5/2`$ state in analogy to the tilted field experiments by investigating the competition between $`E_c`$ and $`E_z`$. Rather than tuning their ratio in a fixed density specimen by tilt, we keep the $`B`$-field perpendicular to the 2DES and employ a variable density specimen. Since for a fixed filling factor, such as $`\nu =5/2`$, $`E_cn^{1/2}`$, whereas $`E_zn`$, increasing electron density modifies the ratio of $`E_z`$ to $`E_c`$. This approach is equivalent to tilting the sample, but it cannot cause a tilted-field induced phase transition . In our density-dependent experiment on the $`\nu =5/2`$ state we observe a strong minimum in diagonal resistance and a developing Hall plateau even at very high electron densities, equivalent to $`B`$-fields as high as 12.6T. Furthermore, the strength of the energy gap varies slowly and smoothly with $`B`$-field. We interpret these observations as evidence for a spin-polarized ground state at $`\nu =5/2`$.

To perform these experiments we fabricate a HIGFET (heterojunction insulated gate field-effect transistor), in which the 2DES density can be tuned from $`n=0`$ to $`7.6\times 10^{11}`$ cm-2 with a peak mobility $`\mu =5.5\times 10^6`$ cm2/Vs. The left inset of Fig. 1 shows the profile of such a HIGFET , with a heavily $`n^+`$-doped GaAs cap layer as the gate. In contrast to modulation-doped GaAs/AlGaAs heterostructures , there is no delta-doping layer in a HIGFET, which allows to reach high mobilities at very high densities. We have investigated as many as 15 HIGFETs. In all specimens we observed qualitatively the same results. In this letter, we present experimental data from the specimen, which allowed for the highest gate-voltages and, hence, highest electron densities.

The right inset of Fig. 1 shows the strictly linear dependence of density $`n=1.542\times V_g0.082`$ on $`V_g`$ where $`n`$ is in units of 1011 cm-2 and $`V_g`$ in units of volt. The mobility $`\mu `$, shown in Fig. 1, rises monotonically to a value as high as $`\mu 5.5\times 10^6`$ cm2/Vs at $`n=5.5\times 10^{11}`$ cm-2. Beyond this density electrons are starting to populate the second electrical subband . This opens an additional scattering channels and leads to a reduction in mobility.

We have measured $`R_{xx}`$ and $`R_{xy}`$ at $`\nu =5/2`$ over the density range from $`n=3.0\times 10^{11}`$ to $`7.6\times 10^{11}`$ cm-2. The high-density limit is set by the breakdown voltage of the HIGFET and for densities much lower than $`n=3.0\times 10^{11}`$ cm-2 the $`\nu =5/2`$ state becomes weak. In Fig. 2(a)-2(c), we show the diagonal resistance $`R_{xx}`$ in the interval $`3>\nu >2`$ at three selected densities, $`n=(3.0,5.3,7.6)\times 10^{11}`$ cm-2. In all cases there exists a strong minimum in $`R_{xx}`$ at $`\nu =5/2`$. This indicates the existence of a $`\nu =5/2`$ FQHE state over the whole range of density from $`n=3.0`$ to $`n=7.6\times 10^{11}`$ cm-2. The $`\nu =5/2`$ FQHE state at 12.6T in Fig. 2(c) represents, to our knowledge, by far the highest $`B`$-field at which this state has ever been observed.

Since we are limited to standard dilution-refrigerator temperatures of $`T>30mK`$, $`R_{xx}`$ never vanishes in our HIGFET data, in contrast to our previous ultra-low temperature measurements . For this reason we cannot determine the energy gap, $`\mathrm{\Delta }`$, from activation energy measurement. In order to quantify our results we employ an earlier method  to attribute a strength, $`S`$, to the $`R_{xx}`$ minimum. $`S`$ is defined as the ratio of the depth of the minimum to the average around $`\nu =5/2`$, $`S=R_{5/2}/R_{ave}`$ (see inset Fig. 3(a)). $`S`$ varies exponentially in temperature, $`S`$ exp($`\mathrm{\Delta }^{quasi}/2k_BT`$) and defines a quasi-energy gap, $`\mathrm{\Delta }^{quasi}`$. Since $`S`$ measures a quantity very similar to $`R_{xx}`$ and is proportional to it, at least at higher temperatures, we can take $`\mathrm{\Delta }^{quasi}`$ to be very similar to $`\mathrm{\Delta }`$.

Fig. 3(a) shows the strength, $`S`$, of the $`\nu =5/2`$ state versus inverse temperature (1/$`T`$) for three selected densities. For $`n<3.0\times 10^{11}`$ cm-2 the $`\nu =5/2`$ becomes too weak to be quantifiable in terms of $`S`$. On this semilog plot the data follow a linear relationship and $`\mathrm{\Delta }^{quasi}`$ is readily deduced. Fig. 3(b) shows the value of $`\mathrm{\Delta }^{quasi}`$ at $`\nu =5/2`$ filling (solid dots) versus $`B`$-field. $`\mathrm{\Delta }^{quasi}`$ hovers around 150mK and there is little variation over this magnetic field range. Our previous, ultra-low temperature measurement  on a fixed density ($`2.3\times 10^{11}`$ cm-2) specimen yielded an activation energy of $`\mathrm{\Delta }`$ = 110mK. A comparison between both data indicates that the value of the $`\mathrm{\Delta }^{quasi}`$ is a good approximation to $`\mathrm{\Delta }`$. In Fig. 3(b) $`\mathrm{\Delta }^{quasi}`$ increases slightly from the smallest $`B`$-field to about 9T, whereupon it decays somewhat to higher fields. The dependence is smooth and no indication for any sharp transition is apparent.

From the data of Fig. 3(b) we conclude that it is highly unlikely that the $`\nu =5/2`$ state is spin-unpolarized. At 12.6T the Zeeman energy amounts to $`E_z`$ 3K, which is more than a factor of 15 larger than the biggest energy gap ever measured for any $`\nu =5/2`$ state, which typically vary from 100-200mK. The overwhelming strength of the Zeeman energy should have overcome any correlation-induced spin-reversal and should have destroyed such a spin-unpolarized FQHE state. The smooth field-dependence of $`\mathrm{\Delta }^{quasi}`$ in Fig. 3(b) also indicates that no phase transition occurs in the $`\nu =5/2`$ state over this field range. Therefore the polarization of the spin system remains unchanged, implying a spin-polarized state for the entire range of Fig. 3(b). For comparison we plot in the same figure the $`\mathrm{\Delta }^{quasi}`$ of the $`\nu =8/5`$ state, measured in the same specimen. The well-studied transition from a spin-unpolarized to a spin-polarized state  expresses itself as a strongly reduced energy gap at $``$7.5T in contrast to the smooth variation of the gap at $`\nu =5/2`$.

The above qualitative argument can be quantified with the help of Fig. 4. Here we plot the dependence of $`E_c=\alpha e^2/ϵl_B`$ and $`E_z=g\mu _BB`$ versus $`B`$-field with $`l_B=(\mathrm{}c/eB)^{1/2}`$ being the magnetic length, $`ϵ=12.9`$ the dielectric constant of GaAs, and $`\mu _B`$ the Bohr magneton. We use $`g=0.44`$ for the $`g`$-factor in GaAs  and $`\alpha =0.02`$ from a numerical calculation of Morf . For a spin-unpolarized state one would take the energy gap to be $`\mathrm{\Delta }=E_cE_z`$ which is represented by the difference between $`E_c`$ and $`E_z`$ in Fig. 4. $`E_z`$ equals $`E_c`$ at a critical $`B`$-field of $`B_{crit}11`$T. For higher fields the Zeeman energy exceeds the correlation energy and, in a simple model, one would expect a spin-unpolarized state to be no longer stable. With the usual uncertainties in the theory and the simple nature of our model 12.6T is insufficiently far from this $`B_{crit}`$ 11T to totally rule out a spin-unpolarized ground state on this basis. However, theoretical energy gaps always exceed experimentally measured energy gaps. The reason is probably an inherent level broadening, $`\mathrm{\Gamma }`$, due to disorder, which subtracts from the theoretical gap. From various experiments  this broadening is believed to be $``$0.5-1K for the standard high-mobility specimens and roughly magnetic-field independent within a given specimen. Such a broadening shrinks the range of observable gaps as indicated in Fig. 4. We have chosen $`\mathrm{\Gamma }=0.6`$K, which best reflects the value of $`\mathrm{\Delta }`$ 0.1-0.2K measured in most samples in the range of 3-5T. This reduces $`B_{crit}`$ to less than 7T, considerably lower than our highest value of 12.6T. Alternatively, in order to reflect the experimental values of $`\mathrm{\Delta }`$ 0.1-0.2K measured in most samples in the range of 3-5T the theoretical value for $`\alpha `$ may be an overestimate and consequently $`\mathrm{\Gamma }`$ may be smaller than $`\mathrm{\Gamma }`$ 0.6K. Yet, inspection of Fig. 4 shows that any such variation would only reduce the upper critical field. A theoretical value larger than $`\alpha =0.02`$ is very unlikely, contrary to the trends in the experimental data and contrary to the results of other few particle numerical calculations . On the basis of these comparisons one can be very confident, that the $`\nu =5/2`$ FQHE state is spin-polarized over the entire density range over which it has been experimentally observed.

There remain some observations that need to be addressed. The gap of a spin-polarized system should follow the correlation energy $`E_cn^{1/2}B^{1/2}`$. The data of Fig. 3(b) initially show such a gradual increase but eventually turn around at $`B`$ 9T. The origin of this behavior is unclear, but is probably related to the increasing confinement of the electrons against the interface due to the rising gate voltage. Electron scattering from residual interface roughness increases, leading to an increase in $`\mathrm{\Gamma }`$ and a deceasing energy gap for higher gate voltages and higher electron densities. (Note that such an $`n`$-dependent $`\mathrm{\Gamma }`$ would further reduce $`B_{crit}`$ in Fig. 4). The decrease is not related to the population of a second subband as seen in the $`B=0`$ mobility data in Fig. 1, although both occur at similar densities. In the range of the $`\nu =5/2`$ state Landau quantization is operative and a simple, $`B=0`$ intersubband scattering model is no longer valid. Furthermore, a comparison between Landau level splitting ($``$22meV) and electrical subband splitting ($``$19meV) at 12.6T indicates that the $`\nu =5/2`$ state resides in the second Landau level of the lowest electrical subband, as in all previously studied, fixed-density specimens. The absence of any discontinuity of $`\mathrm{\Delta }^{quasi}`$ in Fig. 3(b) indicates that this condition persists over the whole density range investigated. In any case, the existence of a $`\nu =5/2`$ FQHE state, residing in the lowest Landau level of the second electrical subband is, a priori, very unlikely. Such a state is equivalent to $`\nu =1/2`$, where CF’s form a fermi sea and not a CF-paired FQHE state.

In summary, in a variable-density specimen we observe a $`\nu =5/2`$ FQHE state at a magnetic field as high as 12.6T. The high Zeeman energy at this $`B`$-field and a detailed comparison with a universal model calculation exclude a spin-unpolarized ground state. The quasi-energy gap of the $`\nu =5/2`$ FQHE varies smoothly over the whole density range from $`n=3.0`$ to $`7.6\times 10^{11}`$ cm-2, which excludes a transition between spin-polarizations. From these findings we conclude that the even-denominator $`\nu =5/2`$ FQHE is spin-polarized, consistent with a Moore-Read paired CF state.

We would like to thank E. Palm and T. Murphy for experimental assistance. A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by NSF Cooperative Agreement No. DMR-9527035 and by the State of Florida. D.C.T. and W.P. are supported by the AFOSR, the DOE, and the NSF.


text:
# 1 Schematic behaviour of the Kertész line for the Ising model

Cluster Percolation and

Pseudocritical Behaviour in Spin Models

Santo Fortunato and Helmut Satz

Fakultät für Physik, Universität Bielefeld

D-33501 Bielefeld, Germany

Abstract:

The critical behaviour of many spin models can be equivalently formulated as percolation of specific site-bond clusters. In the presence of an external magnetic field, such clusters remain well-defined and lead to a percolation transition, even though the system no longer shows thermal critical behaviour. We investigate the 2-dimensional Ising model and the 3-dimensional $`O(2)`$ model by means of Monte Carlo simulations. We find for small fields that the line of percolation critical points has the same functional form as the line of thermal pseudocritical points.

Percolation theory  provides in many cases an elegant interpretation of the mechanism of second order phase transitions: the existence of an infinite spanning cluster represents the new order of the microscopic constituents of the system due to spontaneous symmetry breaking.

In the Ising model, a rigorous correspondence between thermal critical behaviour and percolation was established by Coniglio and Klein ; the clusters are constructed by linking like-sign spins through a temperature-dependent bond probability $`p=1\mathrm{exp}(2J/kT)`$, where $`J`$ is the Ising coupling and $`T`$ the temperature. This result has been recently extended to a wide class of models, from continuous spin Ising-like models  to $`O(n)`$ models . In all cases the system was considered in the absence of any external field. The reason for this is clear: by introducing an external field $`H`$, we explicitly break the symmetry of the Hamiltonian of the system, thus eliminating the thermal critical behaviour of the model. None of the thermodynamic potentials exhibits discontinuities of any kind, since the partition function is analytical for $`H0`$.

On the other hand, the Coniglio-Klein clusters can be built as well when $`H0`$. Because of the field, the system has a non-vanishing magnetization $`m`$ parallel to the direction of $`H`$ for any finite value of the temperature $`T`$. For $`T\mathrm{}`$, $`m\mathrm{\hspace{0.17em}0}`$, and for $`T=0`$, $`m=1`$. This suggests that for a fixed value of $`H`$, the clusters will start to form an infinite network at some temperature $`T_p(H)`$. Varying the field $`H`$, one thus obtains a curve $`T_p(H)`$ in the $`TH`$ plane, the Kertész line .

The Kertész line specifies the usual percolation threshold, just as in the case $`H=0`$, and the percolation variables exibit the usual singularities, leading to a set of critical exponents. In particular, the percolation strength, which is the relative size of the percolation cluster compared to the total volume of the system, remains the order parameter of the percolation transition.

The existence of such a line of genuine critical behaviour also for $`H0`$ could provide a criterium to define different phases and the transition between them in a more general sense . For this purpose, it is important to find out whether there is any relationship between the geometrical singularities at the Kertész line and thermal properties of the system along this line; apparently, rather little is known about such connections -.

In the present work we want to address this aspect in more detail. In particular, we shall investigate the relation between the Kertész line and the so-called ”pseudocritical” or ”crossover” line, defined by the values of the temperature at which the magnetic susceptibility peaks as function of the field. In general, the two lines cannot coincide, since when $`H\mathrm{}`$, the Kertész line leads to a finite $`T`$, while for the pseudocritical line this results in $`T\mathrm{}`$; we return to this point later on. We shall here study two spin models, the Ising model in two dimensions and the $`O(2)`$ model in three dimensions, and calculate in both cases the functional form of the Kertész line and the pseudocritical line in the limit of a small external field.

The Hamiltonian of the models at study can be written in the following general form:

$$\beta =\frac{J}{kT}\underset{i,j}{}𝐬_𝐢𝐬_𝐣\frac{H}{kT}\underset{i}{}𝐬_𝐇𝐬_𝐢,$$
(1)
where $`𝐬_𝐢,𝐬_𝐣`$ are two-dimensional unit vectors for $`O(2)`$ and simple scalars ($`\pm 1`$) for the Ising model; similarly, $`𝐬_𝐇`$ is a unit vector (or a scalar) in the direction of the external field. We now first consider the percolation problem.

In Fig. 1 we show schematically the Kertész line for the Ising model. For a vanishing external field, $`H=0`$, the Coniglio-Klein cluster definition with

$$p_{i,j}p_b=1\mathrm{exp}\{2J/kT\}$$
(2)
as bond weight between two adjacent spins $`i,j`$ assures that we recover the usual thermal threshold of the Ising model. For $`H\mathrm{}`$, all lattice spins will be aligned with the field at any temperature $`T`$. However, the bonds between adjacent spins will be occupied according to the bond weight $`p_b`$, so that now the site-bond problem turns into pure bond percolation. The percolation transition will then take place for that value of the temperature $`T_b`$ for which the probability $`p_{ij}`$ equals the critical density $`p_b(d)`$ of random bond percolation in $`d`$ dimensions,

$$p_b(d)=1\mathrm{exp}\{2J/kT_b\},$$
(3)
i.e., for

$$T_b=\frac{2J}{k\mathrm{log}[1p_b(d)]}$$
(4)
These arguments can simply be repeated for the Kertész lines corresponding to $`O(n)`$ spin models.

As mentioned above, each point of the Kertész line is a standard percolation point. At $`H=0`$, the critical exponents of the percolation variables coincide with the thermal critical exponents of the magnetization transition (Ising, $`O(n)`$). However, for any $`H0`$, they are predicted to switch into the (different) universality class of random percolation in the same dimension. We will verify this for each of the percolation points we shall determine.

Turning now to the pseudocritical thermal behaviour of the models considered, we define the reduced temperature $`t`$ and external field $`h`$ measured relative to the spin-spin coupling strength

$$t=\frac{TT_c}{T_c},h=\frac{H}{J}.$$
(5)
For spin systems, the renormalization group approach leads to the following expression of the magnetization $`m`$ when $`t`$$`h\mathrm{\hspace{0.17em}1}`$

$$m=h^{1/\delta }f_G(t/h^{1/\beta \delta }),$$
(6)
where $`f_G`$ is a scaling function and $`\beta `$, $`\delta `$ are the critical exponents for the magnetization transition at $`h=0`$. We can rewrite Eq. (6) in the form

$$m=t^\beta (t/h^{1/\beta \delta })^\beta f_G(t/h^{1/\beta \delta })=t^\beta F_G(t/h^{1/\beta \delta })$$
(7)
with a new scaling function $`F_G`$; the entire dependence on the field $`h`$ is now put into the function $`F_G`$.

The susceptibility $`\chi `$ is given by

$$\chi \left(\frac{m}{H}\right)_Tt^\beta F_G^{}(x)th^{1/\beta \delta 1},$$
(8)
where $`F_G^{}(x)`$ is the derivative of the function $`F_G`$ with respect to its argument $`x=t/h^{1/\beta \delta }`$. This time we move the dependence on the reduced temperature $`t`$ into a new function $``$ of the argument $`x`$ and obtain

$$\chi h^{1/\delta 1}(t/h^{1/\beta \delta })$$
(9)
To find the equation of the pseudocritical line, we have to determine the temperature $`t_\chi (h)`$ at which the susceptibility $`\chi `$ peaks for a given value $`h`$ of the field. Hence the derivative of $`\chi `$ with respect to the reduced temperature $`t`$ must vanish at $`t=t_\chi `$,

$$\frac{\chi }{t}|_{t=t_\chi }=\mathrm{\hspace{0.17em}0}h^{1/\delta 1}^{}(x)|_{t=t_\chi }(1/h^{1/\beta \delta })=h^{1/\delta 11/\beta \delta }^{}(x)|_{t=t_\chi }=\mathrm{\hspace{0.17em}0};$$
(10)
here $`^{}(x)|_{t=t_\chi }`$ is the derivative of the function $``$ with respect to its argument $`x`$ calculated at $`t=t_\chi `$. The derivative $`\chi /t|_{t=t_\chi }`$ can be zero only if $`^{}(x)|_{t=t_\chi }=0`$. That will occur at some value $`x_\chi `$ of the argument $`x=t_\chi /h^{1/\beta \delta }`$, which provides the relation between $`t_\chi `$ and $`h`$

$$t_\chi /h^{1/\beta \delta }=x_\chi t_\chi h^{1/\beta \delta }.$$
(11)
The procedure is obviously independent of the value of the field $`h`$, so that as long as $`h`$ is small, the pseudocritical line is described by a simple power law.

We stress that, for $`h\mathrm{}`$, the temperature of the susceptibility peak diverges, $`t_\chi \mathrm{}`$, whereas we have seen that the Kertész line has a finite endpoint, given by Eq. (3). Hence the two curves, which start from the same point but tend to two different limits, must certainly differ for sufficiently large values of $`h`$. We want to investigate here, however, what happens in the vicinity of the thermal critical point of the model, i.e. for small fields ($`h\mathrm{\hspace{0.17em}0}`$).

In order to study the percolation transition we need to redefine each given spin configuration as a cluster configuration, by grouping all spins of the lattice in clusters. We do this with the help of the algorithm of Hoshen and Kopelman , using free boundary conditions for the cluster identification and assuming that a cluster percolates if it connects each pair of opposite sides (faces in 3d) of the lattice.

After measuring the size of all clusters, i.e., the number of sites in each cluster, we evaluate the following two quantities:

* The average cluster size  $`S`$,

$$S\frac{_sn_ss^2}{_sn_ss},$$
(12)
where $`n_s`$ is the number of clusters of size $`s`$; the sums exclude the percolating cluster.
* The percolation strength  $`P`$,

$$P=\frac{\text{size of the percolating cluster}}{\text{total number of lattice sites}},$$
(13)
which is the order parameter of the percolation transition.

In the infinite volume limit, the percolation variables are described by simple power laws

$$P(T_pT)^{\beta _p}(T<T_p),S|TT_p|^{\gamma _p}$$
(14)
sufficiently close to the critical threshold $`T_p`$.

In order to determine the critical point of the percolation transition, we exploit the properties of a further very useful variable. For finite lattice size, there can be spanning clusters even for temperatures $`T>T_p`$, and there are spin configurations at temperatures $`T<T_p`$ below the critical threshold without such clusters. The probability of finding a spanning cluster on a finite lattice of linear dimension $`L`$ at a temperature $`T`$ is a well-defined function, the percolation cumulant $`\mathrm{\Pi }(T,L)`$ . For $`TT_p`$ and large $`L`$, it behaves as

$$\mathrm{\Pi }(T,L)=\mathrm{\Phi }\left[\left(\frac{TT_p}{T_p}\right)L^{1/\nu _p}\right].$$
(15)
where $`\nu _p`$ is the critical exponent governing the divergence of the percolation correlation length. The function $`\mathrm{\Pi }(T,L)`$ is not a genuine percolation variable, since it has a non-trivial meaning only on finite lattices. On an infinite lattice it reduces to a step function: it is zero for $`T>T_p`$ and unity for $`T<T_p`$. Nevertheless, its particular features make it a powerful tool to extract information about the critical properties of percolation. In particular, at $`T=T_p`$, $`\mathrm{\Pi }(T_p,L)=\mathrm{\Phi }(0)`$ for any value of $`L`$. That means that if we calculate the percolation cumulant as a function of $`T`$ for different lattice sizes, all curves $`\mathrm{\Pi }(T,L)`$ will cross at the critical temperature $`T_p`$. Moreover, if we plot the different $`\mathrm{\Pi }(T,L)`$ as a function of the variable $`X=[(TT_p)/T_p]L^{1/\nu _p}`$, the resulting functions $`\mathrm{\Phi }(X)`$ must coincide for all lattice sizes. To define the scaling function $`\mathrm{\Phi }(X)`$, one also has to know the value of the exponent $`\nu _p`$. We can therefore determine $`\nu _p`$ by looking for the best scaling form of the percolation cumulant curves.

Let us now first consider the results obtained for the 2D Ising model. Our aim is to determine the functional form $`t_p(h)`$ describing the Kerteśz line. In order to compare it to the pseudocritical line, we have to perform simulations at small values of the reduced field $`h`$. We consider the following four values: $`h`$= 0.002, 0.001, 0.0005 and 0.00025. The Monte Carlo algorithm used in the simulations is the Wolff cluster update extended to the case of a non-vanishing external field . For each iteration we have measured the energy $`ϵ`$ of the configuration, its magnetization $`m`$, and the percolation variables $`S`$ and $`P`$. To determine the percolation critical point with a good accuracy, we work with rather large lattices: $`300^2`$, $`500^2`$, $`700^2`$, $`1000^2`$. This is helpful also for another reason. Since the field is quite weak, the scaling behaviour of both thermal and percolation variables could be perturbed by the “near-by” case of the system without field, which could affect the values of the critical exponents. Such perturbation becomes smaller the larger the size of the lattice.

Fig. 2 shows the percolation cumulants as a funtion of $`\beta =J/kT`$ for different lattice sizes, at the smallest field we have studied, $`h=0.00025`$. The curves cross clearly at the same point ($`\beta _p=0.43933`$). To check the universality class of the percolation exponents in this case, we rescale the percolation cumulants as discussed in Section 3. Using $`J/kT_p=0.43933`$, we consider two options for the critical exponent $`\nu _p`$, the 2-dimensional random percolation value $`\nu _{RP}=4/3`$ and the 2D Ising exponent $`\nu _{IM}=1`$. The results are shown in Figs. 4 and 4; they clearly indicate that the critical behaviour at the Kertész line is determined by the random percolation exponent $`\nu _{RP}`$, as expected.

Using a standard finite-size scaling analysis of the average cluster size $`S`$ near the critical point, we obtain for the ratio of the exponents $`\gamma _p/\nu _p=1.80(2)`$, in accord with the value predicted from the random percolation exponents (129/72=1.7916) and not with that from the Ising model (7/4=1.75). Repeating the percolation cumulant study for the other three values of $`h`$, we find the critical percolation temperatures listed in Table 1 and plotted in Fig. 5. In the Figure, the reduced percolation temperature $`t_p=(T_pT_c)/T_c`$ is shown as function of the external field; here $`T_c`$ is the critical temperature of the Ising model without field ($`J/kT_c=(1/2)\mathrm{log}(\sqrt{2}+1)=0.44068679`$). On the logarithmic scale of the plot the data points fall remarkably well on a straight line. We thus conclude that, for the Kertész line of the 2D Ising model,

$$t_p=a_ph^\kappa ,\mathrm{for}h\mathrm{\hspace{0.17em}0},$$
(16)
where $`a_p`$ is a constant of proportionality. From the slope in Fig. 5 we obtain $`\kappa =\mathrm{\hspace{0.17em}0.534}(3)`$. On the other hand, the exponent of the thermal pseudocritical line is $`1/\beta \delta `$ (see Eq. (11)). For the 2D Ising model, $`\beta =1/8`$, $`\delta =15`$, so that $`1/\beta \delta =8/15=0.533`$. The agreement between this value and the exponent found for the Kertész line is excellent.

We thus ask ourselves whether the Kertész line and the pseudocritical line overlap in the range of $`h`$ values we have chosen. To answer this, we determine some points of the pseudocritical line and compare them to the corresponding points of the Kertész line. Fig. 6 shows the susceptibility as a function of $`\beta `$ for $`h=0.001`$, calculated on the two largest lattices we have studied. The curves are obtained by interpolating the data points using the density of state method . The Figure clearly indicates that there is no noticeable change of the peak between $`700^2`$ to $`1000^2`$, so that we are effectively at the infinite volume limit. The position of the peak in the Figure, however, clearly differs from the corresponding percolation point at $`h=0.001`$, indicated by the dotted line.

To complete the comparison, we determine the positions of the susceptibility peaks for two other values of the field, $`h=0.0005`$ and $`h=0.002`$. As above, we find clear discrepancies with the corresponding points of the Kertész line; the results are included in Fig. 5, where it is also seen that the power law of Eq. (11)

$$t_\chi =a_\chi h^{1/\beta \delta },a_\chi =\mathrm{const}.$$
(17)
continues to be valid with great precision. We conclude that for the 2D Ising model the functional dependence of the reduced temperature $`t`$ on the field $`h`$ for $`h0`$ is for the Kertész line and for the pseudocritical line described by the same function, determined by the thermal critical exponent $`1/\beta \delta `$; but the two curves do not coincide, since $`a_pa_\chi `$, as is evident in Fig. 5.

We now turn to the corresponding study of the 3D $`O(2)`$ model. Here we have determined four points of the Kertész line, for external field $`h`$ values 0.001, 0.00050, 0.00025, and 0.000125, using again the cluster update  as Monte Carlo algorithm, on lattices of size $`30^3`$, $`40^3`$, $`50^3`$, $`60^3`$, $`70^3`$, $`80^3`$, $`100^3`$. The clusters we considered are those used in the Wolff algorithm with no external field (see ).

The positions of the percolation points, determined as before by the crossing of the percolation cumulant curves, are listed in Table 2. As for the Ising model, we have checked that the universality class of the exponents of each percolation transition is that of random percolation, here in three dimensions. We show these points in Fig. 7, again plotting the $`h`$-dependence of the reduced percolation temperature $`t_p=(T_pT_c)/T_c`$. The critical temperature $`T_c`$ of the magnetization transition for the 3D $`O(2)`$ model without field was found to be $`J/kT_c=0.454165`$ in .

Again, the straight line in the log-log plot clearly indicates the power law behaviour of the Kertész line. The slope of the straight line is $`0.601(3)`$, in good agreement with the exponent of the pseudocritical line $`1/\beta \delta =0.5995(6)`$, obtained by using the $`O(2)`$ thermal exponents values determined in .

To compare the two lines directly, we make use of the results of a recent numerical determination of the $`O(2)`$ pseudocritical line . The results are included in Fig. 7. The two lines are parallel in the log-log plot, since the exponents of $`h`$ agree, but they do not coincide.

We have shown that the Kertész lines of the 2D Ising model and the 3D $`O(2)`$ model are in the limit of small fields described by the same functions as the corresponding thermal pseudocritical lines, i.e., a power law with exponent $`1/\beta \delta `$. This feature, which is likely to be general, suggests a relationship between the geometrical percolation singularities which still occur for $`H0`$ and the remaining non-singular thermal properties of the system. In both models we have found that the Kertész line does not coincide with the pseudocritical line even in the scaling region we have explored. The difference in the behaviour of the two lines could well be a consequence of using the conventional $`h=0`$ bond weight (2) also for $`h0`$; introducing a dependence on the field $`h`$ into the Coniglio-Klein factor might well make the two lines coincide. Simple expressions of the form $`1\mathrm{exp}[2\beta (1+h)]`$ would lead to a Kertész line which goes to infinity for $`h\mathrm{}`$, so that an eventual interpretation of the pseudocritical line as a curve of percolation points is not yet excluded. Further work on this is in progress.

Acknowledgements

It is a pleasure to thank J. Engels for helpful discussions. We would also like to thank the TMR network ERBFMRX-CT-970122 and the DFG Forschergruppe Ka 1198/4-1 for financial support.


text:
# 1 Introduction

## 1 Introduction

Munitions containing depleted uranium (DU) have been used by NATO and US forces during the war operations in Iraq (1991), Bosnia (1994), Kosovo and Serbia (1999). Recently some information on 112 sites targeted by DU weapons in Kosovo has been supplied by NATO to the United Nation Environmental Program Balkans Task Force (UNEP BTF); on November 2000 measurements to detect contamination have been undertaken by a UNEP team in 11 among the 112 sites.

Aim of this paper is outlining some aspects of the exposure of people living in an area contaminated by DU, on the basis of official available information and of simulations, looking for main pathways of average and critical exposure.

Individuation of pathways of high exposure could allow to advice to population; average dose assessment, together with measures of DU concentration in soil, will make delimitation of areas to be cleaned up possible.

## 2 Military use of depleted uranium

The Gulf war against Iraq in 1991 was the first one known where DU rounds have been used in large quantity (approximately 300 tonnes) . The consequences on the health of the Iraqi population and of the US veterans are still under study. DU exposure at the moment is not considered the most probable cause of the Gulf War Syndrome experienced by hundreds thousand veterans ; on the other hand, the effects of the DU left over the Iraqi territory are difficult to show, due to the large number of toxic substances dispersed in the environment during the war and the deterioration of the sanitary situation caused by the embargo to which the country is submitted from 1991 (cited work in , app.3).

Reports on potential effects on human health and environment from the use of DU have appeared during the last years: studies on risk assessment for the Jefferson Proving Ground, a US facility for testing DU munitions, have been performed ; the risk for population for the Kosovo conflict and for the Gulf war has been also considered .

DU can be obtained as by-product in the enrichment process of natural uranium for the production of nuclear fuel and for military applications; as the ore extracted natural uranium, DU is associated to a reduced chain of radioactive isotopes, formed by $`{}_{}{}^{238}U`$ and $`{}_{}{}^{235}U`$ decay products having shorter decay times: $`{}_{}{}^{234}Th`$ (24 days), $`{}_{}{}^{234m}Pa`$ (1.17 min) and $`{}_{}{}^{234}Pa`$ (6.7 hours), 231 Th(25.5 hours). DU can also be obtained by the reprocessing of nuclear power plant spent fuel, and so traces of transuranic elements and $`{}_{}{}^{236}U`$ can be present. According to official information, DU used by the U.S. Department of Defence contains approximately 0.2$`\%`$ of $`{}_{}{}^{235}U`$ and traces of $`{}_{}{}^{234}U`$, $`{}_{}{}^{236}U`$. Following the indications in , we will assume the uranium isotopic composition of DU given in Table 1.

DU specific activity is in part due to uranium isotopes (14.9 $`Bq/mg`$, 36%), and for the residual part to beta emitting short-life decay products (64%); among the transuranic elements official information is available only for $`{}_{}{}^{239}Pu`$ (2.4 104 years), whose content is estimated in 11 ppb . DU specific activity is not substantially affected by the declared amount of traces elements.

Metallic uranium has a high density (19 $`g/cm^3`$), is pyrophoric and cheaper than tungsten, and so has been attractive for U.S. Army for the production of armor piercing ammunition since 1960s. Tungsten alloys have been preferred up 1973, when a DU alloy with 0.75$`\%`$ of titanium (U-3/4Ti) was adopted for ammunition made by a thin cylinder in DU alloy encased with lighter material. Systems of DU weapons are owned or under development in different countries (Saudi Arabia, France, United Kingdom, Israel, Pakistan, Russia, Thailand and Turkey) .

Use of DU ammunition causes exposure of people soon and after, because DU is dispersed as aerosol when the projectile strikes a hard target and then falls out on a limited area . Contamination of all environmental matrices takes place and health effects on people living nearby must be taken in account, both for toxicological damage and for radiological risk. Among different isotopes present in DU as declared, $`{}_{}{}^{238}U`$,$`{}_{}{}^{234}U`$ and $`{}_{}{}^{235}U`$ are of concern in risk assessment. For chemical hazard, kidney is identified as the target organ, whatever the path of assumption . Due to prevalent short-range emitted radiation, the risk associated with exposure to ionizing radiation mainly derives from ingestion and inhalation of radioactive material; external irradiation from soil is less relevant.

## 3 Dispersion of DU in the environment and exposure of the population

DU contained in projectiles, spread out as aerosol in air after striking the target, falls out producing environmental and food chain contamination. Possible occurring of chemical hazard and entity of radiation dose must be assessed for people living in the area, taking into account both average and critical group exposure.

DU concentration in the soil is the starting point; while waiting for measurements of contamination in Iraq, Bosnia, Kosovo and Serbia, we present computed radiation doses and associated concentrations for different contaminated soil thickness, as soil mixing will extend the initial superficial deposition to underlying layers in not undisturbed areas. Available soil measured DU concentrations in contaminated sites that we are aware of, are the following:

* at Jefferson Proving Ground area an average $`U`$ concentration of 318 $`Bq/kg`$ was reported ; more recently a lower and an upper bound of the concentration ranging from 592 $`Bq/kg`$ to 13690 $`Bq/kg`$ was also measured ;
* among the areas where the US personnel lived in the Gulf region (outside Iraq) the highest DU concentration (433 $`Bq/kg`$) was measured in the Iraqi Tank Yard (the area where captured Iraqi equipment is stored in Kuwait) ;
* in some sample analyzed by the RFY scientists a specific activity of $`{}_{}{}^{238}U`$ up to $`2.3510^5Bq/kg`$ was detected .

Following the hypothesis assumed in the BTF report, we have assumed as a reference value a DU contamination of $`1000Bq/kg`$ of soil over an area of $`A=10000m^2`$, in the hypothesis of 10 $`kg`$ of DU entirely dispersed in the impact as aerosol of uranium oxides, contaminating 1 $`cm`$ of soil. With the composition given in Table 1 initial activities per $`kg`$ of soil for $`{}_{}{}^{238}U`$, $`{}_{}{}^{235}U`$, $`{}_{}{}^{234}U`$ and $`{}_{}{}^{236}U`$ are respectively 830 $`Bq`$, 11 $`Bq`$, 154 $`Bq`$ and 5 $`Bq`$.

Average effective dose is conservatively assessed using the residential farmer scenario. The following pathways are considered: external irradiation from soil, inhalation from resuspended dust, ingestion of contaminated soil and water, ingestion of plants and animal products grown in site and ingestion of fish grown in a pond contaminated by groundwater. Different pathways are considered for plant contamination due to first root uptake (water independent) and due to secondary root uptake from use of contaminated water (water dependent). Radon inhalation is excluded. RESRAD 5.91  code is used, all parameters default except for the ones given in Table 2. Estimates of dose to individuals and population for risk in contaminated sites have been performed by EPA employing primarily the code RESRAD (for related work see ).

RESRAD default libraries values have been corrected to give effective dose  rather than equivalent effective dose : due to the algorithm used by RESRAD, anyway, values for external irradiation $`E_G`$ in Tables 4 and 5 have been impossible to modify, and are approximate by 10% maximum defect.

In Tables 3, 4, 5 and 6 we show average annual effective doses and corresponding DU concentrations in water and vegetables for three different soil thickness, respectively $`1cm`$, $`10cm`$ and $`20cm`$. The following quantities are given at different times, from the first year to about two hundred years after maximum dose, for main pathways: the total dose ($`E_{tot}`$), the dose from external irradiation from the soil ($`E_G`$), from inhalation of contaminated dust ($`E_I`$), from consumption of edible plants (water independent $`E_P`$, water dependent $`E_P^w`$) and of water ($`E_{H_2O}`$).

The dependence of $`t_{max}`$ and $`E_{max}`$ on some hydrogeological parameters, mainly affecting the water dependent pathways, is shown in Table 7 and 8. The maximum value of the dose is not much affected by most of the parameters considered in Table 8 except $`K_d`$. This parameter is defined as the ratio of the mass of solute species observed in the solids per unit of dry mass of the soil to the solute concentration in the liquids. A wide range has been observed for uranium $`K_d`$ values . For largest value of $`K_d`$ the DU is retained in surface and does not reach at least within the first 1000 years the watertable. A measurement of the local value of this parameter is therefore necessary to reduce the uncertainty on the dose assessment.

Strong dependence of maximum inhalation dose has been found, as expected, on the dust loading parameter, as shown in Table 9.

As already outlined, presented doses and concentrations have been obtained from an average value of soil contamination, in order to assess the average exposure of population. Whatever the average value considered, anyway, highly inhomogeneous soil concentrations must be expected in the contaminated area, both for sparse aerosol deposition and for oxidation of DU fragments: concentrations up to 12% in weight have been reported . In order to assess the dose to critical population group, this must be taken in account, especially if inhalation of soil was the critical pathway: inhalation of 0.1 $`g`$ of soil with maximum reported DU contamination, equal to 12 $`mg`$ DU, corresponds to 1.44 $`mSv`$; ingestion of 1 $`g`$ of soil, equal to 120 $`mg`$ DU, corresponds to 0.08 $`mSv`$.

A scenario, in which permanence in dusting air and ingestion of soil are possible, is the one for children playing with soil. From the presented dose assessment and considerations children playing with soil may be identified as the critical population group, with inhalation and/or ingestion of contaminated soil as critical pathway. Evidently, average and critical doses are somehow competitive, because the higher fraction of DU is dispersed as aerosol, the lower part of it can rest in soil as fragment, being presence of fragments the main cause of hot spots in soil contamination.

It must be outlined that the amount of DU considered in the simulation corresponds to 37 A-10 /GAU-8 ammunitions. According to the available information, a much larger number of projectiles has been fired on each site (between 50 and 2320, average 300) and up to now unknown is the extension of targeted sites. Both for average and critical exposure, anyway, more realistic dose assessment will be possible only when measured contamination data will be known, scaling the values in the tables for the appropriate factor.

Increment of inhalation dose attributable to $`{}_{}{}^{239}Pu`$ presence in DU is officially estimated in 14$`\%`$ : with 11 ppb of $`{}_{}{}^{239}Pu`$ in DU RESRAD gives a maximum dose increment of 0.6%.

## 4 Normative and recommendations framework

Before discussing compliance of average assessed doses and exposure with international standards set to prevent from toxicological damage and limit ionizing radiation risk, we shortly line out an aspect relative to radioprotection system, maybe useful even in wider considerations on risk.

Due to accepted linear-no-threshold model for effects produced by ionizing radiation, justification of a practise has to be the first one posed, that is if the population exposure from military use of DU is justified or not. Comparison between dose estimates in such a scenario and dose limits and dose constraints stated by regulations is anyway useful, for a quantitative perception of risk. In order to assess the need for remediation in contaminated areas, once again the question of justification has to be considered; specific reference levels, linked to the avertable annual dose, have to be defined by national authorities. ”Generic reference levels … should be used with great caution” and their use ”should not prevent protective actions from being taken to reduce … dominant components \[of existing annual dose\]” . We next report a comment to assessed doses, comparing them with radiological and toxicological reference values, in order not to hold the question narrowed to exceeding of dose limits.

Values of annual dose in Tables 3, 4 and 5 show the same temporal shape, with an initial prevalent dose from irradiation by soil and a maximum from ingestion of contaminated drinking water occurring after about five hundreds years, when contamination reaches the acquifer serving the population. Maximum dose, progressively increasing as inventary of DU increases, is always lower than annual population limit (1 $`mSv/y`$), starts to be comparable with EPA cleanup limit criterion (150 $`\mu Sv/y`$, ) for 20 $`cm`$ depth. Exceeding of dose constraint of $`0.1mSv/y`$ indicated in  for longlived isotopes may not be excluded. This in general happen only after long times, due to the low mobility of the uranium oxides (the mean transit times for insoluble uranium in the top 10 cm of soil range from 7.4 to 15.4 years with an average of 13.4 years ; soluble forms have a mean transit times of one month). At the maximum dose time concentration of DU in the water reaches the provisional value of WHO guideline for drinkable water (0.05 $`Bq/l`$ ) already for 1 $`cm`$ depth. The concentration of DU in leafy vegetables at time of maximum dose ranges from 2 to 4 $`Bq/kg`$; no derived limit is defined for consumption of dietary parts.

Inhalation of highly contaminated soil may leave to exceeding of annual dose limit, with possible occurring of toxicological damage: maximum allowed concentration in air for workplaces stated by NRC, 45 $`\mu g/m^3`$ for soluble and 200 $`\mu g/m^3`$ for insoluble uranium forms, would be exceeded if dust loading was more than 1700 $`\mu g/m^3`$, a high but not extreme value. Less important seems ingestion of contaminated soil, due to the lower value of dose conversion factor with respect to the inhalation one. Anyway, ingestion of 1 $`g`$ maximum contaminated soil would result in 120 $`mg`$ DU ingestion, when maximum daily ingestion of uranium, due to toxicological effects, was stated in 150 $`mg`$ by italian legislation till year 2000.

## 5 Conclusions

DU contained in projectiles, spread out in air after striking the target, falls out producing environmental and food chain contamination. Possible occurring of chemical hazard and entity of radiation dose must be assessed for different kind of exposure of people living in the area, taking into account both average and critical group exposure. While waiting for measurements of contamination in Iraq, Bosnia, Kosovo and Serbia, we have computed radiation doses and concentrations for different contaminated soil thickness, as soil mixing will extend the initial superficial deposition to underlying layers in not undisturbed areas.

In order to assess the average exposure of population, doses and concentrations have been obtained from an average value of soil contamination. For the individuation of the critical group inhomogeneous soil concentration has been considered. The presented dose assessment suggests a short term exposure due to inhalation and/or ingestion of contaminated soil and a long term exposure due to ingestion of contaminated water and food; the propagation of the superficial contamination to the watertable critically depends on various hydrogeological parameters to be evaluated on the site.

In sites targeted by DU munitions special cautions have to be adopted to reduce unaware exposures and cleanup must be planned on the basis of the measured concentrations.


text:
# Evolutionary route to diploidy and sex

## I Introduction

The evolution of sex through Darwinian selection in spite of the seeming odds with respect to asexual reproduction through cloning, is a still unresolved problem , which has attracted a great deal of recent activity . We would like to present here results from a very simple model based on unicellular organisms where sex seems to have first emerged, without any complicating factors such as the cost of male-female differentiation, parental care, mating preferences, maternity periods and so on, which we assume only to have arisen in higher organisms. These factors may be responsible for the differentiation, modification or even at times reversal (as in meiotic parthenogenesis ,p.234) of the patterns that may have become accidentally established in the remote evolutionary past. Unicellular organisms are not only simpler to model, but this choice also imposes more stringent conditions on the establishment of a sexual population, since it severly limits the number of viable offspring that a pair of parent cells can have. By confining ourselves to “worst case scenarios” we hope to be able to obtain plausible lower bounds to the feasibility of sex as a reproductive mechanism.

Sexual reproduction as we know it today in higher plants and animals can be regarded as a haploid-diploid cycle with a highly abbreviated haploid phase, where the (haploid) gametes typically do not perform mitosis.  These haploids are not viable unless two of them fuse and once more make a diploid cell, which then multiplies by mitosis and eventually forms the mature individual.

The haploid diploid cycle (HDC) where, under unfavorable conditions the population becomes diploid, is found in many species , and was the motivation behind the hypothesis of Jan, Stauffer and Moseley  (JSM) who proposed that diploidy and sex may have first arisen as a way to escape death, when a simple, unicellular individual is threatened by too many deleterious mutations. The fitness is taken to be a step function of the number of mutations, also setting the threshold for conversion to diploidy and sex. Here, the number of deleterious mutations can be read as the distance from the ideal case, called the wildtype, which can be altered as a result of a change in the environment, so that the individual is not as well adapted as before the change.

In previous work  we showed that the JSM hypothesis indeed leads to a steady source of diploid sexuals, given a population of haploids which multiply by mitosis. This leads to a steady state distribution of coexisting haploids and diploids in a constant population. The premises adopted in  were very restrictive. Two sexual “parent” cells were allowed to have only one offspring, after which they died. In case one considers greater number of offspring, by allowing a greater number of viable gametes to be formed, as is indeed possible in many unicellular organisms , one finds that the diploid, sexual population completely takes over.

In the same paper, we have also considered the situation where conversion to sex occurred with a constant probability $`\sigma `$ (which we varied between 0 and 1), independently of the fitness of the individual. This strategy also gave rise to a sexual population making up a small but macroscopic fraction of the total, so that the JSM threshold mechanism for conversion to sex did not prove to be necessary, although it was more successful. On the other hand the JSM hypothesis has the additional attraction of providing a mechanism which could trigger the fusion of two haploid cells to form a diploid: it is known ,p.149, that extensive genetic damage can lead to gene repair via genetic transmission between two haploid cells; the fusion of two haploids could be seen as an extreme form of such behaviour. An alternative means of forming a diploid from a haploid cell is via endomitosis ,p.230, as in the first step in meiotic parthenogenesis. Endomitosis is the process whereby the genetic material in a cell is duplicated without subsequent cell division as in normal mitosis. Again, it is not implausible that this process originated as a result of grave genetic damage which precluded the successful completion of mitosis. , p.149.

In the present paper we show that, if the diploid cells, once formed, are also allowed to multiply by mitosis, as indeed they do in a haploid-diploid cycle, the whole population is taken over by diploid cells, who perform facultative sex if they are once more threatened by extinction due to too many deleterious mutations. Moreover we show that even an episodic conversion to sex, involving as few as only two individuals who survive to mate, leads to a steady state made up solely of sexual types.

In the next section we briefly describe our algorithms and report the simulation results. In section 3, we provide a discussion and some pointers for future research.

## II Bit-string model for the Conversion to Sex - Algorithms and simulation results

The type of model we consider here is the same as in Ref. . Below and in Section A, we recall the definitions. In Section B and C we introduce new rules to test the autonomous viability of sexual populations.

Each haploid one-celled organism consists, for our purposes, of a 15-bit string of “0”s and “1”s, representing the genetic code in a 16-bit computational word.  We use the bit defining the “sign”, to specify whether the individual is asexual (+) or sexual (-). A mutation consists of flipping a randomly chosen bit except the sign bit, with a constant probability, $`\mathrm{\Gamma }`$, for each individual per generation. Since the genetic difference between individuals of the same species is typically less than $`10\%`$ ,p.52, this rather short string for the genetic code may be considered as a coarse grained model for the complete genome of the individual, which we divide up into different zones, retaining a “0” where there are no mutations, and flipping the bit to “1” if there are one or more mutations in this zone. The wild type is a string of all “0”s. Therefore, at each locus, a “1” corresponds to a deleterious mutation (which we will call “mutation,” for short, where this is not liable to lead to any confusion.) Diploid organisms have two bit strings, which are allowed to be different. The number of deleterious mutations $`m`$ is simply the number of “1”s for a haploid individual. For a diploid, the number of “expressed” deleterious mutations is the number of loci at which both homologous alleles are set to “1,” i.e., we assume that deleterious mutations are recessive.

The total population is fixed at $`N=10^3`$, and we have chosen $`\mathrm{\Gamma }=1/N`$. This corresponds to a mutation rate per allele per generation of $`6\times 10^5`$, which is comparable to the typical mutation rates encountered in eukaryotes ,p.64.

The probability of survival (or fitness) as a function of $`m`$ is given by a Fermi-like distribution , $`P(m)`$,

$$P(m)=\frac{1}{\mathrm{exp}[\beta (m\mu )]+1}.$$
(1)
For large $`\beta `$ (or “low temperatures,” in the language of statistical mechanics), $`P(m)`$ behaves like a step function. Individuals with $`m>\mu `$ die, those with $`m<\mu `$ survive, and those with $`m=\mu `$ survive with a probability of $`1/2`$. In the simulations we confined ourselves to low temperatures ($`\beta =10`$). We chose $`\mu =4`$ which just allows us enough variability without leading to totally unrealistic mutational loads.

We start with a set of $`N`$ asexual (haploid) wildtypes. In each generation about $`\mathrm{\Gamma }N`$ individuals suffer mutations; they are killed off or retained according to the fitness function (1), and the population restored to $`N`$ by once cloning as many randomly chosen survivors as necessary. The population of haploid asexuals settles down to a minimally stable  steady state distribution (see Table I) as shown in Fig. 1, independently of the value of $`\mathrm{\Gamma }`$, for $`\mathrm{\Gamma }1/N`$

For comparison, we have also performed simulations on a diploid population reproducing asexually, according to the same rules as stated in the previous paragraph. We have found that the diploid population reaches a steady state with an $`m`$ distribution peaked at $`m=2`$, rather than $`m=3`$ as found for haploids. (See Fig. 1 and Table I).

### II.1 Conversion to sex

Let us briefly summarise our algorithm for conversion to sex, which we also used in our previous study (see Model A, Ref. ).

We choose the steady state of the haploid population as our initial state from which to start the conversion to sex, as that would be a most likely “natural state” encountered at this stage. Once the haploid population reaches a steady state, we allow those individuals that are threatened by extinction to convert to diploidy and sex, by implementing the following rule at each generation: If an asexual individual with $`m=4`$ has survived, it is converted to an active sexual, by deterministically and irreversibly switching its sign bit to (-). It performs endomitosis and becomes diploid. If there are already sexual, diploid organisms in the population, they will also be made active if $`m`$, the number of their expressed mutations, is $`4`$ (otherwise they do not participate in the reproduction cycle; hence the conversion to sex may be termed non-hereditary). Finally, all the active sexual organisms pair randomly and engage in sexual reproduction, where they each contribute one gamete (formed via one step meiosis) towards a single diploid sexual offspring. If we denote the genotypes of parents as $`\{Aa\}`$ and $`\{Bb\}`$ respectively, then the genotype of the offspring is either $`\{AB\}`$, $`\{Ab\}`$, $`\{aB\}`$ or $`\{ab\}`$. No crossover occurs during this one step meiosis. If there is only one active sexual at a certain time step then it must wait subsequent generations until it either finds a partner or it dies. We keep the population constant  by cloning randomly selected haploid individuals to make up the deficit.

(The consequences of the conversion to hereditary and obligatory sex can be found in Ref. .)

### II.2 Mitotic Diploid Sexuals Win over Haploids

In the previous subsection we described a scenario where too big a mutation load meant conversion to diploidy and sex, if the individual survived. It should be realised that in this scheme the haploids, who multiply by cloning, provide a steady source for the diploid sexuals, whose numbers are halved every time they mate. Now we want to test whether the diploid, sexual population can survive autonomously, if the sexuals are also allowed to perform mitosis themselves.

Once a steady state with coexisting haploid and diploid populations is reached via the algorithm described above, we switch off the conversion of haploid organisms to diploidy and sex. The diploid individuals mate when they face extinction due to too many mutations. We now keep the total population constant by making up for the deficit in the population at each step by cloning randomly selected individuals, regardless of whether they are sexual or asexual. In colonies undergoing a haploid-diploid cycle, it is quite frequently the case  that the diploid phase of the cycle also involves multiplication by mitosis.

The result is that the diploid, sexual individuals completely win over the population. The haploids which now can not compete with the sexuals, become extinct. Here we find that the haploid phase of the HDC becomes abbreviated to the point where haploids appear only as gametes which do not perform mitosis. This is exactly the situation in highly evolved sexual organisms.

### II.3 Episodic conversion to diploidy and sex leads to sexual steady state

One can also think of a scenario in which conversion to sex takes place accidentally over a short period of time after the asexual population settles down to a steady state. In fact this is probably the most realistic situation, given the random nature of the mutations.

The way we actually implemented this in the computer code was by deterministically switching the sign bits of the first two asexual individuals to survive with $`m=4`$, and then turning off the possibility of further conversion. These then form two sexual individuals by endomitosis, and if the first survives long enough so that it can mate with the second, will give rise to one sexual offspring.

The rest of the rules are as explained before; in each generation we clone randomly chosen individuals to make up the deficit in the population, regardless of whether they are asexuals or sexuals. We allow the diploid individuals to mate when they face extinction in the course of their lives.

Suprisingly, the sexuals capture the population in 95 percent of the performed runs. (In the rest, the single diploid, which is still at the threshold of extinction with an $`m`$ value of 4, may not survive until a partner arrives.)

### II.4 Relative Fitness

The steady state $`m`$ distributions for the sexual populations treated in Section II.B and C are much better than the haploid distribution. However, they are identical with the distribution for asexual diploids.(see Table I and Fig. 1) This means that as soon as diploidy is achieved, it is so successful in screening the effects of deleterious mutations that once diploid, the organism practices sex very infrequently, making these three distributions identical.

The diploid $`m`$-distribution in Fig. 1 should be contrasted with the results of Ref.  (Model A, non-hereditary and hereditary; also see Subsection A, above) where cloning of the sexuals is not allowed, but there is a steady influx of new sexuals from the haploids (see Table I). In the non-hereditary model, the $`m`$ distribution has shifted nearer the wildtype than in the purely diploid case (although the peak is again at $`m=2`$); in the hereditary model, the peak has shifted to $`m=1`$. Thus, the frequency with which sex is practiced has a salutary effect on the distribution of expressed deleterious mutations in the population.

On the other hand it should be noted that, with the Fermi-like fitness function (Eq.(1)) for large $`\beta `$, to which we have confined ourselves, the survival probability is $``$ constant for $`m<\mu `$, and therefore does not discriminate between steady state populations that differ solely in the shape of their $`m`$ distributions for $`m<\mu `$. A better measure might be the deathrate, namely, the average number of individuals which are eliminated at each step, under the criterion given in Eq.(1).

Our findings, which are rather revealing, are summarized in Table I. Since the fluctuations were very large, we have only reported order of magnitude results for the death rates. Pure diploids have the lowest death rate, an order of magnitude smaller than pure haploids. Mitotic (including the “episodic” case) sexuals are very close to diploids. For non-hereditary sexuals the rates are comparable to asexuals, but hereditary sexuals have a death rate that is two orders of magnitude larger than asexuals. We infer from this that the success of hereditary sexuals in lowering their mutation load is due to the much greater rate at which they can select-out highly mutated genes through death.

Thus it is seen that the clear advantage of sex over asex (haploid or diploid) in our model can only become manifest with a sufficiently high frequency of sexual reproduction, which is driven by the mutation rate $`\mathrm{\Gamma }`$, and moreover in situations which can discriminate between the different adaptabilities (smaller typical $`m`$) of the different modes of reproduction, such as time varying environments.  This is similar to the finding of Pekalski , who has considered environments (wildtypes) that vary over time, and found that the benefits of sexual reproduction are enhanced by higher birthrates. For lower birthrates, meiotic parthenogenesis, which can be compared to diploid asex, has a slight advantage over sex, a difference from our results which can be ascribed to the exponential fitness function he uses.

## III Discussion

This paper is a culmination of a series of studies where we have considered very stringent rules for the survival and propagation of diploid, sexual individuals in competition with haploids. Here we have finally been able to show that a pair of simple, unicellular organisms who have accidentally converted to diploidy, and which subsequently engage in sexual reproduction, begetting one sexual offspring, can give rise to a population of sexual types which totally take over a finite population, provided they are also allowed to multiply by mitosis, on an equal footing with the haploids in the population. They engage in sex when the going gets tough, that is, when the number of their expressed deleterious mutations exceeds a certain number. This success seems to vindicate the hypothesis of JSM , that sex could have been a mechanism of last resort when simple organisms were faced with extinction.

It is interesting to pose the question of what happens if we do not allow the diploids to clone themselves, but, on the other hand, allow the gametes, under special conditions, to enter a haploid phase where they multiply by mitosis. Eventually these haploids will be allowed once more to fuse and give rise to diploids. Within the present scheme, since the haploids convert to diploidy only at the threshold, $`m=4`$, the gametes of the diploid individuals have $`m4`$, i.e., they are not viable. Therefore the extension of the haploid phase calls for a modification of the rules in such a way as to allow the gametes to survive nevertheless, e.g., by forming “spores.”

One way the parameters of the present model could be modified is to take the threshold for survival to be different (greater than) the threshold for the conversion to (or practice of) sex. Then, if a haploid gamete were to be considered as a member of the haploid population, and subjected to the same rules for selection and reproduction, it would seem as if one would get haploid phases of arbitrary duration. However, upon closer inspection we see that this is not true: Since the gamete comes from a diploid organism which engages in sex due to an excess mutational load, it already has deleterious mutations in excess of (or at least equal to) the sex threshold. (For the diploid organism to have $`m`$ expressed mutations, both of its genetic strings must have at least that many). So, according to the rules which we have adopted, the abbreviated haploid cycle is automatically selected, since within this scheme the gamete would immediately seek a mate. The lengthening of the haplophase would call for further evolution, e.g., such mutations as would lead to the enhancement of the fitness of the haploid gametes with $`m`$ in excess of the sexual threshold (say by the formation of spores, to survive in an environment which has become disfavorable ). Since this will provide a mechanism for the diploids once more to make up for the reduction in their numbers due to the single-offspring mating rule, we expect this HDC to win out over the haploid, asexual population.

To understand the seemingly contradictory of results from similar models investigating the phenomenon of the selection of sex as the dominant mode of reproduction, it is of interest to compare our model with that of Redfield  and Cui et al., . The former employs a more complex system in which there is sexual differentiation between males and females, with the additional feature of a much greater mutation rate for males than for females, eventually making meiotic parthenogenesis safer than sex. Redfield has used three different “selection functions” (survival probabilities), falling off exponentially (no “epistasis,” or correlation between the effects of successive mutations) quadratically (some positive epistasis) or as a step function (extreme positive epistasis) with $`m`$. The fact that the step-function (“truncation”) survival probability alone leads to a much better fitness resulting from sex, tells us that the form of the survival probability which is adopted determines the outcome very strongly. The step function favors sex because it does not punish small deviations from the wildtype, thereby allowing genetic diversity, while it severely penalizes large deviations, helping eliminate highly mutated individuals.

This insight enables us to understand the results of Cui et al. who have considered a diploid population reproducing either asexually or sexually, depending upon a fixed probability. There is a constant mutation rate. The survival probability is chosen to be of the “independent mutation,” i.e., exponential type. These authors find that the diploid population accumulates such a large number of deleterious mutations (which however do not get expressed), that sexual reproduction, i.e., the random pairing of gametes from different parents, results in a disasterous reduction in the fitness of the offspring. Moreover, sexually reproducing diploid populations are susceptible to invasion by asexually reproducing ones, while the converse is not true, i.e., asexually reproducing diploid populations are not invaded by sexually reproducing ones. These results change drastically in favor of sex when they introduce cell senescence, which means a cell can clone itself only a certain number of times before it stops dividing and dies, unless it engages in sexual reproduction, which resets the senescence clock.

We would like to argue that in the model of Cui et al., introducing cell senescence is very much like turning on a “truncation like” survival function, which, as in our model, also controls the switch to sexual reproduction. In the presence of a constant mutation rate, the number of times a cell has cloned is another way to measure its (mean) mutation load, which, in an asexual haploid population of course gets directly expressed, whereas for diploid individuals hides behind the dominant unmutated alleles. Introducing a cell senescence threshold, beyond which the individual is either killed off or has to engage in sexual reproduction, is therefore equivalent, in an average way, to the JSM criterion, applied to a diploid population. Seen in this way, it is very gratifying that the Cui et al. model in fact corraborates our findings regarding the greater fitness of diploid populations that frequently engage in sex, given a step-function like survival probability.

Very recently, two variants of the Redfield model, which assign “harsher selection” to males (rather than relatively higher mutation rates) have been considered . These studies both find that the mutational load on the whole population decreases as a result, and that the relative fitness of sexual females is increased by more than enough to compensate for the twofold cost of sex in anisogamous populations, at least for extremely high values of the average mutation rate. Agarwal finds that this can happen even in the absence of synergystic epistasis, in his model. Wilke et. al  have meanwhile found that high mutation rates select for low replication rates and flat regions of the fitness surface.

It should be noted that the Fermi-type fitness function (1) adopted here extrapolates, for finite temperatures (smaller values of $`\beta `$) between highly synergistic (step function) and independent (exponential) survival probabilities as a function of $`m`$. A recent study by Peck and Waxman indicates  that competition for limited resources can lead to synergy between successive mutations, leading to truncation, or step function-like survival probabilities, which they also find favors sexual populations.

In this paper we have presented numerical results for a model involving unicellular organisms, which may shed light on how sex and diploidy emerged and established a foothold in the protozoon world. Many different mechanisms have so far been proposed via which sex may prove advantageous or otherwise, in more highly evolved organisms. We believe that care should be taken while proposing any single mechanism, such as an adaptation to resist infestation by parasites , for the preferance of sex by organisms that range from the unicellular eukaryotes to trees or human beings . More complex organisms may have elaborated much more complex survival mechanisms and behaviorial patterns, which stabilize or destabilize already evolved traits.

Acknowledgements

We thank Aslıhan Tolun and Canan Tamerler Behar for many useful discussions. One of us (A.E.) gratefully acknowledges partial support from the Turkish Academy of Science.


text:
# Comment on “Density-matrix renormalization-group method for excited states”

## Abstract

In a recent paper (Phys. Rev. B 59, 9699 (1999)), Chandross and Hicks claim to present a new density matrix renormalisation group (DMRG) method for dealing with excited states of quantum lattice models. The proposed improvement to the DMRG—the inclusion of excited state wave functions in addition to the ground state in the density matrix when calculating excitations—is in fact standard pratice, is clearly stated in White’s original papers, and has been used repeatedly by many groups to study excited states. The authors apply the method to the extended, dimerised Hubbard model for conjugated polymers. The criteria for determining whether states are bound or not are assessed. The authors claim that their results show that the optically important “$`1B_u`$” state is bound (excitonic), in contrast to a previous study. However, the discussion is qualitative, and the authors arrive at conclusions on the basis of results for one lattice size only. We show that when Chandross and Hicks’ criterion is developed into a quantitative definition of particle-hole separation, with the finite-size dependence analysed, the implication is that the $`1B_u`$ state is unbound, in keeping with the conclusions of a previous study.

In a recent paper , Chandross and Hicks claim to present a new density matrix renormalisation group (DMRG) method  for dealing with excited states of quantum lattice models. They apply the method to the dimerised, extended Hubbard model for conjugated polymers. They claim that a previous study  of this model is flawed because it uses a “conventional” DMRG method which does not handle excitations correctly. The improvement that they suggest is to form a density matrix not only from the ground state, but from all the states being targeted in the calculation. This is in fact standard practice in DMRG calculations of excited states and the structure of the density matrix required to target excited states is given in White’s original papers on the method . It has been used by many authors to target excitations in a variety of quantum lattice models (see, e.g., ) and was certainly used in  when various excitation energies and correlation functions were calculated for the extended Hubbard model. The comparisons presented in Fig. 1 and Fig. 2 of , between the “conventional” DMRG and Chandross and Hick’s “improvement” are therefore of limited value, as, to the best of our knowledge, all DMRG studies of excited states to date have incorporated the targeted excitations into the density matrix . Unfortunately, a slightly different value for the Coulomb $`V`$ is used in  so a direct comparison with the results (e.g., for energies) tabulated in  is not possible. We have run a DMRG program which uses the algorithm used in  for targeting excited states with the parameters $`U=3t`$, $`V=t`$, $`\delta =0.1`$, used in , and found good agreement for the energies and correlation functions with the results plotted in Fig. 1(a) and Fig. 2(a) of . For instance, we plot the $`1B_u`$ and $`mA_g`$  energies as functions of the lattice size $`N`$ in Fig. 1. The results compare very well with Fig. 1(a) of .

In  Chandross and Hicks also examine criteria for deciding whether a particular excitation is bound (excitonic) or not. They claim that the average particle-hole separation, defined in  in terms of the density-density correlation function, is “too approximate” a quantity to determine whether a state is bound or not. They argue that by inspecting the centered correlation function as a function of distance (together with the profile of doubly occupied sites along the chain), for one particular lattice size ($`N=36`$ sites), one can see that the $`1B_u`$ and $`mA_g`$ states are “different” in that the $`1B_u`$ ($`mA_g`$) has its strongest particle-hole correlations at short (long) distances. However, Chandross and Hicks do not present an alternative quantitative definition of the particle-hole separation, based on this observation. In , on the other hand, it is argued that a definition of particle-hole binding must take into account the way in which correlations scale with lattice size $`N`$. In  it is argued that this scaling is different for bound and unbound excitations, and that the scaling of the average particle-hole separation with $`N`$ is but one manifestation of this.

Suppose we wish to take the average double occupancy of the $`1B_u`$ state (relative to the ground state) along the chain $`(\widehat{n}_i1)^2_{1B_u}(\widehat{n}_i1)^2_{1A_g}`$ as an example (Fig. 2(a) in ). In Fig. 2 we plot this quantity for various lattice sizes $`N`$. We see that, although the concentration of doubly occupied sites is greatest in the middle of the chain, the distribution spreads out as $`N`$ is increased. The area under these curves rapidly converges to a non-zero value $`(0.538)`$ as $`N\mathrm{}`$. This shows that the number of pairs of particles and holes in the $`1B_u`$, relative to the number in the ground state, approaches a constant. Our results could indicate that particle-hole pairs separate as $`N`$ is increased and are hence unbound, or they may simply indicate dispersion of a bound exciton in the $`1B_u`$.

To address this, we again consider the averaged, centered, odd-site correlation function $`C_{1B_u}^N(i)`$, (again relative to the ground state value), defined in  and plotted for $`N=36`$ in the inset to Fig. 2(a) in . In Fig. 3 we plot this quantity for a number of values of $`N`$. We see that, although the correlations are generally strongest at short distances, they become increasingly spread out, and hence the particle-hole pair becomes increasingly separated, as $`N`$ is increased. Indeed, if one utilises $`|C_{1B_u}^N(j)|`$ to define a probability distribution for the particle-hole separation, as in , then one finds that the average particle-hole separation grows linearly with $`N`$, as shown in Fig. 4. We note that any use of the density-density correlation function to describe particle-hole separation and the nature of exciton binding of excited states in the extended Hubbard model is merely plausible rather than rigorous , but Chandross and Hicks  do not offer an alternative quantitative definition of the particle-hole separation to the ones provided in .

To summarise, Chandross and Hicks claim that, because the $`1B_u`$ and $`mA_g`$ have their greatest particle-hole correlations at short and long distances respectively (on the $`N=36`$ lattice), the $`1B_u`$ is bound and the $`mA_g`$ is unbound. We would argue that it indicates that the particle-hole potential is more strongly attractive for the $`1B_u`$ state than for the $`mA_g`$. However, from the plausible, quantitative definition of the particle-hole separation given above, it would appear that the attraction between the particle and hole in the $`1B_u`$ state is not strong enough to bind them, and their separation increases throughout the range of lattice sizes studied.

Finally, we consider the structure of the density matrix when targeting excitations such as the $`mA_g`$ and $`nB_u`$. Chandross and Hicks argue that only four states—the $`1A_g`$ (ground state), the $`1B_u`$, the $`mA_g`$ and the $`nB_u`$—need be included in the density matrix. Our examinations of the dipole moments between the $`A_g`$ states and the $`1B_u`$ indicate that this approach is probably reasonable for the $`mA_g`$ which is well defined. That is, there is a reasonably abrupt jump in the magnitude of the dipole moment $`1B_u|\widehat{\mu }|jA_g`$ at $`j=m`$. As shown in , this coincides with jump in the ionicity (the average number of doubly occupied sites) and in the particle-hole separation. However, the $`nB_u`$ state is less well defined in that there can be a number of $`B_u`$ excitations that have a strong dipole moment with the $`mA_g`$. This can be seen in Table I where we list the dipole moments $`jB_u|\widehat{\mu }|mA_g`$ for $`N=6`$, 10, 14 and 18, for the first five $`B_u`$ states. Note that in no case is the $`nB_u`$ state clearly defined, though there is a general trend whereby the $`2B_u`$ increases its relative dipole strength with the $`mA_g`$ at the expense of the $`4B_u`$. Our contention here, as proposed in , is that, at least in terms of dipole moments or the density-density correlation function, the $`1B_u`$ state is the threshold of unbound states in the $`B_u`$ sector and the “$`nB_u`$” is not well defined for this model.

Calculations were performed at the New South Wales Center for Parallel Computing. This work was supported by the Australian Research Council.