text:
# Modeling the actinides with disordered local moments

## I Introduction

The elemental actinide metals, Pu in particular, exhibit several unique features. They are among the most complex elements in nature, with a rich set of allotropes, of which several have complex low symmetry crystal structures . Pu, for example, is the only element with seven condensed matter phases at zero pressure of which one (the $`\delta `$-phase) demonstrates negative thermal expansion. The understanding of these anomalous properties is a serious challenge. Several of the main features of the actinides can be understood in terms of the progressive filling of the 5f shell. In the light actinides, Th to Np, the 5f electrons are itinerant and participate in the bonding, whereas in the heavier actinides, Am to Cf, the 5f electrons are localized and exhibit behavior more similar to the lanthanides. Plutonium takes a particular place between these extremes. The low temperature $`\alpha `$-phase has been shown to be well described with itinerant bonding 5f electrons whereas the high temperatures phases with their increased volumes suggests a localized or partly localized configuration. This behavior indicates a Mott-like transition of the f-electrons similar to the $`\alpha `$-$`\gamma `$ transition in cerium .

There have been several first principles approaches to calculating the properties of the actinides . Although the light actinides are well described within density functional theory, the local density approximation (LDA) fails in general to describe the 5f localization which occurs in actinide compounds, elemental Pu at high temperature, and in the actinides past Pu. A number of techniques have therefore been developed to include correlation effects beyond LDA in order to describe the partially localized 5f electronic structure. These include LDA+U, orbital polarization, and the Mixed Level Model . Of particular relevance to the work presented here are recent attempts to combine LDA and Dynamical Mean-Field Theory (DMFT). DMFT has been shown to be a very useful approximation at the consideration of strongly correlated electron systems (for review, see ). LDA combined with the DMFT may provide a first principles technique to study correlated electron materials, and there have been several attempts to apply implementations of LDA+DMFT ; for a review, see . Recently this LDA+DMFT approach has been applied to the Pu problem . The LDA+DMFT approach gives an opportunity to describe correlation effects on the electronic structure and properties of $`d`$\- and $`f`$\- electron systems. However, the technique is cumbersome and it is not completely ab initio because of the problem with the choice of $`U`$ (see, e.g. ). In addition, implementations of DMFT are computationally intensive which makes the calculation of complex structures difficult. Sometimes, because of the complexity of the calculations, uncontrollable approximations are made, e.g., using a single $`U`$ parameter instead of a complete interaction matrix .

It is commonly accepted that accounting for Hubbard correlation effects is of crucial importance for $`f`$-electron systems. On the other hand, some of these effects can be taken into account, in an approximate way, in the framework of more traditional density functional techniques. This is the approach taken in the present study where we have modeled the actinides by means of Disordered Local Moments (DLM)  within the Local Spin Density Approximation and the Coherent Potential Approximation  (LSDA+CPA). The purpose of the present work is to investigate and possibly demonstrate that correlation effects beyond the standard local density approximation can be simulated, at least partially, by means of a parameter free first principles DLM approach, and to see if this picture gives an adequate description of the actinides. The DLM picture, even if it is insufficient for a complete description, might lead to some insights for the understanding of the electronic structure of the actinides.

Consider, first, the status of the DLM approach in the many-body lattice models like Hubbard or s-f exchange (“Kondo lattice”) model. The key point is an equivalence between a many-body interacting system with Coulomb on-site interactions and a one-electron system in fluctuating charge and spin fields. This equivalence, which can be proven by the Hubbard-Stratonovich transformation, is a base of spin-fluctuation theories of itinerant-electron magnetism . In a complete theory, the charge and spin fields are dynamically fluctuating both in space and time. However, a “static approximation”, where we neglect the dynamics of the fluctuations, captures an important part of the correlations while greatly simplifying the formalism, and may be sufficient for many problems of interest. In this case the correlated system is described in terms of a DLM alloy and a CPA for this alloy becomes equivalent to the ”Hubbard III” approximation  for the original many-body problem (see ). The Hubbard III approximation, and therefore also the corresponding DLM description, gives qualitatively the correct picture of the electronic structure both in atomic and broad-band limits and can describe the metal-insulator (Mott) transition for half-filled bands when $`U`$ is of the order of the bandwidth . The most important difference with the more elaborate DMFT picture  is an absence of the ”Kondo peak” near the Fermi level at the metallic side of the transition. However, this Kondo peak has a small spectral weight and, generally speaking, it is not obvious that it is of crucial importance for the description of total energy and related characteristics. It is shown in Ref. that the Hubbard III approximation can be rigorously justified, in the s-f exchange model, in the limits of large space dimensionality and classical spins ($`S\mathrm{}`$) whereas for DMFT only the first limit is essential. Note, that the Kondo resonance is essential also for a correct description of the properties of electron states in doped Mott insulators (large -$`U`$ limit near the half filling). On the other hand, for half filling the Hartree-Fock-like descriptions of the band-split states are, in general, quite adequate . Physically, this is the consequence of a freedom in the description of occupied bands in terms of Bloch functions, as well as Wannier functions or any other orthogonal linear combinations. The spatial fluctuations of the exchange on-site field can lead to a splitting of the energy spectra (provided the fluctuations are larger than the bandwidth) which corresponds, in the Hubbard model, to the Hubbard band splitting. The self-energy in the DLM picture is energy dependent and has an imaginary part describing the damping of the electrons on spin fluctuations. This distinguishes the DLM approach not only from magnetically ordered LSDA calculations but also from LDA+U, SIC, and the Hartree-Fock approximation. In this sense, the DLM-CPA approach can be considered as a particular case of the “LDA++” approach with a local, energy dependent, complex self-energy. However, in contrast to schemes taking into account Hubbard correlations, it is completely ab initio. A shortcoming of this approximation is an incorrect description of electron damping near the Fermi level; i.e., the absence of the Kondo resonance, which is a consequence of the quantum character of spin, as well as problems describing localization in systems that do not have half filled (or completely filled) electron shells. However, this shortcoming may be of minor importance for the description of electron energy spectrum at large energy scales as well as for the calculations of total energy and related characteristics such as the equilibrium volume and elastic moduli. Note also that this finite damping, because of electron scattering by spin fluctuations, is a physically correct picture for high enough temperatures.

Comparing the many-body lattice models with the density functional approach one can face with the well-known “Hubbard $`U`$ vs. Stoner $`I`$” problem, i.e., with the inadequacy of the LDA approach near the atomic limit . On the other hand, for moderately correlated systems such as, e.g., 3d metals, the main correlation effects are connected with spin degrees of freedom and can be described, in principle, basing on the LSDA electronic structure. Of course, it is difficult to say a priori where the boundary is between ‘moderately correlated” and “strongly correlated” systems. We will demonstrate that, at least for early actinides a LSDA-DLM description of the electronic structure turns out to be rather successful.

Due to the localized character of the spin moments, any ordered magnetic structure will resemble the result of the DLM description. It has often been noted that localization effects could partly be accounted for by way of exchange spin splitting of the 5f band, and despite an unphysical long range magnetic ordering (absent in the DLM approach), several spin-polarized ferromagnetic (FM) and anti-ferromagnetic (AFM) calculations have been performed. The presented DLM model gives a natural generalization to the paramagnetic state.

First we discuss some calculational details and then present the calculated volumes and bulk moduli for the actinide series. We then look at the results for $`\delta `$-Pu in detail, and the results are compared with different alternative approaches and experiments.

## II Calculational details

All total energies and densities have been calculated self-consistently within the framework of density functional theory , in the local-density approximation (LDA) in the non-magnetic cases, and within the local-spin-density approximation (LSDA) for the spin-polarized systems , with the local exchange-correlation functional by Perdew et al. .

We used the basis set of the $`s`$, $`p`$, $`d`$, and $`f`$ linear muffin-tin orbitals (LMTO) in the tight-binding representation and the atomic sphere approximation (ASA) for the crystal potential . The method was implemented within the scalar relativistic Green’s function technique . Spin-orbit coupling was not included. The disorder of local spin moments was treated within the coherent potential approximation (CPA), and other details relevant for the present calculation can be found in Refs. .

Since the interatomic distance is known to be the most important factor that determines the energetic of an actinide and since we are mainly interested in trends rather than in a detailed quantitative description of all the different phases of the actinides, calculations were performed for the fcc crystal structure only.

## III Disordered local spin moment picture of the actinides

Figure 1 shows a comparison between experimental and calculated Wigner-Seitz radii for the actinide metals. For the non-magnetic LDA calculation, we find a parabolic behavior typical of bonding through a series such as the transition metals. This is in full agreement with the experimental values for the early actinides, but fails do describe the volumes of the heavy actinides. Within the DLM picture on the other hand, we find a substantial improvement. For the light actinides, Th-U, results are identical to the non-magnetic LDA calculation, though for Np we find a slightly increased volume. For the later actinides we find an abrupt volume increase in close agreement with the experimental values with the transition taking place at Pu, which has an intermediate volume close to the value of the high temperature $`\delta `$-phase in the fcc structure. The results in Fig.1 are quite close to the data presented by Skriver et al. .

We find similar results for the bulk moduli. In figure 2 we display a comparison between the experimental bulk modulus for the ground state structures and the calculated bulk modulus in the fcc structure. The LDA calculations gives a fairly good description of the early actinides but fails completely for the later actinides. The DLM picture on the other hand gives a significantly improved result where, for example, the bulk modulus of Pu is reduced by a factor of two, though still somewhat higher than the experimental result for $`\delta `$-Pu.

Given the fact that the results are from parameter free ab initio calculations without considering the exact crystal structures and spin-orbit coupling the agreement between the calculated equilibrium volumes and bulk moduli with experimental values is remarkably good and the results clearly indicate the ability of the LSDA+CPA approach within the DLM picture to model important correlation effects beyond LDA in the 5f band, without incorporating long ranged magnetic ordering.

The mechanism behind the improved description of the equilibrium volume and bulk modulus is the formation of local disordered moments. This gives rise to a band splitting and partial localization of the 5f electrons which reduces the bonding resulting in an increased lattice constant and reduced bulk modulus. For Th and Pa local spin moments are still quenched, and for U only a negligible moment of about 0.07 $`\mu _\mathrm{B}`$ is formed. Within our DLM scheme, neptunium has a moment of 1.69 $`\mu _\mathrm{B}`$. Though this has little effect on the equilibrium lattice constant, the bulk modulus is reduced by a factor of two. For Pu a spin moment of 4.76 $`\mu _\mathrm{B}`$ is found in the DLM calculation, which corresponds to an almost complete spin polarization of the 5f electrons. For Am and Cm a moment of 6.57 and 6.90 $`\mu _\mathrm{B}`$ is formed, respectively. For Bk the disordered local moment is slightly reduced to 5.61 $`\mu _\mathrm{B}`$. Notice that fully relativistic calculations, including the spin-orbit interaction, would reduce the moments. In the case of $`\delta `$-Pu from about 5 $`\mu _\mathrm{B}`$ to less than 2 $`\mu _\mathrm{B}`$  and for Am the total moment can be expected to be $``$ 5 $`\mu _\mathrm{B}`$ . However, the formation of a J=0 atomic ground state comes natural from a completely localized 5f shell, and is consistent with experimental data.

