text:
# Long time relaxation of interacting electrons in the regime of hopping conduction

## I Introduction

The study of electron-electron interactions in the localized regime was initiated by PollakPollak70  and Srinivasansri . Later on Efros and Shklovskii(ES) ES  showed that the single particle density of states (DS) tends to zero as the energy tends to the Fermi energy. This phenomenon, called the Coulomb gap, is due to the long range part of the Coulomb interaction which, in some sense, remains non-screened. In fact, the Coulomb gap results from the Coulomb law and from the discrete nature of the electron charge. In their first works ES claimed that the DS in the Coulomb gap has a universal form, depending only on electron charge $`e`$ and dielectric constant $`\kappa `$. Then DS $`|ϵ|^{D1}\kappa ^D/e^{2D}`$, where $`D`$, is the spatial dimension and $`ϵ`$, the single-electron energy whose reference point is the chemical potential. This expression for the DS is the only combination of the energy and the electron charge which has a proper dimensionality. It has been shown later on pik  that in 2D-case the above universality is exact only for strong disorder. In the 3D-case the question was never studied in detail, however deviations from the quadratic law have been reportedMobius92 .

Simple quantitative arguments which assume that the single-particle excitations are responsible for the variable range hopping (VRH) lead to the so-called Efros-Shklovskii (ES) lawES , which has been observed experimentally in many materials

$$\sigma _c(\gamma e^2/T)\mathrm{exp}\left((T_0/T)^{1/2}\right),$$
(1)
where $`T`$ is the temperature and $`T_0=\beta _0e^2/\kappa a`$, $`a`$ is a localization length of electrons, $`\kappa `$ is an effective dielectric constant of the media above and below 2D-gas. A self-consistent type of percolation approach (See references inlev ) gives $`\beta _0=6.5`$. The hopping length is given by $`R_C(a/4)(T_0/T)^{1/2}`$.

Recentlyourprl , using the so-called lattice model, the ES law has been checked in detail by computer simulation. The Hamiltonian of this model is formulated on a square lattice and has a form

$$H=\underset{i}{}\varphi _in_i+\frac{1}{2}\underset{ij}{}\frac{\left(n_i\nu \right)\left(n_j\nu \right)}{r_{ij}},$$
(2)
where $`n_i=0,1`$ are occupation numbers. The quenched random site energies $`\varphi _i`$ are distributed uniformly within the interval $`[1,1]`$, and the average occupation number $`\nu `$ , is taken to be 1/2. The magnitude of the quenched disorder is enough so as to provide the universal Coulomb gap at all energies which are important within the temperature range under studypik . In what follows the lattice constant is taken to be the unit of length. The nearest neighbor Coulomb energy which is in this case equal to the amplitude of the disorder is considered both as the energy and the temperature unit.

Simulations of the conductivity in the lattice modelourprl  confirm the ES law in all details, i.e., the pre-exponential factor, T-dependence and $`a`$-dependence. It has also been shown that simultaneous transitions of two electrons do not play any role. Arguments have also been given that any many-electron excitations are not important. Therefore simultaneous transitions of electrons were not included in the simulation below.

Glassy properties, due to both randomness and the long range Coulomb interaction, are another interesting manifestation of electron-electron interactions in such a system. Davies, Lee, and Rice daviesprl ,daviesprb  were the first to raise this issue. They coined the phrase “Electron Glass” which is still used, sometimes also referred to as the “Coulomb Glass”. Both these terms stress the relation of the above electron system to a spin glass system.

We believe that, in the same way as in real glasses, the glassy properties in this electron system are due to those states, which have very close total energies but substantially different sets of the occupation numbers $`n_i`$. Such states have been first observed by Baranovskii et al.Baranovskii  during the first computer simulation of this system and have been called “pseudoground states” (PS’s). These authors found that the PS’s have the same universal Coulomb gap and have concluded that the existence of PS’s is not important for the VRH conductivity because transitions between them are very slow. They attributed these slow transitions between the PS’s as well as their difference from the ground state to be a result of a certain amount of many-electron transitions.

Experiments on the glassy properties of such systems, conducted around 25 years later, have confirmed this conclusion, but they have also shown that Baranovskii et al. missed an important feature. Since transitions between PS’s take a long time and conductivity of these states is not exactly the same, they can serve as a basis for memory effects.

Experiments, started by the group of Ovadyahuzvi1 ; zvi2  in 1993, definitely show the relation of this system to the ordinary glasses. In these experiments the difference in the conductivities of the PS’s is not larger than 10-12%. Similar phenomena was observed by the group of Goldmangol  on ultra-thin films of metals near the superconductor-insulator transition. Slow relaxation has been demonstrated by Don Monroe at al.dm  in compensated GaAs.

The properties of the PS’s have been studied recently, mostly by computational methods kogan ; perez98 ; menashe . Perez-Garrido et al.perez98  argued that transitions between PS take a huge time, substantially larger than a time available in any experiment. Menashe at al.menashe  have proposed a different method in order as to study the Coulomb Glass. By completely ignoring the tunneling term in the transition probability, still keeping the activation probability for the electrons, they performed a thermodynamic Monte-Carlo simulation, which transforms a non-ergodic system into an ergodic one. Thereby they were able to obtain all the thermodynamic properties of the ergodic system, since the electrons were able to move across the system in a single transition, leading to effective mixing of all PS’s.

The above method only permits calculation of the thermodynamic values, since the number of the Monte-Carlo steps can not be related to a physical time. Still, these authors claimed that the metal-insulator transition in the Coulomb glass may coincide with the glassy transition, that occurs due to increase of the localization length. Their method menashe  is partially employed in this paper.

The goal of this paper is to understand the origin of the long-time relaxation of conductivity, observed in the experimental papers cited above. To do so we performed Monte-Carlo simulations employing: the lattice model Eq.(2) and the random site model described later on, for two dimensional systems. In our simulation we perturb the system by adding some extra electrons and then trace the relaxation of energy and conductivity with time. One should clearly understand that it is impossible to observe the long-time relaxation by a direct simulation. The general reason for this is due to the long range of the Coulomb interaction and to the fact that in a real sample many transitions take place simultaneously while a computer processor is limited to performing them one at a time as well as updating the resulting changes in the site energies one at a time. A few additional processors cannot help much. Therefore, even using the most sophisticated program, we may simulate no more than 40 $`\mu `$sec of the real time in a system of the size of $`100\times 100`$ lattice sites.

We have shown that during such available times we can reach an apparent saturation of conductivity and energy and that saturation of the DS occurs much earlier. We interpret this “saturation” as a saturation within one PS. Thus, it is rather the transition from the fast relaxation within one PS to a very slow relaxation to PS’s with lower energy, that we can not observe.

Simultaneously with the conductivity we studied the relaxation of energy. We devised an analytical theory for the energy relaxation which fits the simulation data fairly well and serves as a reference point for our understanding of the short time relaxation.

To study the long time relaxation due to transitions between the different PS’s we came up with a new idea. We create different PS’s by relaxation from states with different initial distributions of electrons. Then we study the difference between the saturated values of the conductivities of the different PS’s if these are different, one should expect that a long time relaxation exists and the total change of the conductivity should be of the order of this difference. For the lattice model we obtain no such effects. The conductivities of different PS are the same in the limits of our accuracy (about 1-2%).

Since PS’s with different energies were observed by Menashe et al. as well as in other papers, we arrive at the conclusion that the effect of long-time relaxation of the conductivity is absent due to the universality of the Coulomb gap for the lattice model, in the temperature range which we consider.

We also performed simulations for the random site model. In this model all disorder comes from the random position of sites. Our results for these simulations exhibit a difference in the conductivities of the PS’s whose value is sufficiently large so as to explain the experimental results of Ovadyahu’s group.

The difference between these two models in the 3-dimensional case has been discussed in Ref.discussion . Note that in this case a glassy transition at non-zero temperature has been claimedGrannan . Xue and Lee Xue  performed Monte-Carlo simulations employing the 2-dimensional random site model in which the disorder comes about only through the random position of sites. These authors found evidence for glassy behavior at low temperature but they claimed the absence of the glassy transition at non-zero temperature. The same result was obtained for the Ising modelbinder .

The paper is organized as follows. In the next section we study the relaxation of the energy and the conductivity after an extra amount of electrons were added to the system in the framework of the lattice model. In Sec.III, employing the lattice model, we study the difference in the conductivities of the PS’s, obtained through different relaxation methods. In Sec.IV the results for the random site model are presented and discussed.

## II  Relaxation after addition of extra electrons

In this section we study the relaxation of the energy and conductivity of the system after it has been initially perturbed by the addition of a few percents of extra electrons. In principle this is the same procedure employed in the experiments of Ovadyahu’s group but the time in which we are able to observe the system is very short. We use the lattice model with Hamiltonian given by Eq.(2) and initial filling factor $`\nu =1/2`$. The extra electrons are added randomly on empty sites and the background is adjusted such that the system remains neutral.

### II.1 Analytical theory of energy relaxation

To describe the physics of the short time energy relaxation we devise the following analytical theory. We start by considering a division of the plain with the 2D electrons into regions with a linear size R. If $`\delta n`$ is an average density of additional electrons, the charge $`Q`$ of each region is of the order of $`Qe\sqrt{\delta nR^2}`$, because the average charge is compensated by the background. The extra energy due to electron-electron interaction per a region containing excess charge $`Q`$ is $`Q^2/R`$. The number of regions per area is $`L^2/R^2`$, where $`L`$ is the length of the total system which is square shaped. Thus the extra energy per area due to all regions of size $`R`$ is $`=e^2\delta n/\kappa R`$. The regions with the smallest $`R`$ give the largest contribution. However, due to relaxation they become neutral faster. So, the main contribution at a time $`t`$ comes from the regions in which relaxation has not yet ended. In the 2D case the relaxation goes with the velocity of the order of the conductivity $`\sigma _0`$ of the systemdyakon . Thus, at a time $`t`$ only the regions with $`R>\sigma _0t`$ have an excessive charge. Finally, the energy per area decreases with time as $`e^2\delta n/\kappa \sigma _0t`$. Analyzing results for the auto-correlation functions and density-density correlation functions it is not difficult to find the numerical coefficient in this expression. Let $`n^{}(𝐫)`$ stand for the density of additional electrons with their homogeneous background so that the average value $`<n^{}>`$ is zero. The linearized equations have the form

$$𝐣^{}=\sigma _0𝐄$$
(3)
and

$$e\frac{n^{}}{t}+\sigma _0𝐄=0.$$
(4)
One can see that the diffusion current, omitted here, is more important on the earlier stages of relaxation than the Ohmic current.

In the 3D case the field can be eliminated from these equations resulting in an equation for $`n^{}`$ only. In the 2D case the situation is different and the problem can be solved using the Fourier transformation,

$$n^{}=\frac{1}{L}\underset{𝐪}{}n_𝐪e^{i\mathrm{𝐪𝐫}},\varphi =\frac{1}{L}\underset{𝐪}{}\varphi _𝐪(z)e^{i\mathrm{𝐪𝐫}}.$$
(5)
Here $`𝐫`$ and $`𝐪`$ are two dimensional vectors in the plane of the electrons, $`z`$ is the coordinate in the perpendicular direction and $`\varphi `$ is the scalar potential. The equation for the potential has the form