## IV $`\delta `$-Pu

Figure 3 shows the calculated density of states (DOS) for Pu in comparison with the results of photo-electron spectroscopy . The LDA curve has qualitatively a very different behavior compared to experiments. For example, the peak at the Fermi level seen in the photoemission spectra is absent. The DLM DOS on the other hand, is in qualitative agreement with experimental data. The DOS is similar to the LDA-DMFT spectra by Savrasov et al.  (calculated at 600 K), but the quantitative agreement for the DLM picture is considerably better.

The electronic structure is strongly modified in the DLM picture compared to the LDA calculation. The 5f band is split, as shown in Fig. 4, due to the formation of local disordered moments, leading to an effective energy lowering of the occupied 5f band mass. The same effect is found also for the AFM calculation. This leads to a decreased bonding, increased equilibrium lattice constant and a reduced bulk modulus. The effect is similar to the results of LDA+U, LDA++, and LDA-DMFT. In contrast to some of these models we here have an improved agreement with the photo-electron spectra and the result is similar to what was found in the Mixed Level Model by Eriksson et al. . However, though a direct full comparison between the Kohn-Sham eigenspectra and the photo-emission data may not be possible, the data clearly indicates a resonance at the Fermi level which is not present in the LDA description.

In Fig. 3 we also display the density of states for an antiferromagnetic (AFM-LSDA) calculation. Though, resembling the general DLM 5f band split, the DOS is in less agreement with photoemission data.

## V discussions

There are two principle approaches for describing the $`\alpha \gamma `$ transition in Ce and the $`\alpha \delta `$ transition in Pu: the Mott transition model by Johansson  (for Pu, it has been considered in Refs.) and the Kondo-collapse model of Allen and Martin . In the context of DMFT, the appearance of the Kondo resonance is a feature of the Mott transition. An important question is thus what is more important for the energetics of the transitions: the appearance of the Hubbard gap or the formation of the Kondo resonance? The first feature is taken into account in the DLM scheme whereas the second one is not. It was already mentioned that the Kondo resonance has small spectral weight, so it is not obvious that ignoring it is an essential shortcoming. Moreover, both for Ce and Pu there are some arguments connecting the Mott transition to the peculiarities of the atomic electronic structure (atomic collapse ). These peculiarities can obviously be taking into account in the local spin-density approach and thus in the DLM description. Note that this collapse phenomenon, leading to a sharp dependence on the band width and lattice constant, leads to a first-order phase transition. Therefore the difference between the description of the Mott transition in the DLM or the Hubbard III approximation (continuous transition ) and DMFT (first-order transition, see Ref.  and Refs. therein) which is very essential for the Hubbard model, does not play a serious role for real actinide systems where the transition is essentially of first order for “quasi-one-electron” reasons.

As for the description of magnetic properties at high temperatures $`T`$, the static approximation gives a qualitatively correct picture. The magnetic susceptibility $`\chi `$ is proportional to $`<ϵ^2>/T`$, where $`ϵ`$ is the exchange on-site field. The susceptibility corresponds to the paramagnetic Pauli spin susceptibility provided that the exchange splitting fluctuations are Gaussian near the point $`ϵ=0`$, and to the Curie-Weiss susceptibility provided that spontaneous spin splitting exists. The description of finite temperature magnetism of transition metals in the DLM-CPA approach  is different from the DMFT treatment  only by the consideration of spins in a classical way; in the Heisenberg model, it corresponds to the appearance of the multiplier $`S^2`$ instead of the correct factor $`S(S+1)`$ in the Curie constant.

The DLM picture leads to an expected Curie-Weiss behavior of the magnetic susceptibility for $`\delta `$-Pu. However, recent experiments with Al and Ga stabilized $`\delta `$-Pu shows only a vannishingly small temperature dependence of the magnetic susceptibility, see Fig. 2 in Ref. . This is also consistent with other susceptibility measurements . The possible reason, explaining the absent Curie-Weiss behavior, could either be due to crystal field splitting and the formation of a non-magnetic multiplet  or because of the occurrence of a Kondo resonance that may diminish the temperature dependence of the susceptibility provided that the resonance width (“Kondo temperature”) is larger than $`T`$. At least the first explanation could possibly be included separately into the DLM picture, but the second one can be described only within the framework of DMFT. At the same time, it is important to realize that the intention with the DLM picture is to model some of the main characteristics of the energetics of the actinides, and it does not necessarily describe the magnetic properties correctly.

## VI summary

In summary we have presented a first principles Disordered Local Moment (DLM) method within the local density approximation with the disorder treated within the coherent potential approximation. The DLM picture gives an reasonably good description of bond lengths and bulk modulus for the actinide series. The equivalence between the DLM picture and the Hubbard III approximation and their relation to the DMFT description was discussed and it was argued for that the DLM picture is related to DMFT through a static approximation. The DLM density of states compares well with photoemission on $`\delta `$-Pu, in contrast to that obtained from the LDA or the magnetically ordered AFM configuration. In general, it is found that the DLM picture gives a considerable improvement over the LDA results and quantitatively good agreement with experiments.

## VII Acknowledgement

Discussions with D. C. Wallace, G. H. Landers and M. Colarieti-Tosti are gratefully acknowledged. I. A. A. is grateful to the Swedish Research Council (VR) and the Swedish Foundation for Strategic Research (SSF) for financial support.


format: text

text:
# Coupling time scales for simulation of structure transformation: an attempt to combine molecular dynamics and phase-field theory

## Abstract

A multiscale scheme combining molecular dynamics (MD) and microscopic phase-field theory is proposed to study the structural phase transformations in solids with inhomogeneous strain field. The approach calculates strain response based on MD and atomic diffusion based on the phase field theory. Simulations with the new technique are conducted in two examples. The first involves interface roughening in a Co/Cu thin film, where interfacial undulations due to lattice mismatch is demonstrated. The second example is a study of spinodal decomposition in AgCo/Pt/MgO(001) thin film, where we show that Co atoms are attracted by dislocations in the Pt/MgO interface, producing an interesting nanostructure.

Phase-field theory and molecular dynamics (MD) are both important theoretical tools to investigate the phase transformation in solids. In particular, phase-field theory has been frequently used to simulate the structure evolution of micro/mesoscopic morphology during phase transformations. A key factor to this capability lies in the flexibility of time scale used in the theory: one can always choose to trace the evolution of field variables whose dynamic evolution is “slow” compared to the remaining microscopic degrees of freedom. To incorporate the effect of elastic strains that are responsible for the formation of various mesoscale morphological patterns, continuum elasticity theory is usually used to calculate the strain energy induced by the lattice mismatch. However, a priori assumption on the transformation outcome has to be made because the elastic constants of the outcome phases are needed in the calculation. This may seriously reduce the predictive ability of the phase-field theory (especially the microscopic phase field model). In addition, the general formula of strain energy in various situations is often difficult to obtain. The second shortcoming of the phase-field theory comes from its assumption of coherent microstructure. It is not suitable for describing dislocations and their corresponding strain effect at an incoherent interface between two different phases. On the other hand, the method of molecular dynamics directly simulates the motion of moleculars by solving Newton’s equations. MD is highly effective in relaxing the local atomic distances and thus evaluating the stain energy. MD calculations are able to predict equilibrium and nonequilibrium properties of condensed systems. Despite of these advantages, MD has not been used to model dynamical morphological evolution in phase transformations. The main problem is that the characteristic time scale of MD is too short ($``$1 ns with time step $``$1 fs) compared with that of most phase transformation phenomena.

In this paper, an approach is developed to combine the microscopic phase-field theory and molecular dynamics to overcome the above difficulties. Two simple systems are investigated to illustrate how the new method works.

In MD, the variables to simulate are the atomic positions, {$`𝐫_i`$}, which evolve according to Newton’s equation:

$$m_i\frac{d^2𝐫_i}{dt^2}=_iV(\{𝐫_i\}),$$
(1)
where $`V(\{𝐫_i\})`$ is the empirical potential. For the microscopic phase-field theory, the variables considered are the site occupation probabilities, {$`n_i`$}. The value of $`n_i`$ represents the probability to find a solute atom at the $`i`$-th lattice site. The evolution equation is given as

$$\frac{dn_i}{dt}=\underset{j}{}L_{ij}\mu _j=\underset{j}{}L_{ij}\frac{F(\{n_i\})}{n_j},$$
(2)
where $`\mu _j`$ is the chemical potential at the $`j`$-th site and $`L_{ij}`$ is the kinetic coefficient proportional to the inverse average time of elementary diffusional jumps from site $`i`$ to $`j`$, and $`F`$ is the free energy function.

In our method, both the lattice position and the occupation probability, {$`𝐫_i,n_i`$}, are considered as variables. Here, $`𝐫_i`$ are used to describe local atomic motion such as strain response or dislocation, and $`n_i`$ are used to describe atomic diffusion between different lattice sites. With mean-field approximation, the empirical potential and free energy can be expressed as functions of both atomic position and composition, i.e., $`V(\{𝐫_i,n_i\})`$ and $`F(\{𝐫_i,n_i\})`$. For example, when the pair potential is considered, $`V(\{𝐫_i,n_i\})`$ in the current approach can be written as:

$$V(\{𝐫_i,n_i\})=\frac{1}{2}\underset{i,j}{}\left[n_in_jV_{AA}(𝐫_i𝐫_j)+2n_i(1n_j)V_{AB}(𝐫_i𝐫_j)+(1n_i)(1n_j)V_{BB}(𝐫_i𝐫_j)\right].$$
(3)
The evolution of atomic position {$`𝐫_i`$} is conducted according to Eq. (1) with the new potential $`V(\{𝐫_i,n_i\})`$. When the elastic equilibrium is reached, {$`𝐫_i`$} satisfy the following equations:

$$\frac{V(\{n_i,𝐫_i\})}{𝐫_i}=0.$$
(4)
Since the free energy depends on {$`𝐫_i`$} only via $`V(\{𝐫_i,n_i\})`$, one has

$$\frac{F(\{n_i,𝐫_i\})}{𝐫_i}=0.$$
(5)
Since diffusion is much slower than elastic response, we can assume that the system is in elastic equilibrium during the diffusion process. The diffusion equations are:

$`{\displaystyle \frac{dn_i}{dt}}={\displaystyle \underset{j}{}}L_{ij}\mu _j`$ $`=`$ $`{\displaystyle \underset{j}{}}L_{ij}{\displaystyle \frac{dF(\{𝐫_i,n_i\})}{dn_j}}`$ (6)
$`=`$ $`{\displaystyle \underset{j}{}}L_{ij}\left[{\displaystyle \frac{F(\{𝐫_i,n_i\})}{n_j}}+{\displaystyle \underset{k}{}}{\displaystyle \frac{F(\{𝐫_i,n_i\})}{𝐫_k}}{\displaystyle \frac{d𝐫_k}{dn_j}}\right]`$ (7)
$`=`$ $`{\displaystyle \underset{j}{}}L_{ij}{\displaystyle \frac{F(\{𝐫_i,n_i\})}{n_j}},`$ (8)
where $`dF(\{𝐫_i,n_i\})/dn_j`$ means that {$`𝐫_i`$} is function of $`n_j`$ and the differential with respect to $`n_j`$ should also operat on $`𝐫_i`$. The elastic equilibrium conditions Eq. (5) is used to obtain the final result of Eq. (6). In the numerical scheme, $`m`$ MD steps with time interval $`\delta t`$ are conducted to reach the elastic equilibrium state before every diffusion step with time interval $`dt`$. It is a multiscale scheme since $`m\delta tdt`$.

In order to illustrate how the method works in practice, we first consider interface roughening in a heteroepitaxial thin film system.

In uniform heteroepitaxial thin films, there usually exists large stress arising from the lattice mismatch between thin film and substrate. To reduce the strain energy, surface roughening by mass diffusion occurs during film growth or annealing. The length scale of surface roughening is about a few hundred nanometers, which is too large for the current atomic simulation. We alternatively simulate a similar process, the interface roughening, to demonstrate the effect of strain on the mass transport.

The system studied here is Cu-Co system where the lattice mismatch is about 2%. Cu and Co is immiscible in the bulk with a mixing energy of +13 kJ/g atom. Atomic interaction among Cu and Co atoms is described with the Tight-binding second-moment-approximation (TB-SMA) potential. The potential parameters for the pure species are available in the literature, while the cross parameters (Cu-Co) are evaluated by fitting the mixing energy and the lattice length of the disordered phase within the mean-field approximation. We consider a coherent Co film grown on Cu substrate. The as-grown film (initial state) has flat surface and interface \[Fig. 1(a)\]. The evolution of phase morphology is simulated by the above method, with result shown in Fig. 1(b). It reveals interfacial undulations with varying amplitude and wavelength. The undulation of interfacial increases the interfacial energy while decreases the strain energy. The total energy decreases in such process. In comparison with surface roughening, no groove or cusp is found in the current simulation of interface roughening.

As the next example, we consider the spinodal decomposition of a thin film on a periodically strained substrate.

Adsorbed layers epitaxially grown on single-crystal surfaces are generally strained due to the lattice mismatch. In the weakly incommensurate situation such as Ag layers grown on Pt(111) surface, highly ordered periodic dislocation network is formed to relieve the mismatch strain, hence invalidating the assumption of coherent microstructure usually adopted in the phase field theory. Such dislocation pattern can be used to fabricate interesting metal island arrays in two-dimension. If the second alloy film is grown on this strained film and the second film is unstable to spionodal decomposition, it may be expected that the dislocation network will act as a preferred length scale for the spinodal decomposition. To investigate this issue, we simulate the spinodal decomposition in a AgCo/Pt/MgO(001) thin film system.

For Pt films deposited onto a (001) MgO substrate, a dislocation network with periodicity 4.05 nm is predicted. We consider the growth of AgCo on the strained Pt/MgO thin film and simulate the spinodal decomposition of the AgCo phase. The potential between Ag-Co-Pt is adopted as TB-SMA potential. For simplicity, the interaction between Pt and MgO is approximated as the Lenard-Jones potential by fitting the adsorption energy and atomic distance. MgO atoms are kept fixed (with crystal lattice of 4.21 Å) in simulation. For the initial state, Pt atoms are arranged in their ideal lattice sites (with a lattice parameter of 3.924 Å), and a dislocation network between Pt and MgO naturally forms in the simulation. The initial AgCo phase is assumed to be coherent with Pt and is disordered with an average composition $`n_i=0.5`$. A simulation result of the thin film with 3 monolayer (ML) of Pt and 10 ML of AgCo is shown in Fig. 2. The atomic configuration with dislocation network and the corresponding strain effect are demonstrated in Fig. 2(a) and 2(b) for the initial stage before atomic diffusion occurs. It can be seen that large strain is concentrated near the dislocation cores. The final result of phase transformation is depicted as Fig. 2(c). It shows that Co atoms segregate in strained regions because they are smaller in size. Within a few layers near the AgCo/Pt interface, Ag atoms gather in less strained regions and produce an interesting nanostructure. The surface layers are occupied by Ag atoms since the surface energy of Ag (1.30 J/m2) is smaller than that of Co atoms (2.71 J/m2). It is interesting that some Pt atoms above the dislocation are “pulled up” by Co atoms in Fig. 2(c), which reflects the important interaction between Pt substrate and AgCo film. Such effect can not be described within the coherency approximation.

The method presented here is on the microscopic level. Similar to other atomic scale simulations, one of its disadvantage is that the system size simulated is relatively small compared to mesoscale morphology. In contrast, continuum field theories have more freedom in changing the length scale of simulation. Such theories can be expected to couple with computational mechanics methods such as the Finite Element method by a similar coupling approach as described here. This is not pursued in this paper.

This work was supported by a Max Planck Post-Doc Fellowship for ZL and by US National Science Foundation through Grant CMS-0085569. We thank Prof. Long-Qing Chen at the Pennsylvania State University for helpful discussions. The code of MD simulation in this work comes from the program GONZO of Stanford Multiscale Simulation Laboratory (MSL).


text:
# Shearing of loose granular materials: A statistical mesoscopic model

## I Introduction

Modeling the rheology of granular media using continuum solid mechanics, has reached a high degree of sophistication in terms of constitutive equations Wood . Whatever the complexity of load paths being studied, an accurate account of the experimental stress strain relationship can now be achieved provided enough parameters or internal variables are included in the constitutive laws. However, such approaches are descriptive and leave unanswered questions pertaining to the scale of grain sizes.

In parallel to such phenomenological descriptive theories, a lot of effort has been spent in recent years in developing powerful computer models able to simulate granular systems at the individual grain level. Molecular dynamics approaches LudingMD , or other techniques such as “contact dynamics” Radjai\_Brendel ; Radjaiforce , now offer the possibility of dealing with several thousand particles, and provide extremely realistic pictures of the detailed micro-mechanics.

Such numerical techniques can be used to accurately investigate displacement fields, resolved both spatially and temporallyRadj . The latter reveal an intriguing feature: namely that even in the most simple tests, such as a simple steady shear imposed over large strains, the local displacement appears very unsteady, with short quiescent periods where the displacement field is spatially smooth, separated by sudden changes where the configuration of grains reaches a local instability and undergoes a rapid reorganization through significant displacements at the grain level. This temporal variability manifests itself in giant stress fluctuations observed experimentally when particles and walls are stiff, and when the high frequency part of the stress signal is not filtered outBehr . Numerical simulations indicate More  that instantaneous strain fields consist essentially in localized strains occurring along one or a few shear bands. However, as the strain increases (over the moderate range accessible in the simulation), there seems to be little or no correlation between successive shear bands, so that the time average of the displacement field erases these discontinuities and produces smooth strain fields.

Such fluctuations are obviously ignored in continuum modeling. And indeed it may appear that the identification of these instabilities is relevant only for discussing fine details of microscopic and transient features. However, their relevance can be judged only at a mesoscopic level of modeling, since the microscopic numerical techniques are far from being able to reach the relevant time scales. In fact, in the following we will argue that these instabilities may have a significant impact, both on large scale heterogeneities of the medium itself, and on a systematic slow time evolution of the macroscopic friction angle. A short account of some of our results has appeared in a previous publication TKKRLett .

The paper is organized as follows: in Section II, we will recall some features observed experimentally or numerically that we consider essential, and, in Section III, we progressively introduce the rules of a model whose aim is to describe some statistical aspects of shearing of loose granular media over large strains. In Section IV, we present in details the different quantities studied numerically for this model. We conclude in section V with a summary of our results and a discussion on possible experimental checks.

## II The shear process in loose granular material

We will address here the question of the behavior of granular media subjected to a simple shear for large strains. We restrict ourself to the simplest granular medium one may consider, namely rigid (undeformable) grains with Coulomb friction. This refers experimentally to dry sand subjected to a low confining pressure. We are concerned here with large strains, and thus in order to avoid the problem of boundary conditions which would limit the maximum strain, we consider an annular shear cell. To simplify the problem further, we consider only the case where the problem is invariant along the shear direction. As shown in Figure 1, the displacement is a single function of the coordinate of a radial cross-section $`(x,y)`$, and constant along the orthoradial direction $`z`$ (traditionally this situation is termed “anti-plane”). Moreover, we are interested only in the quasistatic regime, i.e., time as such is irrelevant, and only the total strain matters. Thus, in what follows, what is referred to as “time” is to be seen here as a practical means of parameterizing the total strain being imposed on the medium.

One of the important observations of soil mechanics concerning such media is the concept of a critical state Critstate1 ; Critstate2 . Depending on the preparation of the sample, the behavior under shear may differ considerably. For loose sand, (low density), the deviatoric stress to be applied increases with the total shear strain and simultaneously, a densification is observed CompactExp ; LudingModel . However, as the shear strain increases, the density and shear stress seem to reach a plateau independent of the initial density. This state is called the “critical state”. On the contrary, if the initial density is large, a single shear band forms, while the rest of the medium remains frozen widthsh . The formation of the shear band is preceded by a volume expansion of the medium dil , but after the band is formed, all further properties remain quasi constant. A detailed experimental investigation has revealed Desrues  that inside the shear band, the density tends to approach the critical state density. This concept of the critical state has received considerable experimental evidence over the years, and is implemented in a number of continuum constitutive laws. Experiments however mostly deal with a rather moderate total strain well below unity.