$$\frac{d^2\varphi _𝐪}{dz^2}q^2\varphi _𝐪=\frac{4\pi e}{\kappa }n_𝐪\delta (z).$$
(6)
It can be solved in the regions $`z>0`$ and $`z<0`$. Under the condition of continuity of $`\varphi _𝐪`$ at $`z=0`$ and zero at $`z\pm \mathrm{}`$ the solution is

$$\varphi _𝐪=A_𝐪e^{q|z|}.$$
(7)
The integration of Eq.(6) over an infinitesimally small interval around $`z=0`$ leads to the boundary condition

$$\frac{d\varphi _𝐪}{dz}|_{z=+0}\frac{d\varphi _𝐪}{dz}|_{z=0}=\frac{4\pi e}{\kappa }n_𝐪.$$
(8)
This condition gives

$$A_𝐪=\frac{2\pi e}{\kappa q}n_𝐪.$$
(9)
Eqs.(7) and (9) give

$$𝐄_𝐪=i𝐪\frac{2\pi e}{\kappa q}n_𝐪,(𝐄)_𝐪=i\mathrm{𝐪𝐄}_𝐪=\frac{2\pi eq}{\kappa }n_𝐪.$$
(10)
Now the Fourier transformation of Eq.(4) is

$$\frac{dn_𝐪}{dt}+\frac{n_𝐪}{\tau _𝐪}=0.$$
(11)
where

$$\frac{1}{\tau _q}=\frac{2\pi q}{\kappa }\sigma _0.$$
(12)
This leads to