A simple picture which is consistent with the critical state concept is that both the friction and the dilation angle increase with the density, and that the critical state is the density for which the dilation angle is zero (no change in volume under shear). Retaining the density as the only internal variable is an approximation. Other characteristics of the texture of the medium such as the fabric tensor (which has information about the orientation of contact normals), certainly play a significant role. For the purpose of simplicity, we will in the following only retain one single scalar internal variable governing the friction angle. It could be either simply the density or a combination of density and texture. Nevertheless, in all cases we will refer to this internal variable as “density”, irrespective of its precise meaning.

As mentioned above, numerical simulations seem to reveal More  evidence for the existence of instantaneous shear bands even in loose granular media. On the other hand, in experiments the strain appears to be homogeneous and not localized. The resolution of this apparent paradox is that the shear bands change rapidly, and may visit the entire medium in the process. Thus during an increment of shear, which can be observed experimentally, only a time average over many such shear bands is seen. In our modeling, we introduce a basic time scale for each elementary procedure. This time step is then clearly much shorter than most experimentally accessible time scales. However, our model attempts to achieve a qualitative rather than a precise quantitative mapping. One of the main features of the model is to show that these two apparently unrelated facts: the existence of instantaneous shear bands at early times and its localization at late times, are actually related, with a slow transition between these two limiting cases. This slow dynamics is reminiscent of slow ageing properties encountered in glassy systems, and indeed, we will see that a breakdown of ergodicity does appear in this model.

## III The model

### III.1 Motivation and definition

At every instant the two-dimensional medium is characterized by a single, scalar internal variable, the density $`\varrho (x,y)`$. This represents an average of the density along the orthoradial direction $`z`$. From this density, we deduce a corresponding local friction coefficient $`\mu (x,y)`$. The latter is assumed to be a single monotonically increasing function of the density JanoPull . For simplicity, we may assume a linear relationship in the following although this is inessential.

The strain is imposed on the shear cell through prescribed displacements of the bottom and top planes. As particles are considered rigid, (no elastic deformation), the shear cell can only move if the shear force exceeds a threshold value proportional to the normal pressure. This limit stress is given by the “weakest internal surface”. Indeed, in our anti-plane geometry, the shear strain will localize on the surface (i.e. path in the $`(x,y)`$ plane) which will fail first. The latter is assumed to be given by the following algorithm. For each directed path $`𝒫`$ spanning the entire cross-section along the $`x`$-axis, we compute the maximum shear force it can support according to the local density. Assuming that the local slope of the path is always small, this maximum force $`F(𝒫)`$ is simply proportional to the sum of local friction coefficients, and thus making use of the assumed linear variation of the friction coefficient with the density, $`F(𝒫)`$ is proportional to the sum of local densities,

$$S(𝒫)=\underset{(x,y)𝒫}{}\varrho (x,y),$$
(1)
where the sum runs over the sites along the path. Among all the possible paths, the weakest $`𝒫^{}`$ (for which $`S(𝒫^{})=\mathrm{min}`$) will fail first, and this fixes the value of the shear force $`F=F(𝒫^{})`$. In agreement with the previously mentioned observation, at every basic time step, the shear strain is realized along a single shear band. Away from this shear band, the strain rate is zero, and thus the density is kept constant in time. However, inside the shear band, there is a relative motion of grains, and thus the density is susceptible to evolve.

The next step is now to determine how the density inside the shear band evolves with time. Even though it is observed that large strains are necessary in order for a system to reach its critical state, we argue that at the microscopic level, the evolution of the medium cannot depend on the total imposed strain. Thus the evolution rules for the density within the shear band should be designed in such a way that they do not depend on the past history, but only on the present state (density field). As the density ought to contain the basic information of the local characteristics, we propose that within the shear band, the local density $`\varrho (x,y)`$ is randomly modified. More precisely, in one elementary time step, corresponding to the “life-time” of the shear band in a very loose granular sample, we assume that the density along the shear band acquires random uncorrelated values picked from a statistical distribution $`p(\varrho )`$. The uncorrelated character of the distribution is however justified only on a mesoscopic scale.

After the elementary strain event, we have a new density map $`\varrho (x,y)`$. We now simply reiterate this procedure as long as desired: Namely we identify the new path which minimizes $`S(𝒫)`$, and update the value of the densities along this path randomly. As the purpose of the present article is to illustrate some statistical aspects of this dynamics, we do not try to mimic any specific granular system by imposing a realistic density distribution or initial correlations in $`\varrho `$. We will choose here a simple uniform distribution between 0 and 1 for $`p(\varrho )`$. The mean value and variance of the distribution $`p`$ can be chosen arbitrarily, since a translation and rescaling of $`\varrho `$ does not affect the result.

A key assumption of our model which may appear as precluding the occurrence of a slow evolution toward a critical state is the selection of the density values within a shear band from uncorrelated, smooth distributions. In fact, we will show below that, on the contrary, a collective and purely statistical effect produces a slow increase of the mean density over large strains.

Our model is furthermore discretized on a regular square lattice. We have looked at two different kinds of square lattices to check the robustness of our results. In the first one the value of the density $`\varrho `$ is carried by the bonds. The orientation of the lattice is chosen so that the principal directions lies at $`\pi /4`$ with respect to the $`(x,y)`$ axis as shown in Fig. 2 a). In the other version, density values are assigned to the sites of a square lattice. In this case the minimal path can be connected through the next nearest neighbours too as shown on Fig. 2 b). Both square lattices give exactly the same results so in the following we just refer to them as square lattice realizations.

We mention here that a third type of lattice, the hierarchical diamond lattice ber  was also studied. The numerical results are surprisingly, essentially unchanged by the unusual topology of this recursively constructed lattice. The easier construction of this lattice however allows us to solve the model analytically thus giving us a quantitative picture of the behaviour of the system. These results are presented in TKKRhier  where the intimate relationship of our model to other models of statistical physics is also discussed.

The rules of our model, finding the extremal directed spanning path at every instant, is similar to finding the ground state of a directed polymer in a random potential HalHeaRev . However, in our case this potential is uncorrelated only at the beginning; it changes in time through the process described above, of ascribing new densities to all sites along the minimal path. It is clear from this relationship between the models, that the shear band is self-affine with a Hurst exponent $`\zeta =2/3`$ at the beginning, i.e. the transversal fluctuations of the band grow with the size ($`L_x`$) of the sample width as $`L_x^\zeta `$. We will discuss the time evolution of the roughness later in this paper.

## IV Numerical results

We first show the density map, Fig. 3, of the system at different times, $`t/L`$ ranging from 4 to 4000. The grey scale chosen focuses on the vicinity of 1 so as to highlight the progressive densification. It may appear counter-intuitive at first that the rest of the medium shows a densification at all, when the only dynamics consists of finding a minimal path and updating sites along it randomly. However the reason is simply that this update systematically hunts out the sites with the lowest density values and replaces them.

At early times, we observe an apparently uncorrelated field. However, as time proceeds, it is possible to distinguish preferential channels of high density aligned along the direction ($`x`$-axis) of the minimal paths. These channels however have a significant width which shows that though the minimal path has been confined to this zone, it has enough freedom to explore different neighboring configurations and achieve a significant local densification. We also see within these wide and dense channels, a single path with a lower density. This has been the last active minimal path in the channel. As time passes, the number of channels increases, and so does their width. They get partly interconnected, leaving always the same scars of low density paths. Finally, at the latest time shown on the figure, the average density is quite high, and traces of ancient minimal paths are still visible. Nevertheless what is striking is the occurrence of islands entrapped by these high density channels, where the density map looks like at the very early stage of the time evolution. This signals that these regions have basically not been visited by the minimal path during the entire history of the system. These features are quite generic, and they reveal that the spatial (and temporal) organization of the activity is rather complex. The rest of the study is devoted to a more quantitative account of this activity, of the resulting kinetics of compaction, and the unexpected finite size effects which appear in this problem.

In the following subsection, we will introduce the main measurements performed numerically on the model.

### IV.1 Definitions of numerically measured quantities

#### IV.1.1 Average density

The most important quantity is the average density of the sample that we define as the mean of the density of the inactive sites, i.e. the sites not belonging to the shear band. We denote this by $`\varrho `$. This definition is convenient because $`\varrho `$ increases monotonically by the rules of our model. In experiments one of the most frequently measured quantities is the volumetric strain which is just the change in the inverse of the average density. In our model, when $`p(\varrho )`$ is chosen to be the uniform distribution in the interval $`[0,1]`$, the average density is bounded by ($`\varrho 1`$). It is easy to see that in finite systems the steady state (the asymptotics) cannot be other than a system with maximal densities everywhere except for a path which will be always chosen as the minimal path. This state is equivalent to $`\varrho =1`$.

We are interested in the approach to this asymptotics so we plot the quantity $`(1\varrho )`$ as a function of time. Within the granular medium context, this means that we mainly study a loose initial state and its convergence to the critical state Critstate1 ; Critstate2 . However, we will also present results obtained when one starts from a high initial density later (Section IV.7).

Figure 4 shows that the difference of the average density from its asymptotic value first remains almost constant during a first stage $`t/L1`$, and then it decreases steadily. This first increase of $`\varrho `$ is well captured by a reduced time equal to $`t/L`$. However, as time progresses, the average density increases more and more slowly. Quite strikingly, the larger the system size, the slower the increase in density. Further down this will be interpreted as a breakdown of ergodicity.

#### IV.1.2 Shear band density

It is also natural to define the density of the shear band that we denote by $`\varrho _{SB}`$. This is just the average density of the sites along the minimal path (before updating them). As already mentioned in the introduction, we assume that the maximal static shear force is a single function of the density. Thus the density of the shear band can be related to the shear stress in experiments.

Figure 5 shows the evolution in time of the difference $`(0.5\varrho _{SB})`$. As expected, as time proceeds, the density along the shear band will tend toward the average of the random densities which are used to refresh the sites or bonds along the shear band. Using a uniform distribution of densities between $`0`$ and $`1`$ implies that this average is $`0.5`$. Note that in contrast to the previous case, the reduced time $`t/L`$ accounts nicely for the time evolution of this quantity for all system sizes for $`t/L10^4`$. This is a second puzzle we will try to address further in the following sections as well as in TKKRhier .

#### IV.1.3 Mean Hamming distance

We also calculate the Hamming distance of two successive shear bands which is defined as the number of sites (or bonds) by which two consecutive shear bands differ. We denote this distance by $`d`$. The natural normalization is to divide this distance by the total length of the path, $`L_x`$. As we will see this quantity is very useful in characterizing the time evolution of the localization process.

Fig. 6 shows that the mean Hamming distance is close to unity (i.e. two consecutive paths do not overlap at all) at early times, and decrease towards $`0`$ for $`t/L1`$. We recall that when the distance is equal to $`0`$, then the two consecutive conformations of the shear band are identical, in spite of the total renewal of random densities along them. This indicates that the shear bands have a tendency to remain more and more persistent as the system “ages”. We will analyze further the complete statistical distribution of the Hamming distance later in this paper.

#### IV.1.4 Cumulative shear

An experimentally relevant quantity is the cumulative shearing denoted by $`\sigma _{cum}`$. The numerical procedure we follow to obtain this quantity in our model is the following (see Fig. 7): We mark a line in the $`y`$-direction (see Fig. 2). We measure the total activity $`n_a(y)`$ along the line, i.e. the number of instances when the shear band went through a point $`y`$ on the line. From this, we define

$$\sigma _{cum}(y)=\underset{j=0}{\overset{y}{}}n_a(j).$$
(2)
By definition $`\sigma _{cum}(L_y)=t`$, since at every instant, the shear band has necessarily to pass through one or the other site on a cut along the $`y`$ axis. The fluctuations of $`\sigma _{cum}(y)`$ about its mean value $`(t/L_y)y`$ then reflect the inhomogeneity of the shear process within the sample. In Fig. 8, we track the time evolution of $`\sigma _{cum}(y)`$, after subtracting out the mean value. As indicated, a snapshot of this quantity encodes the history of the process of shearing in this system.

### IV.2 Early time regime

It is apparent from Figures 4, 5 and 6 that the initial behavior of the model is very different from the late stages. In the former regime, the average distance $`d`$ (Fig. 6) is equal to the system size indicating that successive shear bands do not overlap at all. In other words there is an effective strong repulsive interaction between them.

The density on the shear band (Fig. 5) provides an explanation for this behavior. From the directed polymer HalHeaRev  picture it is known that the first shear band has a mean density of $`\varrho _{SB}(t=0)=e^{}0.22`$ on the tilted square lattice. The state of the shear band after the densities of the sites along it have been changed, is independent of its previous state and the mean density is now $`0.5`$ (since it is just an average of $`L_x`$ independent random numbers taken from the uniform distribution ). This implies that at the next instant, the chosen shear band will have a very low probability of sharing bonds/sites with the previous one, since there will still be many spanning paths with density smaller than $`0.5`$. Thus, until all the sites are visited at least once, successive shear bands have few sites/bonds in common ($`d=L_x`$ in Figs. 6). Also in this regime, since the path sweeps all sites through the cumulative shearing appears to be homogeneous (Fig. 8 first plot). As mentioned in the introduction, this uniform shear strain is consistent with the experimental observations that no well-defined shear-bands persist over observable time scales in loose samples.

The characteristic time to build up correlations is set by the sweeping through the sample, i.e., it is given by $`L_y`$. This is the reason why in the early time regime the plots in Figs. 4-6 can be scaled together with $`L_y`$.

### IV.3 Localization

The above described behavior is drastically changed as time goes on and the shear band gets localized for very long times at the same location.

As the average density increases, the probability of choosing a minimal path with density less than $`0.5`$ starts decreasing, and thus for the new path, it becomes more favorable to overlap to a greater and greater extent with the previous one. As a result, the density of the region where the path is located, is further increased, hence trapping the minimal path in canyon-like structures surrounded by extremely high density regions (see Fig. 3 b-d).

The early repulsive interaction is thus now inverted to an attractive one as is shown by the rapid $`1/t`$ like decrease of the Hamming distance (Fig. 6), as well as the decrease of $`0.5\varrho _{SB}`$ and of $`\varrho `$ from its initial value (Figs. 5 and 4 respectively) in this regime.

There are a number of consequences of the localization. First, as time goes on, the shear band (which was not visible at all in the density map of the system (see Fig. 3 a)) becomes more and more apparent until finally it becomes localized at a given position for macroscopic times. This is in accordance with experiments and with the critical state concept Critstate1 ; Critstate2 . Secondly, since the same path for the shear band is chosen most of the time, its density saturates to its asymptotic value $`0.5`$ (Fig. 5). Again this is consistent with the experimental observation that the density within the shear band tends to achieve a well defined value somewhat smaller that the rest of the medium for dense granular media. Simultaneously, the shear stress saturates to a constant value since this is imposed by the shear band itself Wood . However, we note that in our model the global density of the system continues to increase in time, albeit very slowly. This extremely slow trend may well be out of reach experimentally. However, it is known Wood  that the shear stress saturates much faster than the volumetric strain (the average density in our case) which is clearly justified by the numerical results.

On the cumulative shear which is a straight line with small statistical fluctuations in the early time regime, there appear step-like structures, indicators of progressively more persistent localization (Fig. 8). However, this localization is not everlasting since the shear band may perform big jumps to other local minima. This can be seen on the series of cumulative shear curves in the form of certain steps disappearing and others becoming more prominent. This prediction could easily be tested experimentally.

### IV.4 Systems with different aspect ratios

All changes take place along the shear band, which is aligned along the $`x`$ direction, and thus we may anticipate that the $`x`$ and $`y`$ direction will play different roles. Therefore in this section, we study the influence of the width and length of the system. In what follows, we use the terms long for samples with $`L_y<L_x`$ and wide in the opposite case ($`L_y>L_x`$).

With very long and wide systems we are able to separate the two kinds of dynamics described in Section IV.3. If one considers a long sample ($`L_y/L_x`$ is small), we expect the time evolution to be independent of $`L_x`$ for all quantities of interest since the lattice can be split into subparts placed in series. Thus one may expect the large jumps to disappear and the average density to scale solely with $`L_y`$ (Fig. 9). A wide system, on the other hand, might be expected a priori to behave like a number of competing subsystems.

#### IV.4.1 Long systems

On Fig. 9 we have plotted the time dependence of the density in long samples with $`L_y=4`$ (lower curves) and $`L_y=8`$ (upper curves). The $`t/L_y`$ scaling is excellent in both cases. However, the densification obeys a different time evolution for different $`L_y`$. The rate at which the density increases is slower as the width increases.

#### IV.4.2 Wide systems: Breakdown of ergodicity

Wide samples can be considered as subsystems placed next to each other and coupled in parallel. In contrast with the previous case, we will see that the evolution of the different subsystems cannot be accounted for by a simple average.

In the case of the wide systems the same plot as before (Fig. 10 a) shows no $`t/L_y`$ scaling for small system sizes.

We could imagine the following construction: suppose we split a given wide system into two subsystems of size $`L_x\times (L_y/2)`$. Provided $`L_y`$ is large enough, we can ignore the interaction between the two subsystems, and thus we could study independently the time evolution of both subsystems. Now, if we merge them again, we realize that the only reason why the resulting densification could differ from the measurement on the separate subsystems is that the time $`t_1`$ the shear band has stayed in subpart $`1`$ is far from being equal to $`t/2`$. In other words, the breakdown of the data collapse of $`1\varrho `$ vs. $`t/L_y`$ is a breakdown of ergodicity. This is naturally associated to what we termed “localization” earlier.

In order to make this concept more explicit, let us consider an extreme version of such a breakdown of ergodicity. During a first stage, up to $`t/L_y`$ of order 1, the activity is evenly spread over the system. Then we assume that after such a time, the activity remains confined in a subsystem of size $`L_x\times \mathrm{}`$. All other subsystems are assumed not to be visited by the shear band, and thus their density is quenched at their value reached at the onset of localization, $`\varrho _0`$, at time $`t_0=\theta L_y`$. The global density will thus obey

$$1\varrho (t)=\frac{(1\varrho _0)(L_y\mathrm{})+f_{\mathrm{}}(t\theta (L_y\mathrm{}))\mathrm{}}{L_y}$$
(3)
where $`f_{\mathrm{}}(t)=1\varrho (t)`$ describes the densification of the representative cell of size $`L_x\times \mathrm{}`$ where the activity is confined, and thus $`1\varrho _0=f_{\mathrm{}}(\theta \mathrm{})`$.

In this crude scenario, we note that the global density does not converge to $`1`$ as time goes to infinity, but rather remains stuck at a value such that $`(1\varrho (t))f_{\mathrm{}}(\theta \mathrm{})(1\mathrm{}/L_y)`$. In more quantitative terms, we tried to carry out such a procedure, and indeed for a fixed $`L_x`$, it is possible to account for the time evolution of systems of different width using $`\varrho _0`$ as free parameter and calculating $`\mathrm{}`$ from the $`L_y`$ dependence. It turns out that $`\mathrm{}`$ changes for small values of $`L_x`$ but becomes constant ($`\mathrm{}30`$) above the system size of $`L_x30`$. The test of this analysis can be seen on Fig. 11. However, the asymptotic density turns out to depend on $`L_x`$.

The conclusion is that although, such an extreme modeling of the localization is able to capture part of the strong size effects observed numerically, it is too crude to provide a quantitative account of the densification. The hierarchical lattice provides us with a convenient case where an analytical investigation of this breakdown of ergodicity can be made. It is shown in Ref. TKKRhier , that the local “age” distribution assumes a multifractal distrbution whose spectrum can be computed exactly. This property can then be used to provide an expression of the density evolution in time.

### IV.5 Time evolution of the Hamming distance distribution

We study here the distribution of Hamming distances as a function of time, $`P(d,t)`$ for both lattices. This quantity is the analogue of an “avalanche distribution”, such as is usually studied in self-organized critical(SOC) systems. As we will see below, this quantity does indeed decay as a power-law like in SOC systems. However, since a steady state is never reached, the power-law decay is multiplied by a time-dependent prefactor.

The distribution $`P(d,t)`$ is shown in Fig. 12. At early times, this quantity is peaked around the maximum value ($`d=L_x`$) while in the localized regime, it becomes peaked at the minimum value ($`d=0`$). This corresponds to the transition from repulsive to attractive effective interaction between consecutive path conformations.

At fixed (large) times, the distribution $`P(d,t)`$ decays as a power-law of $`d`$, as can be seen in Figure 13. The measured exponent is

$$P(d,t)d^3$$
(4)
in addition to which there exists a peak at $`d=0`$ the amplitude of which varies significantly with time. The decay of average $`d`$ as $`1/t`$ (see Fig. 6) found earlier, implies that the time dependence of the $`d0`$ part is $`p(d,t)1/(td^3)`$. Thus, including the different scalings with $`L_x`$ and $`L_y`$, we obtain finally the asymptotic form

$$P(d,t)\frac{L_xL_y}{td^3}$$
(5)
Here again, the hierarchical lattice allows us to compute this distribution analytically (for $`L_x=L_y`$), and we find that a similar behavior is obtained TKKRhier .

### IV.6 Roughness exponent of the shear band