$$n_𝐪(t)=n_𝐪(0)e^{t/\tau _q}.$$
(13)
The energy of the fluctuations is

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \varphi (𝐫,z)en(𝐫)\delta (z)d^2r𝑑z}={\displaystyle \frac{\pi e^2}{\kappa }}{\displaystyle \frac{|n_𝐪(0)|^2}{q}e^{2t/\tau _q}\frac{d^2𝐪}{(2\pi )^2}}.`$ (14)
If the initial state corresponds to randomly distributed electrons,

$$n^{}(𝐫,0)n^{}(𝐫^{},0)=\overline{\mathrm{\Delta }n}\delta (𝐫𝐫^{}),$$
(15)
then

$`|n_𝐪(0)|^2`$ $`=`$ $`{\displaystyle \frac{1}{L^2}}{\displaystyle e^{i\mathrm{𝐪𝐫}}n^{}(𝐫,0)𝑑𝐫e^{i\mathrm{𝐪𝐫}^{}}n^{}(𝐫^{},0)𝑑𝐫^{}}`$ (16)
$`=`$ $`{\displaystyle \frac{\overline{\mathrm{\Delta }n}}{L^2}}{\displaystyle e^{i𝐪(𝐫𝐫^{})}\delta (𝐫𝐫^{})𝑑𝐫𝑑𝐫^{}}=\overline{\mathrm{\Delta }n}.`$
The substitution of Eq.(16) in Eq.(14) gives

$``$ $`=`$ $`{\displaystyle \frac{\pi e^2\overline{\mathrm{\Delta }n}}{\kappa }}{\displaystyle \frac{1}{q}e^{2t/\tau _q}\frac{d^2𝐪}{(2\pi )^2}}={\displaystyle \frac{e^2\overline{\mathrm{\Delta }n}}{8\pi \sigma _0t}}`$ (17)
This result gives the value of the numerical coefficient for the above estimate. One should understand that Eq.(17) is valid at $`t<L/\sigma _0`$. At such times the energy relaxation has finished due to the finite system size. It is easy to take the size effect into account qualitatively. The result is that energy as a function of $`t^1`$ saturates at $`tL/\sigma _0`$. It is important also that the system is assumed to be ergodic, i.e., it does not contain different PS’s with slow relaxation between them.

### II.2 Simulation results for short time relaxation of energy and conductivity

All the time-dependent simulation results in this paper are obtained by our modification of the kinetic Monte- Carlo method that is presented in the Appendix. For the lattice model our computer units (CU) are as follows: the length unit is taken to be the lattice constant $`a_0`$, the units of energy and temperature are given by the Coulomb interaction at the lattice constant $`e^2/a_0`$, the time unit is chosen as the reciprocal transition rate due to phonons $`\gamma ^1`$ see Eq.(1), and the unit of the 2-dimensional conductivity is given by $`a_0\gamma `$. The numerical values for the time given below, are calculated using the assumption that $`\gamma =10^{12}s^1`$.

The simulation temperatures were chosen to be $`T=0.1,0.2`$ which corresponds to the region of hopping conduction in the case were the localization radius $`a`$ is of the order of 1 CUourprl . Decreasing the value of $`a`$ we may switch from the VRH to the nearest-neighbor hopping regime. The size of the system was taken to be $`100\times 100`$ lattice sites.

The simulation was performed using the following steps: In the first stage the system was brought to “thermal equilibrium”, inside a single PS, using the thermodynamic method employed in the work of Menashe et al.menashe . In the second stage, which is the reference point for the time from which the relaxation of the DS, energy and conductivity were studied, the electron concentration was slightly changed by adding (or removing) electrons randomly. The time dependences were averaged over $`10^3`$ different sets of randomly distributed extra electrons (holes). Due to electron-hole symmetry of the Hamiltonian Eq.(2) at $`\nu =1/2`$ the relaxation of extra electrons and extra holes is the same. The data below describes the addition of electrons. It should be noted that the insertion of electrons into the system can be simulated in many ways ,e.g., adding electrons in a confined spatial region simulating the way electrons might enter the sample experimentally, or adding electrons to the most probable states. It is clear however that since the unperturbed system is saturated within one PS the difference between different electron insertion scenarios should only be notable within very short time scales i.e. of the order of our simulation times. Thus the choice of how to simulate the perturbation of the system by adding charge does not in any way effect the long time relaxation of the system which is due to transitions between different PS’s.

The energy relaxation is shown in Fig.1 and Fig.2. The reference point for the total energy $`E_{ex}`$ is the saturated value for the longest time measured. At time $`t=0`$ the filling factor $`\nu `$ has been increased from 0.5 to 0.52, such that the total number of extra electrons equals $`N_{ex}=0.02L^2/2`$. Fig.1 shows the regime of the VRH with $`a=1`$ at two different temperatures. In Fig.2 we compare relaxations of energy in the VRH regime at $`a=1`$ with the relaxation in the nearest neighbor hopping at $`a=0.2`$. The values of the VRH conductivity at $`a=1`$ are $`\sigma _0=0.0048,0.021`$ in CU for $`T=0.1,0.2`$ respectively. These values of conductivity are also obtained as the saturated value at the longest time measured. Since Eq.(17) shows that at large $`t`$, the energy $`E`$ is a function of $`\sigma _0t`$, we use this product as a reciprocal argument in Figs.2,3. The straight line in the inset of Fig.1 shows the slope as given by Eq.(17). The inset shows the final stage of the relaxation where the curves at different temperatures coincide and obey the time dependence given by Eq.(17). The saturation at even larger $`t`$ is due to a size effect. Note that the values of the conductivity $`\sigma _0`$ differ by almost five times for the two curves presented. The observed saturation of energy relaxation is connected to the finite size effect as discussed in a previous section.

The relaxation of the charge fluctuations given by Eq.(17) not only has the universal form for the different values of conductivities, but it is also independent on the conduction mechanism. In the Fig.2 we present the relaxation curves for the system for VRH conductivity and for nearest-neighbor hopping where the localization radius of the electrons is much smaller than the distance between neighboring sites. The time region where the two curves coincide within our accuracy is even broader despite the fact that the values of the conductivity differs by a factor of $`3000`$. The main source of the errors (see Appendix) arises from calculating the values of the conductivity $`\sigma _0`$ rather than measuring the total excessive energy of the system.

Next we present the results for relaxation of the conductivity after adding extra electrons to the system. The procedure for adding the electrons is the same as stated above. The time dependence of the conductivity of the system is shown on Fig.3 for two different temperatures. These were obtained in the same simulation as the data for the relaxation of the energy in Fig. 1. One can see that for both temperatures the conductivity decreases with time finally reaching a kind of saturation, which may also be interpreted as a transition to a substantially slower rate of relaxation. We consider this ’saturation’ as the end of the short time relaxation in our finite system. The value of the conductivity at the largest time we considered, is denoted by $`\sigma _0`$. One can see that as a function of $`1/\sigma _0t`$ the results for both temperatures coincide even in an even wider time region than for the total energy of the system. The reason that the conductivity decreases is that, due to the non-equilibrium the extra electrons occupy states above the Fermi level therefore providing higher current than in thermal equilibrium. Note that the saturation of conductivity occurs at the same value of $`\sigma _0t`$, that is $`4`$ as the saturation of energy. It is reasonable to think that the two processes are connected and the saturation of the conductivity is also a size effect.

The time corresponding to the saturation $`t=4/\sigma _0`$ is very short. For T=0.1 and $`a=1`$ it is $`0.8ns`$. Even if we assume an infinite system and extrapolate the law $`\sigma \sigma _0(\sigma _0t)^1`$ to much larger times, we find that in the microsecond range all relaxation which can be observed will be over.

The saturation of the conductivity averaged over different initial distributions of extra electrons, which we have presented in this section, happens at very short times. It is a Maxwell-type relaxation of the extra charge. We have checked that the relaxation of the Coulomb gap inside the energy interval that is responsible for the VRH occurs even faster than the relaxations of the energy and conductivity. These results do not support the idea that long time relaxation is due to the slow formation of the Coulomb gapyu .

However in the experimental data zvi1 ,zvi2  the long time relaxation happens on time scales of the order of seconds and hours. These time scales are 9-11 orders of magnitude longer than those which we are able to simulate. Thus, if a long time relaxation exists in this system, the change of the averaged value of the conductivity is negligible during the physical time scales we have simulated. This explains the apparent saturation of the conductivity, we observe. If the long time relaxation results from the transitions between different PS’s, the conductivity we have observed should be considered as the conductivity within a single PS. The encouraging result obtained in this section is that the saturation of the conductivity of a single PS can be achieved during the time scale available for our computation. Based on this result another approach to the problem of long time relaxation, is proposed in the next sections.

## III Conductivity of different pseudoground states in the lattice model

The idea behind our new approach is rather simple. We want to compare the conductivity of the system in the different PS’s. If the values of the obtained conductivities are different for different PS’s then the long time relaxation of the conductivity can be attributed to the slow transition between those PS’s. We are unable to measure time scales of the order of such transitions still we can explain its effect and magnitude which are reflected through the difference in the conductivities of different PS’s. In the experiments of the group of Ovadyahu zvi2  such differences are of the order of 10%.

To observe the conductivity of different PS’s, we measured the conductivity of the same sample with the different initial distributions of electrons during the longest time we are able to simulate. The sample is characterized by the total set of random energies $`\varphi _i`$. Starting with different initial distribution of electrons the system relaxes to the different PS’s. If the saturated values of conductivities, measured as in previous section, is different, one should expect that the system will have a long time relaxation due to transitions between different pseudoground states.

The results for the lattice model simulations performed at the lowest temperature are presented in Fig.4. The time evolution of the conductivity averaged over the time of measurement is shown for different initial distributions of electrons for the same sample. As one can see there is no appreciable difference in the values of the conductivities within our accuracy which is about $`12\%`$. We obtained the same result for all higher temperatures even with greater accuracy.

We therefore conclude that there is no apparent difference in in the conductivities of different PS’s for this model down to the lowest temperatures we are able to simulate. We think that the reason for this is the following: It has been shownourprl  that the VRH conductivity of the system is provided by the single-electron excitations. The properties of these excitations is defined by the structure of the Coulomb gap. In the lattice model at large $`A`$ the DS in the Coulomb gap has a universal formpik . It is independent of the properties of the system and is the same for all PS’s. For our temperature range, the case $`A=1`$ can be considered as a large disorder. Thus all PS’s have the same conductivities and no long time relaxation can be observed.

We realize that the lattice model with large disorder can not account for the effect of the long time relaxation of the conductivity which is observed experimentally. Thus, to observe the difference in the conductivities of different PS one should take a system with smaller disorder, where the Coulomb gap is not universal. Unfortunately, it is difficult to find such a regime in the lattice model, due to the Wigner crystallization at low temperatures.

## IV Conductivity of different pseudo ground states in the random site model

In this section we present the results for the random site model. The Hamiltonian of the model has the form

$$H=\frac{1}{2}\underset{ij}{}\frac{\left(n_i1/2\right)\left(n_j1/2\right)}{r_{ij}}.$$
(18)
It differs from Eq(2) in two important respects. The first: it is formulated on sites i,j, which have random positions on the plane. The second: any random energies $`\varphi _i`$, that are not correlated with the interaction, are absent. Thus, the random positions of the sites is the only source of disorder.

We consider the case $`\nu =1/2`$ for which each state of the system has an exact twofold degeneracy: the total energy is invariant with respect to changing all of the occupation numbers $`n_i1n_i`$. It is this symmetry for half filling which allows this model to be mapped on to a spin glass modelGrannan . Probably, this fact is also important for formation of PS’s that consist of fragments of both states.

This model is unusual for electronic systems. For example partially occupied donors located in a plane with a gate electrode above the plane are not described by this model. In this case the above symmetry is absent. However, the neutral system of random donors with a large “negative U” can be described by the Hamiltonian Eq.(18), if 1/2 of them have double occupation with charge -1 and the other half are not occupied with the charge 1.

For the simulation of the random site model we use a similar computational algorithm which is described in the Appendix. Unfortunately, it is more time and memory consuming than the algorithm for the lattice model. In this model the unit of length is $`a_0=n^{1/2}`$, where $`n`$ is the concentration of sites.

Fig.5 shows the temperature dependences of the VRH conductivities both in the random site model and in the lattice model. Unfortunately for the random site model we were unable to check the importance of simultaneous many-electron transitions on the VRH conductivity. One can see a deviation from the ES law in the case of the random site model, that might be a result of deviation of the DS from the universal DS in the Coulomb gap.

Fig.6 confirms this point of view. It shows DS of different PS’s for both lattice and random site models. As in the previous section the different PS’s have been obtained by simulations starting from different initial distribution of electrons but with the same disorder. The latter condition means the same set of $`\varphi _i`$ or the same set of random sites depending on the model. In each case the DS is calculated by two methods. The first is the kinetic Monte-Carlo, which gives the DS in one PS due to the time constraint, and the second is the thermodynamic Monte-Carlo method which gives the DS averaged over all PS’s. To get to the thermodynamic regime one should ignore the tunneling exponent in the transition rate so that transition at any distance becomes possible. However the energy dependence of the rate should be strictly preserved to get the correct ergodic thermodynamic result. Menashe et al.menashe  have shown that this method provides an effective thermalization including fast transitions between the PS’s.

Fig.6 contains two important results: (i) The DS of the random site model strongly differs from the DS of the lattice model and it does not have the standard energy dependence typical for the two-dimensional Coulomb gap, which is known to be very robust in the lattice model for A=1. For the random site model the DS is quadratic rather than linear. We have checked that its curvature is $`T`$-independent. From dimensionality considerations this DS can only be of the order of $`ϵ^2/e^6n^{1/2}`$, which also explains the deviation from the ES law for the VRH conductivity shown in Fig.5. (ii) The relative difference between the DS for different PS is much larger for the random site model than for the lattice model. For the random site model the time fluctuations of DS near the Fermi level are $`10`$ times smaller than the difference between the DS’s of different PS’s. In the lattice model those fluctuations are so close to the difference itself that the difference is not a reliable measure. The thermodynamic DS is between the DS’s obtained for the different PS’s as it should be for an average function.

In fact results (i) and (ii) are connected to each other. At large $`A`$ the universal behavior of the Coulomb gap can be obtained from conditions $`ϵ_jϵ_i1/r_{ij}>0`$ for every empty site $`j`$ and occupied site $`i`$. These conditions are important near the Fermi energy only, for which the total energy becomes irrelevant. That is the reason why all PS have similar DS near the Fermi level. For the random site model these conditions are also necessary but there should be some other restrictions that make the DS smaller. Those restriction should be connected with the average distance between the sites and therefore they are related to the total energy. Therefore it is interesting to study the conductivity of different PS in this model.

The simulation results for the time evolution of the conductivities of different PS’s for the model with the random spatial site distribution are shown in Fig.7 by filled symbols at the same temperature as for the lattice model. In this case the conductivity value of each individual PS saturates within 1% during the time of the simulation. However, the saturated values of the conductivities differ by 12% for the different initial distributions of electrons in the same sample. This is by an order of magnitude greater than the simulation accuracy. The study of the conductivities in the random site model shows that the slow transitions between PS’s may result in the long time relaxation of the conductivity observed in the experiment. An important question now is whether the obtained result is due to the randomness in the site distribution or to the absence of the disorder which is not correlated with the interaction.

In order to answer this question we consider a “compound” model. Namely we have added to the Hamiltonian of the random site model Eq.(18) the first term in Eq.(2) with $`A=1`$. One can see from Fig.7 that the difference in the conductivities of different PS’s disappears at $`A=1`$ within the simulation accuracy. Thus, the randomness of the sites is not important, but the value of $`A`$ is very important.

Another important question is whether the values of the conductivity of the different PS’s are correlated with the energies of these PS’s. The total energy of the system in a given PS is averaged over time and plotted versus the value of the conductivity for the same PS’s in Fig. 8. The data shows that the energy dependence of the conductivity is close to a linear behavior. PS’s with higher average total energy have also a larger value for the conductivity which is reasonable, since in the states with lower energy the electrons are in positions where they are more tightly bound and therefore their conductivity is lower. The same behavior is apparent in the experimental results of Ovadyahu zvi2 . The linear dependence is also reasonable because the energy difference between PS’s is small and conductivity as a function of energy may be expended into a Taylors series, taking the first term. Similar to Anderson et al.anderson  one can show that possibility of a Taylor’s expansion leads to the relaxation law $`\delta \sigma \mathrm{log}t`$.

## V Conclusions

A novel computational algorithm is presented which permits one to simulate the energy, density of states and conductivity of a system with localized interacting electrons during times of the order of 40$`\mu `$s. We argue that during this time the relaxation of the system to some pseudoground state is completed. An analytical theory of the energy relaxation which is in a good agreement with the computational data is presented. Our computational results for the conductivity exhibit two very distinct time scales: the first is a very short time scale corresponding to the average value of the conductivity, the second a very long time scale defined by the long time relaxation of the conductivity. We attribute these two scales to the following physical picture, in which the relatively short time scale is a consequence of the relaxation of the system within one PS, well described by our analytical theory, as well as by our simulations, whereas the long time relaxation of the conductivity, is related to transitions between different pseudoground states. The microscopic origin of this huge time scale separation can be attributed to the fact that whereas many-electron transitions are not important for the VRH conductivity within one PS they play an important role in slowing down transitions between PS.

Current computational resources are not able to confirm this theory by directly observing the long relaxation processes, rather they are limited to the range of the short relaxation time scales, which can be simulated. In order to check our theory we have studied the conductivities of the different pseudoground states to see whether or not they are different. We employed two different models: the lattice model and the random site model, with random sites and no external disorder. We have shown that these two different models which correspond to different realizations of disorder lead to different physical effects. For the lattice model no difference in the conductivities of the different pseudoground states has been found. We understand this result in terms of universality of the Coulomb gap and ES hopping conduction.

For the model with random distances between sites and no external disorder we have found the difference of the conductivities to be within 10-12% which is large enough to explain the experimental data. We have shown also that the density of states in this model is not universal and that hopping conductivity does not obey the ES law. We think that similar effects might be observed in the lattice model with $`A=1`$ as well as at lower temperatures than those which we are able to achieve. With increasing $`A`$ this temperature range should become lower. Thus we think that universality of the Coulomb gap, that manifests itself in the ES law for the VRH, suppresses the long time relaxation because in this case the conductivities of different pseudoground states are very close to each other. The work has been funded by the US-Israel Binational Science Foundation Grant 9800097. The computations have been made in CHPC of the University of Utah. A.E. is grateful to Aspen Center for Physics, where this paper has been presented for the first time, and to Boris Shklovskii and Zvi Ovadyahu for important comments.

## VI Appendix. Computational algorithm.

To perform the simulation of the transport and thermodynamic properties on the finite array $`L\times L`$ for both the lattice model and the model with random spatial site distribution we use the periodic boundary conditions on a torus. In fact this means, that for the pair of sites $`i`$ and $`j`$ the distance between them is given by $`r_{ij}=[(\mathrm{\Delta }x_{ij})^2+(\mathrm{\Delta }y_{ij})^2]^{1/2}`$, where

$$\mathrm{\Delta }x_{ij}=min(|x_ix_j|,L|x_ix_j|)\mathrm{\Delta }y_{ij}=min(|y_iy_j|,L|y_iy_j|).$$
(19)
Here $`x_i`$ and $`y_i`$ are the sets of the site coordinates, which form a lattice in one model and are random in the other.

To simulate the conductivity one should add a term $`_iEx_i`$ to the Hamiltonian Eq.(2) where $`E`$ is a weak electric field. Due to the field the current flows around the torus in $`x`$ direction. It is convenient to calculate the total dipole moment due to electron transitions in the direction of electric field and obtain the conductivity from the equation

$$\sigma =\frac{1}{EL^2}\frac{dP}{dt},$$
(20)
On average the dipole moment $`P`$ increases linearly with time.

To find $`dP/dt`$ one needs a kinetic Monte Carlo(MC) algorithm, that connects the number of MC steps with a real time $`t`$. Note, that for any thermodynamic calculations the time is irrelevant. The starting equation is a transition rate for a single electron hop from site $`i`$ to site $`j`$, that has a form

$$\mathrm{\Gamma }_{ij}=\gamma \theta _{ij}\frac{\mathrm{exp}(2r_{ij}/a)}{1+\mathrm{exp}(ϵ_{ij}/T)},$$
(21)
where $`r_{ij}`$ is the distance between the sites $`i`$ and $`j`$, $`ϵ_{ij}`$ is the energy difference between the two configurations $`ϵ_{ij}=ϵ_iϵ_j1/r_{ij}E\mathrm{\Delta }x_{ij}`$, $`ϵ_i=\varphi _i+_j\frac{1}{r_{ij}}\left(n_j\nu \right)`$, $`\mathrm{\Delta }x_{ij}`$ is given by Eq.(19) and $`\theta _{ij}`$ is equal to $`1`$ if the site $`i`$ is occupied and site $`j`$ is empty or $`0`$ otherwise. The transition rate should have the dimensionality of frequency, it is written in a dimensionless form, assuming our time unit is $`\gamma ^1`$. The MC process can be started with any initial set distribution of the occupation numbers $`n_i`$ which evolves during the simulation.

There are two different algorithms developed for this type of computer simulations. The first one implies the calculation of all the transition rates $`\mathrm{\Gamma }_{ij}`$ in the system at each Monte Carlo (MC) step. Then the probability that the next transition to occur is $`ij`$ is given by

$$\frac{\mathrm{\Gamma }_{ij}}{_i_{ji}\mathrm{\Gamma }_{ij}}.$$
(22)
Then at each MC step the code chooses the transition with the above probability and performs it. This means changing the occupation numbers $`n_i`$ and $`n_j`$, calculating the contribution of this transition into the total dipole moment $`P`$ and recalculating all site energies and transition rates. After that the code comes to the next MC step.

For the above algorithm the physical time $`\mathrm{\Delta }t`$ for each MC step is

$$\mathrm{\Delta }t=(\underset{i}{}\underset{jj}{}\mathrm{\Gamma }_{ij})^1,$$
(23)
because in the real system all processes run simultaneously. Note, that $`\mathrm{\Delta }t`$ depends on the configuration of the system and varies during the simulation.

In the second algorithm at each MC step a pair of sites $`(i,j)`$ is chosen with equal probability from all possible sets. Then the transition $`(ij)`$ is accepted with the probability $`\mathrm{\Gamma }_{ij}`$. If the transition is rejected then the MC step is over and code proceeds to the next MC step. If it is accepted the transition is performed. This means that code changes occupation numbers $`n_i`$ and $`n_j`$, calculates the contribution of this transition into the total dipole moment $`P`$, and recalculate all the site energies $`ϵ_{ij}`$. This is the end of the MC step. In this case, the physical time per one MC step is constant and equal to

$$\mathrm{\Delta }t=\frac{1}{N_{tr}}=\frac{1}{L^2(L^21)},$$
(24)
$`\mathrm{\Delta }t=1/N_{tr}`$, where $`N_{tr}`$ is the total number of different transitions in the system.

The advantage of the first algorithm is that any MC step is successful: as a result of each step one electron moves from one site to another. The disadvantage is that for the interacting systems at each MC step the computer should recalculate $`N_{tr}`$ values of $`\mathrm{\Gamma }_{ij}`$. Therefore each MC step is very time consuming as compared to the MC step of the second algorithm. However, the disadvantage of the second algorithm is that at strong dispersion of the transition rates (small $`a`$ or $`T`$) the probability of the rejection is very high. In other words, the physical time $`\mathrm{\Delta }t`$ per one MC step is much smaller than in the first algorithm. This can be seen from Eqs. (23,24). Indeed, the double sum in Eq. (23) contains $`N_{tr}`$ terms. However, the majority of these terms are very small.

In this paper, we used a mixed scheme which combines the best features of both algorithms discussed above. We show that it is very efficient in the VRH regime for the interacting electrons. The original idea for this algorithm belongs to O.Biham ofer .

The transition rate Eq. (21) for the VRH can be written as a product $`\mathrm{\Gamma }_{ij}^T\mathrm{\Gamma }_{ij}^A`$, where the ’tunneling’ part of the transition rate $`\mathrm{\Gamma }_{ij}^T=\mathrm{exp}(2r_{ij}/a)`$, while $`\mathrm{\Gamma }_{ij}^A=\theta _{ij}/(1+\mathrm{exp}(ϵ_{ij}/T))`$ reflects activation. It is important now that $`\mathrm{\Gamma }^T`$ is independent of the configurations of electrons and should be calculated only once. Since the probabilities of tunneling and activation are independent we may apply the first algorithm with $`\mathrm{\Gamma }^T`$ and the second one with $`\mathrm{\Gamma }^A`$. Practically it means that we choose a pairs $`(i,j)`$ with the probabilities

$$\frac{\mathrm{\Gamma }_{ij}^T}{_i_{ji}\mathrm{\Gamma }_{ij}^T}$$
(25)
and accept the transition $`ij`$ with the probability $`\theta _{ij}/(1+\mathrm{exp}(ϵ_{ij}/T))`$. If the transition is rejected then the MC step is over and the code proceeds to the next MC step. If it is accepted, the transition is performed. To finish this step the code changes the occupation numbers $`n_i`$ and $`n_j`$, calculates the contribution of this transition into the total dipole moment $`P`$, and recalculate all the site energies $`ϵ_i`$. In this case, the physical time per one MC step is constant and equal to

$$\mathrm{\Delta }t=(\underset{i}{}\underset{ji}{}\mathrm{\Gamma }_{ij}^T)^1.$$
(26)
For the lattice model the sum $`_{ji}\mathrm{\Gamma }_{ij}^T`$ is independent on $`i`$ thus $`\mathrm{\Delta }t=(L^2_{ji}\mathrm{\Gamma }_{ij}^T)^1`$. Using this result the conductivity in units of $`\gamma a_0`$ can be written in a form

$$\sigma =\frac{1}{EL^2}\frac{dP}{dt}=\frac{P_{ji}\mathrm{exp}(2r_{ij}/a)}{EN_{MC}},$$
(27)
where $`P`$ is the total dipole moment due to the electron transitions after $`N_{MC}`$ steps.

For the model with random spatial site distribution the $`_{ji}\mathrm{\Gamma }_{ij}^T`$ is different for each site $`i`$ and the conductivity is given by

$$\sigma =\frac{P_i_{ji}\mathrm{exp}(2r_{ij}/a)}{EN_{MC}L^2}.$$
(28)
Now we compare the efficiency of all three algorithms. The efficiency is the number of operations which are necessary to simulate a physical process during a time $`t`$. The most important is how the efficiency depends on the number of sites in the system $`N=L^2`$. The results are shown in the Table 1. In the table $`\mathrm{\Delta }t`$ is the physical time correspondent to one MC step, $`N_{OP}/N_{MC}`$ is the number of operations per MC step, and the efficiency is given by $`Eff=N_{OP}/N_{MC}\mathrm{\Delta }t`$.

The time efficiency of the mixed algorithm is the same for both lattice and random site models. However, the memory requirements are much harder for the random site model. In this model one need to calculate all $`N^2=L^4`$ tunneling terms and have access to all of them at each MC step, because the transition at each step is chosen with the above weights. While in the lattice model there is only $`N=L^2`$ different tunneling terms $`\mathrm{exp}(2r_{ij}/a)`$ that have to be stored. In fact, this constraint does not allow us to simulate a system in which the number of sites exceeds 5000 employing the random site model.

One can see from the Table 1 that the efficiency of the mixed algorithm is the best. The use of this algorithm allowed us to simulate the macroscopic conductivity. In fact this algorithm has been used in Refourprl .


format: text

text:
# Semiclassical Distorted Wave Model Analysis of Backward Proton Emission from (𝑝,𝑝'𝑥) Reactions at Intermediate Energies

## 1 Introduction

Preequilibrium processes in nuclear reactions at intermediate energies are known to be dominated by multistep direct (MSD) processes. Among the models proposed in the years, the semiclassical distorted wave model  has proved its efficiency to describe MSD, especially for reactions at intermediate energies. It is based on the DWBA expansion of $`T`$-matrix elements and the cross section formula has no free adjustable parameter allowing a simple intuitive interpretation. The SCDW model has been applied to analyses of inclusive cross sections and spin observables in $`(p,p^{}x)`$ and $`(p,nx)`$ reactions over incident energy region of 60 to 400 MeV and wide range of mass number . Recently, several extensions and improvements of the SCDW model have been made. They include calculations up to three-step processes , incorporation of a single-particle (s.p.) wave function of target nuclei by using the Wigner transform of a one-body density matrix (OBDM) , introduction of phenomenological effective mass $`m^{}`$ of a nucleon in the target nucleus  and analyses of spin observables  which contain important information on an effective interaction in nuclear medium. As a result, SCDW model calculations can well reproduce the experimental data of double differential inclusive cross sections (DDX’s) for $`(p,p^{}x)`$ and $`(p,nx)`$ reactions. The model, however, underestimated the cross sections at very backward angles. This suggests that improvement of the SCDW model might be achieved by the inclusion of the high momentum components of the target nucleons.

The problem of backward proton production in proton-nucleus collisions has been a subject of much interest for several years. The main difficulty arises from the fact that the reaction involves a very large momentum transfer with only a small energy transfer. One of the models which have been applied to solve the problem is based on the correlated two-nucleon-cluster picture . A good description of the energy spectra even at the very high energy backward protons is achieved with this model. The main reason is that a phenomenological momentum distribution with high momentum components was employed . These components were found to be quite sensitive to the problem of backward proton production and were considered to be an effect due to many-nucleon correlations. Within the model  single scattering and two-nucleon cluster mechanisms at backward angles have been considered. Although there are some similarities between these contributions and corresponding one- and two-step processes, the advantage of the SCDW model is that the contribution from three-step process can be also obtained. Moreover, in SCDW model the distortion effects for the projectile and observed protons are carefully considered within DWBA while the calculated single scattering mechanism in  has evaluated the distortion effect approximately by the Glauber theory.

In Ref.  the SCDW model was modified so that realistic s.p. wave functions in a finite range potential can be used in terms of the Wigner transform of OBDM. Single-particle models with harmonic oscillator and Woods-Saxon potentials were used to describe the nuclear states for 90Zr$`(p,p^{}x)`$ reactions at 80 and 160 MeV instead of the local density Fermi gas (LFG) model used in previous analyses . It was concluded in  that the SCDW model with realistic s.p. wave functions gives larger one-step cross sections at backward angles, which result in better agreement with the experimental data over a wider angular range than the model with LFG. The calculations still somewhat underestimated the cross section at very large angles. Therefore, apart from the necessity of involving contributions of higher MSD processes the inclusion of higher momentum components of target nucleons becomes apparent. In the present paper we extend the SCDW model in terms of the Wigner transform of OBDM by considering the high-momentum components of nucleon momentum distributions deduced from several theoretical methods. This extended SCDW model is applied to 12C$`(p,p^{}x)`$ reaction at 150 MeV , 40Ca$`(p,p^{}x)`$ reaction at 186  and 392 MeV  and 90Zr$`(p,p^{}x)`$ at 160 MeV  and comparison with the experimental data is made.

The aim of our work is to analyze the angular distributions of the backward scattered protons through testing of different models for the nuclear states which account for nucleon-nucleon (NN) correlations. Such a systematic study could lead to a further refinement of the SCDW model and to a better agreement of the DDX with the experimental data at backward angles due to a proper inclusion of high-momentum components of the nuclear Fermi motion.

The paper is organized as follows. A brief formulation of the SCDW model and the proposed WT within different theoretical models are given in Sections II and III, respectively. The results of the calculations are presented and discussed in Section IV. The summary of the present work is given in Section V.

## 2 Outline of the SCDW model

The formulation of the SCDW model used in the present analysis has been described in details elsewhere , therefore only the outline is mentioned here.

Let us consider the one-step process in $`(N,N^{}x)`$ reactions in which a target nucleon is excited from an initial single particle state $`\alpha `$ with the energy $`\epsilon _\alpha `$ to a final one $`\beta `$ in the continuum with the energy $`\epsilon _\beta `$. A final expression of the inclusive $`(N,N^{}x)`$ double differential cross section for the one-step process is given by

$$\frac{^2\sigma ^{(1)}}{E_f\mathrm{\Omega }_f}=\left(\frac{A}{A+1}\right)^2𝑑𝐫\frac{k_f/k_f(𝐫)}{k_i/k_i(𝐫)}|\chi _i^{(+)}(𝐫)|^2|\chi _f^{()}(𝐫)|^2\underset{\tau }{}\left(\frac{^2\sigma }{E_f\mathrm{\Omega }_f}\right)_𝐫^\tau \rho _\tau (𝐫),$$
(1)
where $`A`$ is the target mass number, $`k_c`$ and $`k_c(𝐫)`$ ($`c`$=$`i`$ or $`f`$) the wave number at infinity and the local wave number in the initial ($`i`$) and final ($`f`$) channels, $`\chi _i^+`$ ($`\chi _f^{}`$) the distorted wave in the initial (final) channel, and $`\rho _\tau (𝐫)`$ ($`\tau `$=$`p`$ or $`n`$) the nucleon density for proton ($`p`$) or neutron ($`n`$) as the struck nucleon.

The local average nucleon-nucleon (NN) scattering cross section for collision between a leading particle, $`N`$, and a target nucleon, $`\tau `$, is given by

$`\left({\displaystyle \frac{^2\sigma }{E_f\mathrm{\Omega }_f}}\right)_𝐫^\tau `$ $`=`$ $`{\displaystyle \frac{2}{(2\pi )^3}}{\displaystyle \frac{1}{\rho _\tau (𝐫)}}{\displaystyle \frac{k_f(𝐫)}{k_i(𝐫)}}{\displaystyle 𝑑𝐤_\alpha 𝑑𝐤_\beta f_h^\tau (𝐤_\alpha ,𝐫)[2f_h^\tau (𝐤_\beta ,𝐫)]}`$ (2)
$`\times `$ $`\left({\displaystyle \frac{\sigma }{\mathrm{\Omega }}}\right)_{N\tau }\delta (𝐤_f(𝐫)+𝐤_\beta 𝐤_i(𝐫)𝐤_\alpha )\delta (\epsilon _\beta \epsilon _\alpha \omega ),`$
where $`𝐤_\alpha `$($`𝐤_\beta `$) is the nucleon momentum for a single-particle state with the energy $`\epsilon _\alpha `$($`\epsilon _\beta `$), $`\omega `$ the energy transfer, ($`\sigma /\mathrm{\Omega }`$) the two-nucleon scattering cross section and $`f_h^\tau (𝐤,𝐫)`$ is the Wigner transform of the OBDM for the hole states of proton and neutron. In Eq. (2) the sum rule of the WT for the hole and particle states is used:

$$f_h(𝐤,𝐫)+f_p(𝐤,𝐫)=2,$$
(3)
where $`f_p(𝐤,𝐫)`$ stands for the WT for the particle states. The factor $`f_h^\tau (𝐤_\alpha ,𝐫)[2f_h^\tau (𝐤_\beta ,𝐫)]`$ in Eq. (2) represents the ”probability” of the momentum state $`\mathrm{}𝐤_\alpha `$ being occupied and $`\mathrm{}𝐤_\beta `$ being unoccupied. The Pauli blocking effect is reflected in this factor. The difference between the single-particle energies in the second $`\delta `$-function in the right-hand side of Eq. (2) is given by

$$\epsilon _\alpha \epsilon _\beta =\frac{\mathrm{}^2}{2m^{}(𝐫)}[k_\beta ^2(𝐫)k_\alpha ^2(𝐫)],$$
(4)
using the effective mass of a nucleon, $`m^{}(𝐫)`$ .

The extension to higher-step processes is straightforward and the same final expressions as in Ref.  with the local average NN scattering cross sections given by Eq. (2) for successive collision points are also deduced for two- and three-step processes.

## 3 Wigner transform of one-body density matrix using different theoretical models

We introduce formulae for the Wigner transform (WT) derived in some theoretical methods, namely calculated with WS s.p. wave function (denoted as WS), with natural orbitals (NO) obtained from the Jastrow correlation method (JCM) which accounts for short-range correlations (SRC) in the case of 40Ca (denoted as JCM), from the approach based on the local density approximation (LDA) including the effects of SRC (denoted as LDA+SRC) and from the coherent density fluctuation model (denoted as CDFM).

### 3.1 WT with WS

For a single-particle potential, for instance Woods-Saxon potential, the WT for the hole states of proton or neutron is given by

$$f_h^\tau (𝐤,𝐫)=\underset{nlj}{}\frac{2j+1}{2l+1}_0^{\mathrm{}}𝑑𝐬e^{i𝐤𝐬}\underset{m}{}\varphi _{nlmj}^\tau (𝐫+𝐬/2)\varphi _{nlmj}^\tau (𝐫𝐬/2),$$
(5)
where $`\varphi _{nlmj}^\tau (𝐫)`$ is the s.p. wave function for the target nucleon and the sum runs over all the occupied orbits $`nlj`$ of protons or neutrons.

### 3.2 WT with JCM

A model independent way to define a set of single-particle wave functions and occupation probabilities uniquely from the correlated OBDM $`\rho (𝐫,𝐫^{})`$ is to use its natural orbital representation (NOR)

$$\rho _{NOR}(𝐫,𝐫^{})=\underset{\alpha }{}n_\alpha \psi _\alpha ^{}(𝐫)\psi _\alpha (𝐫^{}).$$
(6)
The normalized eigenfunctions $`\psi _\alpha (𝐫)`$ of $`\rho (𝐫,𝐫^{})`$, the so-called natural orbitals, form a complete orthonormal set. The associated eigenvalues $`n_\alpha `$ called natural occupation numbers define the probability ($`0n_\alpha 1`$) of the natural orbital $`\psi _\alpha (𝐫)`$ occupation in the ground state $`\mathrm{\Psi }`$. Within the NOR the Wigner transfrm of $`\rho (𝐫,𝐫^{})`$ can be derived as follows:

$$f(𝐤_c,𝐫)=\underset{c}{}n_ch_c^{NO}(𝐤_c,𝐫)=\underset{c}{}n_c_0^{\mathrm{}}𝑑𝐬e^{i𝐤_c𝐬}\psi _c(𝐫+𝐬/2)\psi _c^{}(𝐫𝐬/2)$$
(7)
We would like to note that the summation in Eq. (7) is taken over all hole- as well as particle-state NO included in the calculation of the Wigner transform. In the present work we use the results for the natural orbitals and for the occupation numbers in 40Ca obtained in Ref.  within the JCM in its low-order approximation. As it has been already done in  for the s.p. wave functions, in the present work each NO $`\psi _c(𝐫)`$ is expanded in terms of the Gaussian-type basis function

$$\psi _c(𝐫)=\underset{\upsilon =1}{\overset{N}{}}a_\upsilon ^{(c)}\mathrm{exp}(\kappa _\upsilon ^2r^2),$$
(8)
where $`N`$ is the number of the basis functions, $`\kappa _\upsilon `$ is given by a geometrical progression , $`\kappa _\upsilon =\kappa _1(\kappa _N/\kappa _1)^{(\upsilon 1)/(N1)}`$. The $`\{N,\kappa _1,\kappa _N\}`$ are the input parameters. The expansion coefficients $`a_\upsilon ^{(c)}`$ are determined by linear fitting of the basis functions (8) to the numerical data for the natural orbitals. Then each partial Wigner transform is calculated in accordance with Eq. (2.39) in Ref.  and after applying Eq. (7) the total Wigner transform for the NOR of $`\rho (𝐫,𝐫^{})`$ is obtained.

### 3.3 WT with LDA+SRC

It is well known that the inclusion of correlations in nuclear matter modifies the occupation probability predicted by the local density Fermi gas model. According to Ref.  one can introduce Wigner transform in a general way:

$$f(k_F(𝐫),k)=\mathrm{\Theta }(k_F(𝐫)k)+\delta f(k_F(𝐫),k),$$
(9)
where $`\mathrm{\Theta }(k_Fk)`$ corresponds to the WT of the LFG model, while $`\delta f(k_F,k)`$ is entirely due to the effects of dynamical correlations induced by the NN interaction. The local Fermi momentum $`k_F(r)`$ is related to the mass density through the relation

$$k_F(r)=\left[\frac{3}{2}\pi ^2\rho (r)\right]^{1/3}.$$
(10)
By definition of $`k_F(r)`$ one has $`\delta f(k_F(𝐫),k)𝑑𝐤=0`$. It was shown in  that for a finite nucleus the separation between mean-field contribution and correlation effects can be performed in an analogous way. For convenience more phenomenological procedure based on the results of the lowest order cluster (LOC) approximation developed in  has been followed to evaluate explicitly the correlated term. Choosing a correlation function of the form

$$f(r)=1e^{\beta ^2r^2},$$
(11)
the LOC gives for $`\delta f(k_F,k)`$

$$\delta f(k_F(𝐫),k)=\left[Y(k,8)k_{dir}\right]\mathrm{\Theta }(k_F(𝐫)k)+8\left\{k_{dir}Y(k,2)[Y(k,4)]^2\right\},$$
(12)
where

$$c_\mu ^1Y(k,\mu )=\frac{e^{\stackrel{~}{k}_+^2}e^{\stackrel{~}{k}_{}^2}}{2\stackrel{~}{k}}+_0^{\stackrel{~}{k}_+}e^{y^2}𝑑y+sgn(\stackrel{~}{k}_{})_0^{|\stackrel{~}{k}_{}|}e^{y^2}𝑑y$$
(13)
with

$$c_\mu =\frac{1}{8\sqrt{\pi }}\left(\frac{\mu }{2}\right)^{3/2},\stackrel{~}{k}=\frac{k}{\beta \sqrt{\mu }},\stackrel{~}{k}_\pm =\frac{k_F\pm k}{\beta \sqrt{\mu }},sgn(x)=\frac{x}{|x|}.$$
(14)
The quantity

$$k_{dir}=\frac{2k_F^3}{3\pi ^2}(f(r)1)^2𝑑𝐫=\frac{1}{3\sqrt{2\pi }}\left(\frac{k_F}{\beta }\right)^3$$
(15)
is the direct part of the Jastrow wound parameter.

### 3.4 WT with CDFM

The CDFM has been suggested in  as a model for studying characteristics of nuclear structure and nuclear reactions based on the local density distribution as a variable of the theory and using the essential results of the infinite nuclear matter theory. The model is introduced using the main ansatz of the generator coordinate method for the many-body function and the delta-function approximation for the overlap and energy kernels of the corresponding integral equation for the weight function. In the CDFM the Wigner distribution function can be written in the form:

$$f(𝐫,𝐤)=_0^{\mathrm{}}|f(x)|^2\mathrm{\Theta }(x|𝐫|)\mathrm{\Theta }(k_F(x)|𝐤|)𝑑x,$$
(16)
where the weight function $`f(x)`$ in the generator coordinate method is determined under the condition

$$_0^{\mathrm{}}|f(x)|^2𝑑x=1.$$
(17)
In the case of monotonically-decreasing density distributions $`(d\rho /dr0)`$ one can obtain a relation of the weight function $`f(x)`$ with the density distribution:

$$|f(x)|^2=\frac{1}{\rho _0(x)}\frac{d\rho (r)}{dr}|_{r=x},$$
(18)
where

$$\rho _0(x)=3A/4\pi x^3$$
(19)
and the generator coordinate $`x`$ is the radius of a sphere containing all $`A`$ nucleons uniformly distributed in it (the so-called ”flucton”).

## 4 Results of calculations and discussion

The input data used in the calculations of Wigner transforms can be summarized as follows. The single-particle wave function in Eq. (5) were calculated using a Woods-Saxon potential with the radius parameter $`r_0`$=1.27 fm and the surface diffuseness parameter $`a`$=0.67 fm, including the isovector term and the Coulomb potential. The set $`\{N,\kappa _1,\kappa _N\}`$=$`\{20,1,0.1\}`$ is taken for the WS potential and for the WT with NOR. For the correlation factor in Eq. (11) we adopt the same value $`\beta `$=1.1 fm -1 as in Ref.  which is taken from the microscopic nuclear matter calculations  but reproduces also well the results for the nucleon momentum distribution $`n(k)`$ of finite nuclei. This value leads to $`k_{dir}`$=0.23 \[Eq. (15)\]. In the CDFM method as well as in the LDA approach with SRC we use the nucleon density distribution $`\rho (r)`$ of the Woods-Saxon shape with Negele’s geometrical parameters . The proton and the neutron densities are assumed to be given by $`(Z/A)\rho (r)`$ and $`(N/A)\rho (r)`$, respectively. The local Fermi momentum $`k_F(r)`$ was obtained from Eq. (10) using the same density distribution.

As for the other input data for the SCDW calculations, we use basically the same ones as in Ref. . The global optical potentials based on Dirac phenomenology of Hama et al.  for protons and Ishibashi et al.  for neutrons were used. For a simplicity of the numerical calculations we neglect the spin-orbit coupling in the distorting potentials. The nonlocality correction for distorted waves and WS s.p. wave functions was taken into account using the Perey factor with range 0.85 . An effective NN interaction in terms of $`G`$ matrix parameterized by the Melbourne group  was employed. The two-nucleon scattering cross section in Eq. (2) is calculated using this $`G`$ matrix. In addition to , an effective mass of a target nucleon $`m^{}`$ is used. We assume a simple WS form for the $`r`$-dependence of $`m^{}`$. The latter at $`r`$=0 was chosen to be $`m^{}=0.8m`$ with a bare mass $`m`$.

Figure 1 shows the SCDW double differential cross sections with CDFM, LDA+SRC and WS for the 12C$`(p,p^{}x)`$ reaction at 150 MeV incident and 80 MeV outgoing energies compared with two sets of experimental data . A good overall agreement with the data including backward angle region is obtained in the case of CDFM. The SCDW cross section calculated with WS drops down at backward angles, while those with LDA+SRC overestimates the experimental data. As will be discussed later, the main reason for this behavior is the different account for the high-momentum components of the nuclear Fermi motion. In Fig. 1 the contributions of individual multistep processes when CDFM is used for the same reaction are also plotted. One can see that proton emission via one-step process is dominant in the angular region from 20 to 50. Contributions of two- and three-step processes become appreciable with increasing angles. The role of higher-step processes has been already seen in previous calculations .

Comparisons between the cross sections calculated by using the Wigner transforms obtained with all theoretical methods considered and the measured cross sections for $`(p,p^{}x)`$ on 40Ca at 186 MeV and 392 MeV incident energies are given in Fig. 2 and Fig. 3, respectively. Both figures confirm the fact that the cross sections at large angles depend strongly on the theoretical model used to describe the nuclear states. Fig. 2 shows a fair agreement of the CDFM and JCM results with the data for emission energy of 110 MeV and to less extent for energy of 130 MeV. At the same time, for emission energy of 150 MeV the WS result is able to reproduce the angular distribution at backward angles. Similar behavior can be seen from Fig. 3, where the double differential cross sections calculated with CDFM and that with JCM are in a good agreement with the experimental data. However, measurements over an extended range of backward angles can test the various methods and, hence, the role of NN correlations.

The improved description of the backward proton emission with the SCDW model by the use of realistic models for nuclear states can be explained by the difference in the momentum distribution of the target nucleons. Figure 4 shows the momentum distributions $`n(k)`$ of 12C and 40Ca with CDFM, LDA+SRC, WS and JCM. One sees that $`n(k)`$ calculated with all methods but not with the WS contain high-momentum components at momenta larger than 1.5 fm-1. The LDA+SRC  and Jastrow method  account for short-range correlations, while the CDFM takes into account mostly long-range correlations (LRC). The larger high-momentum component due to the stronger SRC when using LDA+SRC causes the larger cross sections at backward angles. Although $`n(k)`$ calculated in this method fits well the experimental $`y`$-scaling data for 12C  at large momenta, the corresponding cross sections overestimate the data. Here we would like to note that the $`n(k)`$ of 40Ca obtained after integration of the total WT $`f(𝐤,𝐫)`$ with JCM over $`𝐫`$ is not given in Fig. 4 and it underestimates the one obtained in  at $`k2`$ fm-1. This is particularly due to the truncation of the particle-state NO’s included in the calculation of the WT of 40Ca. On the other side, due to the lack of enough high-momentum components the angular distributions calculated with WS underestimate the data at backward angles. The nucleon momentum distributions obtained with CDFM for both nuclei presented in Fig. 4 also exhibit high-momentum components, though smaller than the ones corresponding to methods accounting for SRC. Fig. 4 clearly illustrates the strong dependence of backward proton production from $`(p,p^{}x)`$ reactions on the momentum distribution of target nucleons.

In Fig. 5 the SCDW cross sections for 90Zr$`(p,p^{}x)`$ at 160 MeV and emission energy of 120 MeV calculated with the CDFM, LDA+SRC and WS models are compared. Again a very good agreement with the experimental data over the forward, intermediate and backward angles is achieved when using the Wigner transform obtained from CDFM. In Fig. 6 the contributions of individual multistep processes for the same reaction and corresponding to the same models are plotted to show their variation with the scattering angle. It is seen that for this relatively large emission energy the one-step cross section calculated with LDA+SRC exceeds the experimental cross section for large angles. On the other hand, the one-step cross section obtained with WS drops at angles less than 140 underestimating the experimental cross section at backward angles.

From the analyses of high energy heavy ion collision and deep inelastic lepton-nucleus scattering, the following form of the nucleon momentum distribution has been proposed :

$$n(k)=N[\mathrm{exp}(k^2/p_0^2)+\epsilon _0\mathrm{exp}(k^2/q_0^2)+\epsilon _1\mathrm{exp}(k^2/q_1^2)],$$
(20)
where $`N`$ is a normalization factor and the values of the parameters are determined to be $`p_0=\sqrt{\frac{2}{5}}k_F`$, $`\epsilon _00.03`$, $`q_0\sqrt{3}p_0`$, $`\epsilon _10.003`$ and $`q_10.5`$ GeV/c. The first term of Eq. (20) corresponds to single-particle contribution to $`n(k)`$ in the mean-field approximation (MFA). It was noticed in  that the second term in Eq. (20) is due to the long-range correlations (L-type high momentum component) while the third term was supposed to be due to short-range correlations (S-type high momentum component). Moreover, as it was also noticed in Ref.  only the high-momentum component generated by the Hartree-Fock correlation should be included in the single scattering calculation. In order to reveal better the role of different NN correlations for the correct description of the backscattered proton spectra, we make in Fig. 7 a comparison between the nucleon momentum distributions calculated with CDFM and JCM and those calculated by using Eq. (20) in the case of 12C and 40Ca nuclei. Particularly, our interest is in the comparison of the high-momentum components which the different methods contain. It turns out from Fig. 7 that for $`k>1.5`$ fm-1 the curves corresponding to the inclusion only of the L-type term and CDFM have similar behavior (especially for 12C), while the results when including only the S-type term and Jastrow method are close to each other. It is seen also from Fig. 7 that the contribution of the S-type term dominates in the total $`n(k)`$ calculated by Eq. (20). Therefore, we suggest that the long-range correlations play an important role to the problem of backward proton production. This is confirmed by the good agreement obtained in  and by our results obtained with the use of Wigner transform with CDFM which lead to best description of the experimental cross sections over wide target muss number and emission energy regions.

## 5 Summary

The results of the present work can be summarized as follows:
i) Double-differential cross sections of 12C$`(p,p^{}x)`$, 40Ca$`(p,p^{}x)`$ and 90Zr$`(p,p^{}x)`$ reactions are calculated, using SCDW model and Wigner transform of the OBDM obtained with different approaches, in the incident energy range of 150-392 MeV. A good overall agreement with the experimental data is achieved including the backward angles region. These analyses show that the present SCDW model can describe successfully the MSD processes of $`(p,p^{}x)`$ reactions over the wide range of target mass number and incident energies.
ii) It is found that the high-momentum components of the nucleon momentum distribution are most responsible for describing the angular distributions of the backscattered protons. Since CDFM leads to the best agreement with the experimental cross sections at wide emission energy region, it is demonstrated that the long-range correlations which are related to the collective nucleon motion are rather important than the short-range ones generated in the other correlation methods.

A more detailed study on the backward proton emission in the framework of the CDFM is in progress.

## 6 Acknowledgments

The authors wish to thank Professor A.A. Cowley for providing us the experimental data from Ref. . One of the authors (MKG) is grateful to the JSPS’s Postdoctoral Fellowships for foreign researchers. This work was partly supported by a Grant-in-Aid for Scientific Research of the Ministry of Education, Science and Culture (No. 14002042).


text:
# Testing boundary conditions efficiency in simulations of long-range interacting magnetic models

## I Introduction

Boundary conditions are a central issue in almost every statistical mechanics calculation. In particular, periodic boundary conditions are the usual way of diminishing finite size border effects in translationally invariant systems. This point is of special importance in numerical simulations, where an appropriated handling of finite size effects may determine completely the quality of the results. While both the interpretation and implementation of periodic boundary conditions are straigthforward in systems with short range microscopic interactions, their usage in systems with long-range interactions is not so easy, with several subtleties to take into account. In this work we will consider magnetic systems where the iteractions between pairs of spins decay as $`1/r^\alpha `$, $`r`$ being the distance between spins and $`\alpha >0`$ (this includes, among the most conspicuous examples, Coulomb and dipolar interactions). Then, every spin interacts with any other spin in the system, with a slow decaying distance-dependent intensity. Finite size effects are much more stronger in these kind of systems than in systems with short range interactions, because every spin feels directly the influence of the border. Hence, boundary conditions are a topic of central importance in numerical simulations of these kind of systems.

Spin systems with the above type of interactions may present a very rich phenomenology, even in one dimension, and have attracted a lot of attention in the last years, both because its applications to ultrathin magnetic films systemsDeBell -Gleiser , glassy systems Grousson  and to the study of statistical mechanics fundamentalsCannas1 -Tamarit , among many others.

There are two alternative representations of a finite system with periodic boundary conditions. Let us think on a spin system defined on a hypercubic $`d`$-dimensional lattice with $`N=L^d`$ sites. On one hand, we can visualize the system as a d-torus (closed ring in $`d=1`$; torus in $`d=2`$, etc.), where a given spin interacts with its closer neighbors defined by the topology of the d-torus, up to a certain range of interaction. On the other hand, we can think that we have an infinite system, where the original finite system has been replicated infinite times in all the coordinate directions. In other words, we can think that the infinite lattice has been partitioned into cells of size $`N=L^d`$ and that we have chosen only the periodic solutions of the problem, with periodicity $`L`$ in all the coordinate directions.

While the difference between both views is just a matter of interpretation for systems with short-range interactions, the situation changes for long-range interacting systems. According to the second scheme, a given spin will interact with infinite replicas of everyone of the rest of the spins and we can express the effective interaction between two spins inside the system as an infinite sum over replicas (this includes a “self-interaction” term, that accounts for the interaction of every spin with its own replicas). We will call these infinite periodic boundary conditions (IPBC). On the other hand, in the closed topology of the first scheme it is consistent to consider that every pair of spins inside the system interacts only through its minimal distance over the hypertorus. This corresponds to the minimal truncation of the infinite series of the previous scheme and it is sometimes called the first image convention; we will refer to it as first image periodic boundary conditions (FIPBC). While both type of boundary conditions give the same result in the thermodynamic limit $`L\mathrm{}`$, they are clearly different for finite systems.

Now, the question is how good are both type of boundary conditions to wipe out finite size effects in systems of moderate sizes, such as those available in practical numerical calculations? This is of particular interest, considering that the numerical simulations are much more costly than in the case of short range interactions (the computational complexity of the algorithms is typically $`𝒪(N^2)`$). Although both IPBC and FIPBC have been used in several works, including both numericalDeBell -Cannas1 Luijten  and analyticalLee ; Campa  calculations, to the best of our knowledge the only comparison of its relative efficiency in magnetic systems was reported in an old work by Kretschmer and Binder Kretschmer , in the case of dipolar interactions in a three-dimensional system. In that work the authors observed a very poor performance of the FIPBC compared to the IPBC when the system exhibits a ferromagnetic phase. Since then, several authors assumed this result as a general rule, independent of the dimension $`d`$, the range of the interactions (i.e., of the exponent $`\alpha `$) and the symmetry of the low temperature phase. In this work we compare the relative efficiency of both type of boundary conditions by considering the equilibrium critical behaviour of ferromagnetic, antiferromagnetic and competing interactions in $`d=1`$ and $`d=2`$ and for different values of the exponent $`\alpha `$, using a single spin-flip Monte Carlo algorithm (Metropolis).

While the implementation of the FIPBC is straightforward, the implementation of the IPBC is much more cumbersome and needs care, because the infinite series for the effective interactions are usually slowly convergent. This is another reason for studying its relative efficiency. The usual way of handling numerically the IPBC in $`d=2`$ and $`d=3`$ is to adapt the Ewald sums technique, originally derived for systems of interacting charged particlesAllen ; Frenkel , to magnetic modelsDeBell ; Kretschmer . In section II we compare the critical temperature of one dimensional Ising models with ferromagnetic, antiferromagnetic and competing interactions, for different system sizes and $`1<\alpha <2`$ using FIPBC and IPBC. In this case we introduce a numerical technique for implementing the IPBC that is easier to implement in $`d=1`$ for arbitrary values of $`\alpha `$ than the Ewald sums technique. We also analized the critical behaviour of the antiferromagnetic model for $`0\alpha 1`$. In section III we compare the specific heat numerical results obtained around the critical region for ferromagnetic long range interactions with $`\alpha =3`$, and competing ferromagnetic short range interactions with antiferromagnetic (dipolar) long range interactions, using the Ewald sums for the IPBC. A general discussion of the results is presented in section IV.

## II One-dimensional Ising models

We considered first an Ising model of the type

$$H=J_0\underset{(i,j)}{}f(r_{ij})S_iS_jS_i=\pm 1$$
(1)
where the sum $`_{(i,j)}`$ runs over all distinct pairs of sites $`(i,j)`$ ($`i,j=1,2,\mathrm{}`$) in an infinite chain; $`r_{ij}|ij|`$ is the distance between sites and

$$f(r)=\{\begin{array}{ccc}\frac{1}{r^\alpha }\hfill & if\hfill & r1\hfill \\ 0\hfill & if\hfill & r=0\hfill \end{array}$$
(2)
This model presents a phase transition at a non-zero critical temperature $`T_c(\alpha )`$ for $`1<\alpha 2`$ Dyson ; Frohlich . In the ferromagnetic version, $`J_0>0`$, the critical temperature behaves asymptotically as $`k_bT_c/J_02\zeta (\alpha )2/(1\alpha )`$ when $`\alpha 1^+`$ Cannas4 ; Luijten , where $`\zeta (\alpha )=_{n=0}^{\mathrm{}}\frac{1}{n^\alpha }`$ is the Riemann zeta function. When $`1<\alpha 1.5`$ the critical exponents are classical. The thermodynamic limit in the ferromagnetic model is not defined when $`0\alpha <1`$ and the system presents non-extensive behaviourCannas1 ; however, with an appropriated regularization procedure it can be shown that the system still presents a phase transition where all the thermodynamic functions (not only the critical exponents) obey mean-field behaviourCannas2 ; Lee . These results generalize to ferromagnetic models in arbitrary dimension $`d`$ with $`\alpha =d+\sigma `$: the thermodynamic limit is well defined for $`\sigma >0`$. The critical exponents are classical for $`0<\sigma <d/2`$ and assume continuously varying values for $`d/2\sigma 2`$. When $`\sigma >min(d,2)`$ the system exhibits short range critical behaviourAinzenman .

In this section we considered a finite chain with $`1<\alpha 2`$; the finite size version of the Hamiltonian (1) is

$$H=J_0\underset{(i,j)}{}W(r_{ij})S_iS_jS_i=\pm 1$$
(3)
where now the sum $`_{(i,j)}`$ runs over all distinct pairs of sites $`(i,j)`$ ($`i,j=1,2,\mathrm{},L`$) in a chain of lenght $`L`$, $`r_{ij}|ij|=0,1\mathrm{}L1`$ and $`W(r_{ij})`$ is the effective interaction between spins $`(i,j)`$ resulting from the different types of periodic boundary conditions.

The FIPBC effective interaction $`W_F`$ is then simply defined as

$$W_F(r)=f\left(\mathrm{min}[r,|rL|]\right)$$
(4)
while the IPBC effective interaction $`W_I`$ is defined as

$$W_I(r)=\underset{n=0,\pm 1,\mathrm{}}{}f(r+nL)=\frac{1}{L^\alpha }\underset{n=0,\pm 1,\mathrm{}}{}\frac{1}{|n+r/L|^\alpha }forr=1,2,\mathrm{}L1$$
(5)
and

$$W_I(0)=\underset{n=\pm 1,\pm 2,\mathrm{}}{}f(nL)=\frac{1}{L^\alpha }\underset{n=\pm 1,\pm 2,\mathrm{}}{}\frac{1}{|n|^\alpha }=2\frac{\zeta (\alpha )}{L^\alpha }$$
(6)
where $`n=\pm 1,\pm 2,\mathrm{},\pm \mathrm{}`$ numerates the infinite replicas of the system $`n=0`$. Equation (5) is a slowly convergent series and cannot be truncated in a straitforward way. In order to obtain a controllable approximation we rewrite Eq.(5) as:

$`W_I(r)`$ $`=`$ $`{\displaystyle \frac{1}{r^\alpha }}+{\displaystyle \frac{1}{L^\alpha }}{\displaystyle \underset{n=\pm 1,\mathrm{},\pm \mathrm{}}{}}{\displaystyle \frac{1}{|n+r/L|^\alpha }}`$ (7)
$`=`$ $`{\displaystyle \frac{1}{r^\alpha }}+{\displaystyle \frac{1}{L^\alpha }}\left[{\displaystyle \underset{n=\pm 1,\mathrm{},\pm m}{}}{\displaystyle \frac{1}{|n+r/L|^\alpha }}+{\displaystyle \underset{n=m+1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{|n+r/L|^\alpha }}+{\displaystyle \frac{1}{|nr/L|^\alpha }}\right)\right]`$
Since $`r<L`$, for $`\alpha >1`$ and $`m1`$ we can approximate $`W_I(r)W_m(r)`$ for $`r=1,\mathrm{},L1`$, where

$`W_m(r)`$ $`=`$ $`{\displaystyle \frac{1}{r^\alpha }}+{\displaystyle \frac{1}{L^\alpha }}\left[{\displaystyle \underset{n=\pm 1,\mathrm{},\pm m}{}}{\displaystyle \frac{1}{|n+r/L|^\alpha }}+2{\displaystyle \underset{n=m+1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n^\alpha }}\right]`$ (8)
$`=`$ $`{\displaystyle \frac{1}{r^\alpha }}+{\displaystyle \frac{2\zeta (\alpha )}{L^\alpha }}+{\displaystyle \frac{1}{L^\alpha }}{\displaystyle \underset{n=\pm 1,\mathrm{},\pm m}{}}\left({\displaystyle \frac{1}{|n+r/L|^\alpha }}{\displaystyle \frac{1}{|n|^\alpha }}\right)`$
The finite approximation $`W_m(r)`$ can be calculated once for ever at the beginning of the simulation. In order to control the accuracy of the approximation, we can calculate the differences

$$\mathrm{\Delta }W_m(r)W_{m+1}(r)W_m(r)=\frac{1}{L^\alpha }\left[\frac{1}{|m+1+r/L|^\alpha }+\frac{1}{|m+1r/L|^\alpha }\frac{2}{(m+1)^\alpha }\right].$$
(9)
Noting that

$$\mathrm{\Delta }W_m(r)<\mathrm{\Delta }W_m(L)=\frac{1}{L^\alpha }\left(\frac{1}{(m+2)^\alpha }+\frac{2}{m^\alpha }\frac{2}{(m+1)^\alpha }\right)$$

we can estimate the error of the approximation by choosing $`m`$ such that $`\mathrm{\Delta }W_m(L)<ϵ`$; since for $`m1`$

$$\mathrm{\Delta }W_m(L)=\frac{1}{L^\alpha (m+1)^\alpha }\left[\frac{\alpha (\alpha +1)}{(m+1)^2}+𝒪\left(m^4\right)\right]$$

we have

$$m>\left(\frac{\alpha (\alpha +2)}{ϵL^\alpha }\right)^{\frac{1}{\alpha +2}}1$$

.

In all the simulations we choose $`m`$ such that $`ϵ<10^6`$ for every value of $`L`$ and $`\alpha `$.

We first considered the ferromagnetic case $`J_0>0`$ for $`1<\alpha 2`$. We calculated the critical temperature using both FIPBC and IPBC by means of the so called Binder cumulantBinder

$$U_L1\frac{M^4_L}{3M^2_L^2}$$

where $`M`$ is the order parameter (magnetization per spin). The critical temperatures were estimateded as the point of intersectionBinder  of the curves $`U_L^{}`$, $`U_{L\mathrm{"}}`$ vs. $`T`$ (in what follows we choose $`k_B=1`$), for pairs of values $`(L^{},L\mathrm{"})`$ chosen around some average $`L`$, with a dispersion around $`L`$ of less than 10%. The results for boh type of boundary conditions were calculated using the same single spin-flip Monte Carlo dynamics (Metropolis). The results are shown in Fig.1, where the exact asymptotic behaviour when $`\alpha 1`$ is also shown for comparison. For $`\alpha >1.5`$ the results are very similar, with an obtained difference between both types of boundary conditions of less than 5% and decreasing when $`\alpha 2`$. On the other hand, for $`\alpha <1.5`$ the difference in performance increases dramatically as $`\alpha `$ decreases. In table 1 the results for the best estimate using IPBC for $`\alpha 1.5`$ are compared with the best estimates in the literature obtained by Luijten and BlöteLuijten , using a specially designed block Monte Carlo algorithm for ferromagnetic systems with long range interactionsLuijten2  (IPBC are implicit in this algorithm) and chains of very large sizes ($`L=150000`$). While the IPBC gives a rather accurate result, even for a modest size of $`L100`$ (see table 1), the FIPBC with $`L1000`$ give an error that reaches around one order of magnitude for $`\alpha =1.02`$ (see Fig.1).

We next considered the antiferromagnetic case $`J_0<0`$ for $`1<\alpha 2`$ and system sizes $`L=300,500`$ and $`1000`$ using both type of boundary conditions, and for $`0\alpha 1`$ using FIPBC (IPBC cannot be implemented in this case). While the ferromagnetic version of this model has been extensively studied, to the best of our knowledge the only existing results for the antiferromagnetic case for $`1<\alpha 2`$ have been obtained by Romano (see Ref.Romano  and references therein). In particular, no previous results are known for $`0<\alpha 1`$, while the $`\alpha =0`$ case has been exactly solvedAl-Wahsh . We found that the system presents a low temperature ordered antiferromagnetic phase for every value of $`0<\alpha 1`$ below a finite critical temperature. At variance with the ferromagnetic case, the antiferromagnetic model shows a well defined thermodynamic limit for $`0<\alpha 1`$; this fact has also been observed in an antiferromagnetic model with long range interactions defined in the hypercubic cellFranco .

For this system the critical temperature was estimated as the maximum in the staggered susceptibility, defined as the fluctuations in the staggered magnetization. The results are shown in Fig.2. We see that both type of boundary conditions give in this case the same result (within the error bars) for every system size when $`1<\alpha 2`$. The results for $`L=1000`$ are in good agreement with previous estimationsRomano  using IPBC for $`\alpha =2`$: $`T_c/|J_0|0.219`$, and converge to $`T_c=0`$ for $`\alpha 0`$, reproducing the exact solutionAl-Wahsh .

Finally, we considered an Ising model with competing short range ferromagnetic interactions and long range antiferromagnetic interactions defined by the Hamiltonian

$$H=J_0\underset{<i,j>}{}S_iS_j+J_d\underset{(i,j)}{}f(r_{ij})S_iS_jS_i=\pm 1,$$
(10)
where $`J_0,J_d>0`$, the sum $`_{<i,j>}`$ runs over nearest neighbor paris of spins and $`f(r)`$ is given by Eq.(2). This is a one dimensional generalization of a model for an ultrathin magnetic film with competing ferromagnetic exchange interactions and dipolar interactions that will be described in the next section. The finite temperature phase diagram of this model when $`1<\alpha 2`$ is similar to the one corresponding to the two dimensional case. At low temperatures it presents different modulated phases (with zero global magnetization) of period $`2h`$, depending on the ratio $`\delta =J_0/J_d`$, with a finite critical temperature $`T_c=T_c(\alpha ,\delta )`$. The period of the modulation increases with $`\delta `$, ranging from $`h=1`$ (antiferromagnetic) for low values of $`\delta `$, up to the system size for very large values of $`\delta `$, passing by all the possible values commensurable with the system lenght $`L`$. The calculation of the global phase diagram of this model is in progress and the details will be published elsewhere.

For the comparison of the boundary conditions we choose some representative values of $`\alpha `$ and $`\delta `$; other values gave similar results. The results are shown in table 2. No significative differences were observed between both type of boundary conditions.

## III Two-dimensional Ising models

In this section we considered Ising models in a square lattice of $`L^2`$ sites with interactions that decay as $`1/r^3`$. IPBC were implemented in this case using the Ewald sums techniqueDeBell . We started with a Hamiltonian of the type (1) with $`\alpha =3`$ and $`J_0>0`$. The choice of the vale $`\alpha =3`$ obeys two reasons: first, it allows the usage of the Ewald sums for the implementation of the IPBC and second, because of its close relationship with dipolar interactions. To perform the comparison between boundary conditions we calculated in this case the specific heat

$$C(T)=\frac{1}{NT^2}\left(H^2_LH_L^2\right)$$

as a function of $`T`$ for different system sizes $`L`$. In Fig.3 we show the specific heat calculations for both types of boundary conditions for system sizes $`L=24,32,48,64`$ and $`96`$. Although FIPBC display larger finite size effects for small system sizes than the IPBC, both type of boundary conditions give similar results when $`L=96`$, as can be appreciated in Fig.4, where we compare the results of FIPBC and IPBC for the smallest and the largest sytem sizes. For $`L=96`$ the curves obtained wiht both kind of boundary conditions show clearly the emergency of a discontinuous shape, consistently with the expected classical behaviourLuijten .

Finally we considered the two-dimensional version of Hamiltonian (10). This model represents an ultrathin magnetic film where the axis of easy magnetization is oriented perpendicular to the plane of the filmDeBell . In this case the first term in Hamiltonian (10) represents ferromagnetic exchange interactions while the second represents dipolar interactions. The system presents a low temperature ordered phase below some critical temperature $`T_c=T_c(\delta )`$, composed by ferromagnetic stripes of width $`h=h(\delta )`$, so that spins belonging to adyacent stripes are antialigned. The overall known features of the finite temperature phase diagram of this model are described in Refs.MacIsaac -Gleiser . For the comparison of the boundary conditions we choose the value $`\delta =3`$, for which the low temperature phase corresponds to an $`h=4`$ striped phase, since it is representative of the global behaviour for intermediate values of $`\delta `$ and was previously analized in detail by Booth et alBooth  using Monte Carlo with IPBCnote . The results are shown in Fig.5 for $`L=24,32`$. We see that both types of boundary conditions give almost the same results even for rather small sizes, except very near the critical point (left peak; see Ref.Booth  for details) where the FIPBC give a little larger peak than the IPBC.

## IV Discussion

We have compared the finite size effects introduced by FIPBC and IPBC in one and two dimensional Ising models with different types of long range interactions. This allowed us to analize a large variety of situations of interest and the influence of different parameters. We have shown that, when the interactions are such that the system presents a low temperature ordered phase with zero global magnetization the finite size effects produced by both types of boundary conditions are almost the same. This can be understood if we consider that, in this case, the infinite sums in the local fields calculations converge fast and hence the finite size effects associated with the long range character of the interactions are wiped out quickly; hence, the truncation associated with the FIPBC is enough to take it into account. A similar effect was observed by Kretschmer and BinderKretschmer  in the antiferromagnetic phase of the three dimensional Ising model with nearest neighbor exchange interactions and long range dipolar interactions. On the other hand, when the system presents low temperature ferromagnetic ordering, large differences in the finite size effects produced by the different periodic boundary conditions can be expected, depending on the values of $`d`$ and $`\alpha `$. Moreover, our results suggest that the large differences appear when the ratio $`\alpha /d`$ is such that the system enters into the classical regime $`\alpha /d3/2`$ ($`\sigma <d/2`$), while for $`\alpha /d>3/2`$ both type of boundary conditions produce similar finite size effects, with a decreasing difference as $`\alpha /d`$ increases. This is also consistent with Kretschmer and Binder result Kretschmer  for $`d=3`$, $`\alpha =3`$.

We have also shown that in the ferromagnetic one dimensional model the finite size effects given by the FIPBC becomes very large when $`\alpha 1^+`$. This suggests that the same problem could be observed in the non-extensive regime $`0\alpha 1`$, where the FIPBC are the most direct form of implementing periodic boundary conditions. Moreover, in this regime when $`\alpha <d`$ the implementation of the IPBC is more involved because the infinite sums of the type (5) do not converge. Since in the $`\alpha =0`$ case we lost the boundary completely (this corresponds, with an appropriated re-scaling, to the Curie-Weiss model) we can expect that in the limit $`\alpha 0`$ the influence of periodic boundary conditions will be also negligible. Hence, we may expect the largest finite size effects when using FIPBC in the region $`\alpha /d1`$.

Finally, as a by-product of this work, we introduced a simple method for implementing the IPBC in one dimensional models for arbitrary values of $`\alpha >1`$.

This work was partially supported by grants from Consejo Nacional de Investigaciones Científicas y Técnicas CONICET (Argentina), Agencia Córdoba Ciencia (Córdoba, Argentina), Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba (Argentina), Fundação VITAE (Brazil) and CNPq (Brazil).

We wish to thank R. D az-M ndez and Roberto Mulet for their help in the implementation of the Ewald sums.


text:
# Comment on ”Estimate of the vibrational frequencies of spherical virus particles”

## Abstract

This comment corrects some errors which appeared in the calculation of an elastic sphere eigenenergies. As a result, the symmetry of the mode having the lowest frequency is changed. Also a direction for calculating the damping of these modes for embedded elastic spheres is given.

In a recent articleFordPRE03 , L. H. Ford discusses the normal modes of vibration of spherical virus particles using a liquid drop model and an elastic sphere model.

However the analysis errs in calculating the energy of the $`n=0`$ spheroidal mode with the elastic sphere model. This same error is recurrent in the literature eringen ; Fujii91 ; Tanaka93 ; Kuok03 . Some authors have since corrected themerratumFujii91 . The original Lamb paperlamb1882  was correct and an explanation is available elsewheresaviot96 . This error results in additional energies which do not correspond to any vibration eigenmode. It gives a wrong energy for the lowest $`n=0`$ spheroidal mode comparable to the one of the lowest $`n=2`$ spheroidal mode. For materials with positive Poisson ratio, it is impossible to have the energy for the fundamental $`n=0`$ mode smaller than the energy for the fundamental $`n=2`$ one.

Equation 3 in the paper is valid for $`n`$ different from zero. For $`n`$ equal to zero, $`S_{\mathrm{rp}}(0,x,y)=0`$ should be used instead (see equation 4). Moreover, expressions given in equations 4 and 7 are invalid: there is no square dependence and the sign of $`S_{\mathrm{ts}}`$ is wrong. With correct calculations, the lowest frequencies for the elastic sphere model are reached with $`n=2`$ modes and frequencies are a bit lower (9.2, 9.4 and 4.6 GHz for Nylon, Polystyrene and Polyethylene respectively).

We also would like to point out that these normal modes are damped when the virus particle is embedded inside a liquid. Estimations of damping can be made using the ”complex frequency” approach described elsewheredubrovskiy81 ; saviotPRB03 . When the virus is inside water for example, there is not much acoustic impedance mismatch at the surface. For this reason the normal modes will be broad and have short lifetimes. Therefore, the objective of killing viruses by sending out sound waves that resonate and destroy them is probably unworkable in such a configuration.