From the directed polymer analogy we know that the shape of the shear band is self-affine with an exponent of $`\zeta =2/3`$ for infinitely large systems. It is an interesting question whether this property of self-affinity is conserved in the time evolution of our system. We have investigated this question and have found self-affine scaling for all times, albeit with a time dependent Hurst exponent (Fig. 14).

We have estimated the value of $`\zeta `$ by measuring the width of the shear band, $`w_L(t)`$, for different system sizes, $`L\times L`$. The width is defined as the standard deviation of the $`y`$ coordinate of the active path. The latter is expected to scale as $`w_L(t)L^\zeta `$ for a self-affine object. To estimate $`\zeta `$, the roughness exponent, we have computed the ratio of two such widths for lattice sizes differing by a factor of 2, and used the following estimate

$$\zeta =\frac{\mathrm{log}\left(w_{2L}(t)/w_L(t)\right)}{\mathrm{log}(2)}$$
(6)
where $`w_L(t)`$ is the width of a shear band in a $`L\times L`$ lattice at time $`t`$. The results obtained for the tilted square lattice can be seen on Fig. 14. It starts from $`\zeta (t=0)=2/3`$ as expected from the directed polymer result and has an asymptotic value of $`\zeta 0.8`$.

### IV.7 Systems with high initial densities; relevance of initial conditions

We have so far studied the situation when a very loose granular medium compactifies under shear, while simultaneously a shear band gets quasi-localized in the system. It is also of interest to study samples with higher initial density where it is known experimentally that a shear band is localized from the very beginning. In this section we study the interesting crossover to that state from the previously described dynamics.

We choose initial densities from the interval $`[\varrho _{init}:1]`$ with a uniform distribution and varying $`\varrho _{init}`$. In Fig. 15 the time evolution of the difference of the average density from its asymptotic value is plotted using different initial conditions. Since all previous arguments hold we assume that in finite systems the asymptotic value of the average density is $`1`$ irrespective of the initial density.

The striking result of these simulations is that if $`\varrho _{init}<0.5`$, all curves coincide in the decreasing regime. However, if $`\varrho _{init}>0.5`$ a different time evolution is observed for large times. Thus we can assume that in the early time regime, when the shear band is swapping uncorrelatedly, it visits all sites that have a value less than the expectation value of the refreshing density distribution $`=0.5`$. After the first regime, as all small values are eliminated, the system effaces the initial condition almost entirely.

On the other hand, starting from initial conditions with $`\varrho _{init}>0.5`$ we largely eliminate the possibility of big jumps. The shear band fluctuations are now very small, involving changes in a very few sites, and hence, the density change of the sample is extremely slow.

### IV.8 Summary of the numerical results

We have seen that the system densifies with time so as to approach a unit density, i.e. the maximum available density from the distribution used to refresh the sites. The kinetics of the densification is slow (slower than any power-law). Moreover, after a first transient where the reduced time $`t/L`$ accounts for the $`L`$ dependence, the compaction process depends on the system size in a non-trivial way. The width of the system is the parameter which really controls this anomalous behavior, signaling that the competition between parallel paths may somehow play a key role in this breakdown of system size rescaling. This competition is a subtle one however, lying somewhere in between complete localization of the path (which, as we saw in section IV B, is too crude to mirror the actual scenario) and complete delocalization (which, as mentioned earlier, accounts only for the early time behaviour).

The density maps display an interesting organization of “canyon-like” paths with density much lower than their immediate surroundings (where the density approaches 1 quite uniformly). Moreover, large regions are left quiescent, being systematically avoided by the minimal paths. This contrast of high and low activity within the same system is at the heart of the breakdown of ergodicity observed after an initial transient. As remarked however, the distribution of Hamming distances between consecutive paths displays a somewhat simpler behavior where the role of the width and the length of the system can be simply accounted for.

## V Conclusion and discussion

Though very simply defined, our model seems to capture some essential features of granular shear and provides at the same time several predictions. The model demonstrates the self-organized mechanism of the localization of the shear band in loose granular materials. As the sample ages, very high fluctuations in density appear where we can observe some kind of screening effect: more resistant regions of higher density, protect the looser ones. This model also gives an insight into a dynamics that exhibits very non-trivial system size effects.

Our stochastic model makes an attempt to describe the large strain behavior of sheared loose granular matter on a mesoscopic level. The rules of the model do not include any dependence on the total amount of shear imposed on the medium, nevertheless, a constant friction angle and slow densification is observed — a property referred to as “ageing” — which reproduces the experimental results qualitatively CompactExp . By construction, the strain takes place through local shear bands which initially travel throughout the medium homogeneously (and hence produce a uniform shear), but which progressively become more permanent giving rise to more steady shear bands, a feature also observed experimentally widthsh . Our model reproduces further features seen in experiments and numerical simulations, including the high frequency fluctuations of the local shear Behr .

In addition, we predict a complex self-organization of these shear bands, displayed in the inhomogeneities in the local density. This feature can be studied experimentally, in particular through the use of X-ray tomography, to access the local density of a sheared medium. The use of tracer particles could also be helpful in identifying the inhomogeneous ageing and localization of the shear bands as well as their sudden changes.

Most of the results presented here for the Euclidean lattice are closely mirrored by the results on the hierarchical lattice studied in TKKRhier . The recursive topology of this lattice allows a quantitative analytical understanding of many of the quantities studied numerically in this paper. This includes elucidating the mechanism for the breakdown of ergodicity, the slow density evolution, as well as the behaviour of the of the Hamming distance at late times.

## Acknowledgments

This work was supported by EPSRC, UK, OTKA T029985, T035028.


text:
# Measurements of Compton Scattered Transition Radiation at High Lorentz Factors

## 1 Introduction

Space-borne cosmic ray experiments require the capability to measure the energies of particles with Lorentz factors $`\gamma 10^5`$ with detectors that are relatively large yet lightweight. NASA’s proposed Advanced Cosmic Ray Composition Experiment for Space Science (ACCESS) mission , for example, requires a transition radiation detector (TRD) capable of measuring the energies of cosmic rays up to 100 TeV/nucleon for particles with charge $`Z>3`$. Such experiments require that the range of existing TRDs must be extended upward by an order of magnitude or more, requiring designs modified for use at these higher energies.

Transition radiation (TR) is produced when a charged particle crosses the interface between two materials with different dielectric constants, resulting in the rapid rearrangement of the particle’s electric field as it passes from one material to the next . For highly relativistic particles ($`\gamma =E/mc^21`$) the radiation is emitted at x-ray frequencies. The spectrum produced depends on the plasma frequencies and thicknesses of the two materials as well as the energy of the particle. Typically, the materials used are a low atomic number solid such as plastic with plasma frequency $`\omega _1`$, and a gas or vacuum with plasma frequency $`\omega _2`$. Radiation is emitted up to a frequency $`\gamma \omega _1`$, beyond which the spectrum is suppressed. In the usual case, where $`\omega _1\omega _2`$, the total intensity produced from a single interface is proportional to $`Z^2\gamma \omega _1`$, where $`Z`$ is the atomic number of the incident particle.

The intensity from a single interface is weak. Therefore, in practical applications, a radiator is constructed with a large number $`N`$ (typically $`N501000`$) of thin foils of thickness $`l_1`$ separated by a distance $`l_2`$ (or fiber or foam radiators with equivalent average $`N,l_1,`$ and $`l_2`$) with radiation produced at each of the $`2N`$ interfaces. Interference effects from the superposition of the amplitudes produced at each interface give rise to pronounced minima and maxima in the spectrum, with the last (highest frequency) maximum near

$$\omega _{\text{max}}=\frac{l_1\omega _1^2}{2\pi c}(1+\rho ),$$
(1)
where $`\rho `$ is 1 for a metal and 0 for a nonconductor. As the particle energy increases, the total radiated intensity increases up to a Lorentz factor

$$\gamma _\text{s}\frac{0.6\omega _1}{c}\sqrt{l_1l_2(1+\rho )},$$
(2)
above which saturation sets in due to the interference. We have included here the possibility of a nonzero conductivity which introduces an imaginary part to the wave vector and leads to an effective plasma frequency $`\omega _1\sqrt{(1+\rho )}`$ . The saturation energy and characteristic frequency can be tuned by varying the radiator foil material, thickness, and separation.

An x-ray detector appropriate for absorbing the TR x-rays must be placed after the radiator. The radiation is emitted at an angle $`\theta 1/\gamma `$ with respect to the incident particle direction, so for high Lorentz factors the x-rays are spatially inseparable from the ionization energy deposited in the detector by the particle itself. Therefore, in conventional applications, the detector must be made thin in order to minimize the ionization signal, yet with sufficient stopping power to absorb the x-rays. For $`\omega _{\text{max}}`$ less than about 40 keV, gaseous detectors (e.g. Xenon-filled wire chambers) are typically employed. In order to improve statistics and provide redundancy, a complete TRD consists of multiple layers of radiators and x-ray detectors.

Such TRDs have been used successfully both at accelerators and in space. In most cases, the TRD is employed as a threshold device to identify particle types: For example, a meson or hadron may be accompanied by a small TR signal, while an electron of the same energy but larger $`\gamma `$ is characterized by a large signal. In the case of the Space Shuttle CRN experiment , a fiber TRD was used to measure the energies of cosmic ray nuclei with $`\gamma 3\times 10^3`$. Ref.  gives an extensive review of TR applications and radiator configurations. A brief listing of cosmic ray experiments incorporating TRDs is given in .

In order to increase the maximum particle energy $`\gamma _\text{s}`$, one must increase the plasma frequency (or equivalently, density), thickness, and/or spacing of the foils (Eq. 2). In a space instrument, the overall thickness will be constrained, putting a limit on $`Nl_2`$ (assuming $`l_2l_1`$). Increasing $`\omega _1`$ by using metal foils instead of plastic, for example, and/or $`l_1`$ results in a hardening of the x-ray spectrum produced (Eq. 1). Metal foils have been used in early accelerator tests , and in particular lithium foils have been used in order to minimize the absorption at low x-ray frequencies . In the case of very high energies, though, with $`\gamma _\text{s}10^5`$ and a typical spacing $`l_2=0.11`$ cm, $`\omega _{\text{max}}0.4\gamma _\text{s}^2c/l_2`$ can be in excess of several hundred keV. Gas detectors are then no longer efficient in detecting these hard x-rays. Although  describes the use of gas detectors near $`\gamma _\text{s}10^5`$ by optimizing the radiator design, scintillators such as NaI or CsI provide an efficient alternative at these Lorentz factors and corresponding high x-ray energies. The higher density of the scintillators leads to an increase in the ionization energy deposited by the particle as it traverses the detector. However, as the TR spectrum hardens, Compton scattering in the radiators becomes important, becoming the dominant photon interaction above $`40`$ keV. A significant portion of the x-rays produced are scattered out of the path of the incident particle. Thus, a detector that is segmented or positioned outside of the beam can efficiently detect the TR signal spatially separated from the ionization.

We describe here the test of a scintillator-based Compton Scatter TRD for high Lorentz factor particles, including the use of metal foils, based on the results of accelerator measurements with high energy electrons at the CERN SPS. Comparisons of the measured results with detailed simulations will also be presented.

## 2 Experimental setup

A scintillator-based Compton Scatter TRD was designed to investigate the predicted increase in saturation energy obtained by using thick, dense radiator materials including Mylar ($`\rho =1.4`$ g/cm3), Teflon ($`\rho =2.0`$ g/cm3), and aluminum ($`\rho =2.7`$ g/cm3). Metal foils are of particular interest because of the characteristic enhancement in the signal expected due to the nonzero conductivity .

For the plastics, radiators of $`N=50`$ foils were constructed by attaching 19.1 cm $`\times `$ 18.4 cm plastic foils to 3.4 mm thick wood frames and stacking them together. For the aluminum, each radiator consisted of seven 2.7 cm thick honeycomb panels bundled together and aligned with the cells perpendicular to the particle beam. The honeycomb was a composite material chosen both for its dimensions and its adaptability as a combined detector-plus-structure for a space instrument (Fig. 1). Particles passing through the structure passed through either 1) a section of foils perpendicular to the beam in which two 3 mil sheets glued together form a foil with an effective $`l_1=6`$ mil, $`l_2=5.2`$ mm, and $`N=35`$ foils along the particle trajectory or 2) a section of foils at a $`41^{}`$ angle with respect to the particle beam, resulting in an effective foil thickness $`l_1=3\mathrm{mil}/\mathrm{sin}41^{}=4.6`$ mil, $`l_2`$ ranging from 0 to 5.2 mm, and $`N=70`$. The yield from a composite material (e.g., a foam) with average values $`l_1`$, $`l_2`$, and $`N`$ has been shown to be essentially the same as from a regular foil radiator with the same $`l_1`$, $`l_2`$, and $`N`$ . We therefore calculate the Al honeycomb effective parameters as averages of configuration 1 weighted by 46% (to account for the fraction of the area perpendicular to the beam covered by configuration 1) and configuration 2 weighted by 54%. The resulting average honeycomb radiator parameters are given in Table 1, along with the parameters of the plastic radiator configurations tested. The total length of each radiator was 19 cm.

Each radiator was viewed by three x-ray detectors, each consisting of a 19 cm $`\times `$ 19 cm $`\times `$ 5 mm thick NaI(Tl) crystal hermetically sealed between a 6.4 mm thick glass optical window on one flat face and a 0.75 mm thick aluminum entrance window on the other face. An ultraviolet transmitting Lucite lightguide was coupled to the glass window using optical grease, reducing the aperture to a 13 cm diameter circle. An Electron Tubes 9390KB 130 mm photomultiplier tube with a standard bialkali photocathode was mated to the lightguide with optical grease. The lightguide was wrapped in aluminum foil and the whole assembly wrapped with black tape to make it light tight.

Six identical modules were constructed, with each module containing a radiator and the three NaI(Tl) detector assemblies, one on each side of the radiator and one above the radiator outside of and parallel to the beam (Fig. 2). The modules were positioned one behind the other along the beam and aligned such that the particle beam travelled down the center of the modules. Only x-rays scattered at large angles away from the beam were then detected.

For the first accelerator run, the 18 PMT signals were fed into CAMAC-based, 11-bit CAEN C205A charge ADCs and read out with the CERN CMS H2A DAQ computer. For the second accelerator run, the signals from the PMTs were fed into custom-built 8-channel front end modules which contained a charge integrator, peak detect and hold circuit, and gain-adjustable amplifier. Analog-to-digital conversion was performed with a 64-channel, 12-bit National Instruments PCI-6071E DAQ board running in a PC and controlled under LabVIEW 6i.

The instrument was exposed to high energy electrons at the CERN SPS H2A test beam site in August/September 1999 and again in August/September 2001. Beam energies ranged from 7 to 150 GeV, covering the range of Lorentz factors $`\gamma =1.4\times 10^42.9\times 10^5`$. A set of scintillators in the beam upstream of the TRD provided event triggering. The trigger rate was kept to about 1 kHz to avoid deadtime in the DAQ system. Beam definition scintillators in front of and behind the TRD (S1 and S2 in Fig. 2) flagged events for which the electrons showered within the radiator stacks. A Pb shower counter (S3) was placed downstream of S2 to flag pions present as a contaminant in the higher energy beams.

Energy calibration runs were performed both immediately before and after the beam runs using radioactive 133Ba (81, 303 and 356 keV) and 137Cs (662 keV) x-ray sources. In order to account for bremsstrahlung and other background produced by the electrons in passing through the radiators and upstream material, a background run was performed for each radiator configuration in which the radiators were replaced by solid blocks with the same material and thickness (in g/cm2) as the radiators.

## 3 Results

For each material, a background run and a foil run were made for each electron energy used. Figure 3 shows spectra obtained for aluminum foil and background runs at 150 GeV. Several conclusions can be drawn immediately: first, Compton scattered transition radiation is being detected away from the path of the incident electron, at levels well above that of the background; second, the detected TR x-ray spectrum peaks near 100 keV, with some x-rays detected at energies $`>200`$ keV; and third, the measured intensity increases as the particle moves downstream through the set of radiator/detector modules.

Calculated spectra are produced by a Monte Carlo routine based on the description in , in which the differential intensity per unit solid angle per unit frequency is expressed in terms of the coherent sum of the complex amplitudes from the individual interfaces . The effect of the metal foils is included by incorporating an effective absorption cross section (i.e., the imaginary part of the wave vector) that depends on the foil conductivity. Individual x-rays in the range $`21000`$ keV are produced at random locations along the trajectory and followed through the geometry of the radiator stacks and detector modules taking into account the effects of photoelectric absorption, Compton scattering, fluorescence, and escape in both the detectors and radiators, and photoelectron statistics and electronic resolution in the scintillators, photomultipliers, and electronic readout. Examples of the pure TR spectra, showing the characteristic interference pattern and a maximum in the predicted spectrum near 200 keV, are shown in . The calculated TR plus background spectra are shown here: The triangles in Fig. 3 show the result of convolving the measured background in Module 1 with the expected TR signal; the diamonds show the result in Module 5. The predicted signals agree well with the measured TR-plus-background spectra.

Figure 4 shows the total number of photons detected per NaI detector summed over the x-ray energy range $`35500`$ keV as a function of electron energy. The points show the measured data; the curves show the results calculated as described above. The observed saturation Lorentz factors for the thick Teflon and aluminum honeycomb are $`10^5`$, as expected from the calculated values in Table 1. The calculations reproduce the differences in detected yield between different radiator materials, the effects due to different thicknesses of the same material, and the dependence on electron energy to an accuracy of $`20\%`$.

As shown in Fig. 3, the number of x-rays measured in each detector depends on the position of the module. Most of the TR x-rays produced in the beginning of the first module will pass through that module without interacting ($`73\%`$ probability in an Al honeycomb radiator at 100 keV). But as they encounter the radiator material in successive modules, they can Compton scatter and be absorbed in a downstream detector. For the Al honeycomb, there is an 80% probability that a 100 keV x-ray created at the front of the first module will Compton scatter in the radiators before leaving the last module. This feedthrough effect enhances the number of photons detected in modules downstream and can be used to advantage in designing a practical detector . For 150 GeV electrons and Al radiators over the x-ray energy range $`35500`$ keV, the ratios of x-rays detected in the downstream modules compared to Module 1 are $`1.97\pm 0.12`$, $`3.34\pm 0.19`$, $`3.63\pm 0.21`$, $`4.00\pm 0.22`$, and $`4.31\pm 0.25`$ for Modules 2–6, respectively. The corresponding predicted ratios are $`2.31\pm 0.04`$, $`3.30\pm 0.05`$, $`3.95\pm 0.05`$, $`4.24\pm 0.04`$, and $`3.74\pm 0.13`$.

## 4 Conclusion

A new Compton Scatter Transition Radiation Detector capable of measuring the x-rays produced by particles with Lorentz factors near $`\gamma =10^5`$ has been built and successfully tested. Compton scattered TR x-rays were detected outside of the particle beam using relatively thick NaI scintillator detectors, effectively isolating the TR signal from the ionization signal. For thick Teflon and Al honeycomb radiators, the detected x-ray spectrum peaks near 100 keV with some x-rays of energy $`>200`$ keV detected, and saturation Lorentz factors near $`10^5`$ were achieved. The detected yields for most radiator configurations agree with detailed simulations, including the enhancement expected from metal foils. The measurements tend to saturate at slightly higher particle energies than predicted, and the measured number of photons at saturation are $`1020\%`$ higher than predicted for the thin Mylar and Teflon radiators. These discrepancies are presumably due to nonuniformities known to be present in the radiator material. Likewise, the measured peak of the TR plus background spectrum from the aluminum honeycomb radiator occurs approximately 30 keV higher than expected. This discrepancy presumably reflects the approximations inherent in treating an irregular honeycomb structure as a periodic stack of foils. In the case of hard x-ray energies, photon feedthrough enhances the signal in the downstream detectors, again as predicted by the simulations (within $`15\%`$).

## 5 Acknowledgments

This work was supported by NASA grant NAG5-5177 and NASA/Louisiana Board of Regents grant NASA/LEQSF-IMP-02. The authors wish to thank A. Aranas, S. Apewokin, T. Brown and O. Shertukde for their many hours of work during the construction and testing of the radiators and detectors; the staff at CERN for their superb cooperation and assistance, particularly D. Lazic and G. Bencze; and J. Anderson, J. Marsh and C. Welch for assisting with the data analysis.


text:
# On the irreducibility of the two variable zeta-function for curves over finite fields

## 1 Introduction

Let $`X`$ be a proper, smooth, geometrically connected curve of genus $`g`$ over the finite field $`𝔽_q`$. The zeta-function of $`X/𝔽_q`$ can be written as a power series

$$Z(t)=\underset{D}{}\frac{q^{h^0(D)}1}{q1}t^{deg(D)}.$$

Here the sum is over $`𝔽_q`$-rational divisor classes of $`X`$ and
$`h^0`$($`D`$):= dim$`_{𝔽_q}`$ $`H^0(X,𝒪(D))`$. Writing $`b_{nk}`$ for the number of divisor classes of degree $`n`$ and with $`h^0(D)`$=$`k`$ this becomes

$$Z(t)=\underset{n0}{}\underset{k1}{}b_{nk}\frac{q^k1}{q1}t^n.$$

In \[P\] R. Pellikaan observed that the classical proof of rationality and the functional equation for $`Z(t)`$ go through when $`q`$ is treated as a variable in this expression. He thus introduced the following power series in \[P\], Def. 3.1:

$$Z(t,u):=\underset{n0}{}\underset{k1}{}b_{nk}\frac{u^k1}{u1}t^n.$$

This is called the two variable zeta-function of the curve. We will denote by $`h`$ the class-number of $`X/𝔽_q`$, i.e. $`h=`$$`|`$Pic0($`X`$)$`|`$. Then Pellikaan proved:

###### Theorem 1

We have

$$Z(t,u)=(1t)^1(1ut)^1P(t,u)\text{ with }P[t,u].$$

Furthermore:
1) degt $`P`$=$`2g`$, degu $`P`$=$`g`$.
2) In the expansion $`P(t,u)=_{i=0}^{2g}P_i(u)t^i`$ one has $`P_0(u)=1`$,
degu $`P_i(u)i/2+1`$ and $`P_{2gi}(u)=u^{gi}P_i(u)`$ for $`0i2g`$.
3) $`P(1,u)=h`$.

Here degu and degt denote the degree of a polynomial in the indicated variable. The above results are all taken from \[P\], Prop. 3.5 and we copied only those needed later on. Note that the statement degu $`P_i(u)i/2`$ in \[loc. cit.\] is a misprint. Indeed, we will see below that one always has degu $`P_1(u)=1`$ (unless $`g=0`$). As expected we have $`P(t,u)=1`$ in case $`g=0`$.

In \[GS\] G. van der Geer and R. Schoof used analogies from Arakelov-theory to define a two variable zeta-function for number fields along the above lines. As the number field case will serve only as a motivation in this note we refer to the original sources \[GS\] and \[LR\] for definitions and to \[D\] for a comparison between them. Suffice it to say that in \[LR\], section 8 we find an entire function $`\xi _{}(w,s)`$ of two complex variables which for $`w=1`$ equals Riemann’s $`\xi `$-function. In particular, the zeroes of $`\xi _{}(1,s)`$ are precisely the non-trivial zeros of the Riemann zeta-function. One is thus led to study the zero-locus of $`\xi _{}(w,s)`$. J. Lagarias and E. Rains ask whether it might be the closure of a single irreducible complex-analytic variety of multiplicity one. The corresponding question in the geometric case seems to be whether the zero-locus of $`P(t,u)`$ is irreducible.
This is indeed the case:

###### Theorem 2

In the above situation, $`P(t,u)`$ is irreducible in $`(u)[t]`$.

As an illustration we discuss the cases $`g=1`$ and $`g=2`$:
For $`g=1`$ setting $`N:=|X(𝔽_q)|`$ we have

$$P(t,u)=1+(N1u)t+ut^2,$$

c.f. \[P\], example 3.4. This polynomial is reducible in $`(u)[t]`$ if and only if $`N=0`$ in which case we have $`P(t,u)=(1t)(1ut)`$. But it is well-known that a curve of genus one over a finite field always has a rational point, i.e. $`N0`$.
In case $`g=2`$ let the usual zeta-function of $`X`$ be

$$Z(t)=(1t)^1(1qt)^1L(t)$$

with

$$L(t)=1+at+bt^2+qat^3+q^2t^4$$

for certain $`a,b`$. As $`X`$ is hyperelliptic, Prop. 4.3. of \[P\] can be used to compute

$$P(t,u)=1+((a+q)u)t+((q(q1)+aq+b)(a+q1)u)t^2+((a+q)u)ut^3+u^2t^4.$$

This will not be used in the sequel and we omit the proof.
In order not to lead intuition astray we point out that in general $`Z(t,u)`$ is not determined by $`Z(t)`$, see \[P\], example 4.4.
After a lengthy computation with discriminants one sees that a neccessary condition for this $`P(t,u)`$ to be reducible is

$$b+a(q+1)+(q^2+1)=0.$$

But this expression equals $`L(1)=h0`$ !
The fact that $`h0`$ enters in the general proof of theorem 2 precisely through the condition $`\beta 0`$ of lemma 3 below. Note, however, that condition 2) of this lemma cannot be dropped. So one needs one more result on $`P(t,u)`$, contained in Prop. 4, which follows from Clifford’s theorem.

I would like to thank C. Deninger for posing the above problem and for useful discussions on the topic.

## 2 Proof

We will use the following criterion for irreducibility:

###### Lemma 3

Let $`k`$ be a field, $`Fk[u,t]`$ and assume:
1) $`F`$ is monic in $`t`$.
2) the leading coefficient of $`F`$ as a polynomial in $`u`$ is irreducible in $`k[t]`$.
3) there are $`\alpha ,\beta k`$, $`\beta 0`$ with $`F(u,\alpha )=\beta `$.

Then $`F`$ is irreducible in $`k(u)[t]`$.

This lemma will be applied to

$$F=\stackrel{~}{P}(t,u):=t^{2g}P(t^1,u)[u,t].$$

Note that the irreducibility of $`\stackrel{~}{P}`$ in $`(u)[t]`$ will imply the irreducibility of $`P`$ because $`P(0,u)=P_0(u)=10`$ by theorem 1. The advantage of $`\stackrel{~}{P}`$ is that it is monic in $`t`$ and so satisfies condition 1) of lemma 3. Also 3) is satisfied (with $`\alpha =1`$, $`\beta =h`$) according to theorem 1, 3).

Proof (of lemma 3) Assume to the contrary that $`F=fg`$ in $`k(u)[t]`$ with $`f`$ and $`g`$ of positive degree and monic. One knows, c.f. for example \[E\], Prop 4.11, that the coefficient of $`f`$ and $`g`$ are integral over $`k[u]`$ and as $`k[u]`$ is integrally closed we have $`f,gk[u,t]`$. So we can consider the decomposition $`F=fg`$ as polynomials in $`u`$ and infer from 2) that the leading coefficient of $`f`$ as a polynomial in $`u`$ lies in $`k[t]^{}=k^{}`$ (upon exchanging $`f`$ and $`g`$ if neccessary). In particular $`n:=`$degu $`f(u,t)=`$degu $`f(u,\alpha )`$. Substituting 3) gives $`\beta =f(u,\alpha )g(u,\alpha )`$ in $`k[u]`$. As $`\beta 0`$ we get $`n=0`$, i.e. $`f`$ is constant in $`u`$ hence $`fk^{}`$, contradiction.

$`\mathrm{}`$

We are left with verifying condition 2) of lemma 3 for the given $`\stackrel{~}{P}`$, i.e. the leading coefficient of $`\stackrel{~}{P}`$ as a polynomial in $`u`$ is irreducible in $`k[t]`$. We will in fact determine this coefficient:

###### Proposition 4

For $`g1`$: $`\stackrel{~}{P}(t,u)=(1t)u^g+O(u^{g1})`$.

Proof We already know degu $`\stackrel{~}{P}=g`$. Also the assertion is clear for $`g=1`$ from the formula for $`P(t,u)`$ recalled in the introduction. We assume $`g2`$. Looking at

$$\stackrel{~}{P}(t,u)=t^{2g}+P_1(u)t^{2g1}+\mathrm{}+P_g(u)t^g+uP_{g1}(u)t^{g1}+\mathrm{}+u^gt^0$$

and using the bound degu $`P_i(u)i/2+1`$ we see that $`u^g`$ can only occur in the last three terms: $`u^{g2}P_2(u)t^2+u^{g1}P_1(u)t+u^g`$. So the proof is completed by the following result on $`P_1`$ and $`P_2`$.

$`\mathrm{}`$

###### Proposition 5

For $`g1`$:
1) degu $`P_1(u)=1`$ and the leading coefficient is $`1`$.
2) degu $`P_2(u)1`$.

Proof This is again clear for $`g=1`$. We assume $`g2`$ and write $`P_i(u)=_k\alpha _{ik}u^k`$, $`\alpha _{ik}`$. As we already know degu $`P_1(u)1`$ and degu $`P_2(u)2`$ we need to show $`\alpha _{11}=1`$ and $`\alpha _{22}=0`$. Recalling the notation $`b_{nk}`$ from the introduction we have the following

Claim: $`b_{12}=\alpha _{00}+\alpha _{11}`$ and $`b_{23}=\alpha _{00}+\alpha _{11}+\alpha _{22}`$.

Granting this we observe that Clifford’s thorem, c.f. \[H\], IV, thm 5.4, gives $`b_{12}=b_{23}=0`$. Recalling also $`\alpha _{00}=1`$, because $`P_0(u)=1`$, and substituting gives indeed $`\alpha _{11}=1`$ and $`\alpha _{22}=0`$.

$`\mathrm{}`$

To prove the above claim we write the rational expression for $`Z(t,u)`$ in theorem 1 in terms of coefficients:

$$Z(t,u)=\underset{n0}{}\underset{k1}{}b_{nk}\frac{u^k1}{u1}t^n=(\underset{i0}{}t^i)(\underset{j0}{}(ut)^j)(\underset{l0}{}\underset{k0}{}\alpha _{lk}u^kt^l).$$

This gives

###### Lemma 6

1) for $`\nu `$,$`\alpha 0`$ : $`_{k\alpha +1}b_{\nu k}=_{\mu ,i0,\mu +i\nu }\alpha _{i,\alpha \mu }`$
2) for $`\nu 0`$, $`\mu 1`$: $`b_{\nu \mu }=_{i=0}^\nu (\alpha _{i,\mu \nu 1+i}\alpha _{i\mu })`$

We omit the details of this straightforward computation except to say that for 1) one uses $`\frac{u^k1}{u1}=1+\mathrm{}+u^{k1}`$ and 2) follows from 1) by a telescope-summation. In the formulation of the lemma it is understood that $`\alpha _{nk}=0`$ whenever $`k<0`$. We get from 2):

$$b_{23}=\alpha _{00}\alpha _{03}+\alpha _{11}\alpha _{13}+\alpha _{22}\alpha _{23}$$

and

$$b_{12}=\alpha _{00}\alpha _{02}+\alpha _{11}\alpha _{12}.$$

But we know $`\alpha _{02}=\alpha _{03}=\alpha _{12}=\alpha _{13}=\alpha _{23}=0`$ because degu Pi(u)$`i/2+1`$.
This concludes the proof of the claim, hence of theorem 2.

Mathematisches Institut der WWU Münster
Einsteinstr. 62
48149 Münster
Germany
e-mail: naumannn@uni-muenster.de