text:
# Calculation of the incremental stress-strain relation of a polygonal packing

## I Introduction

The non linear and irreversible behavior of soils has been described by different constitutive theories . Here the stress - strain relation is postulated using a certain number of material parameters that are measured in experimental tests. In practice, the constitutive relations can be constructed directly, taking samples with the same macroscopical state, and measuring in each one the incremental strain that results from the application of a specific stress increment. However, such tests are difficult to perform, because they require the fabrication of many samples with identical material properties.

Numerical simulations result as an alternative to the solution of this problem. They allow to develop different tests on identically generated samples and provide detailed information about micro-mechanical rearrangements during the loading process. Usually, disks or spheres are used to capture the granularity of the sample . Although the simplicity of their geometry allows to reduce the computer time of calculations, they do not provide a detailed description of realistic granular textures.

We present here a two dimensional discrete model that takes into account the diversity of shapes of the grains in the soils. The granular samples consist of randomly generated polygons. As presented by Tillemans , The interaction between the polygons could be handled letting the polygons interpenetrate each other and calculating the force as a function of their overlap. This approach has been successfully applied to model different processes, like fragmentation  and strain localization . A suitable contact force law is introduced in Sect.II A, that attempts to combine the Hertz contact law with the Coulomb friction criterion for quasi-static deformation.

The incremental stress-strain relation was calculated performing different stress increments on the same sample, and measuring the corresponding strain response. The stress is applied on the boundary through a flexible membrane that surrounds the sample. The modeling of such a membrane, whose details are given in Sect. II B, results more complex than rigid walls, However, it results more advantageous than walls, because it allows to implement a stress-controlled condition without any restriction in the deformation of the boundary.

## II Model

The polygons of this model are generated using a simple version of the Voronoi tessellation: First we set a random point in each cell of a regular square lattice. Then each polygon is constructed assigning to each point that part of the plane that is nearer to it than to any other point. Each polygon is subjected to interparticle contact forces and boundary forces that are inserted in Newton’s equation of motion.

### A Contact force

Usually, the interaction between two solid bodies in contact is described by a force applied on the flattened contact surface between them. Given two polygons in contact, such surface is obtained from the geometrical construction shown in Fig. 1. The points $`C_1`$ and $`C_2`$ result from the intersection between the edges of the polygons. The contact surface is taken as the segment that lies between those points. The vector $`\stackrel{}{S}=\stackrel{}{c_1c_2}`$ defines an intrinsic coordinate system at the contact $`(\widehat{t},\widehat{n})`$, where $`\widehat{t}=\stackrel{}{S}/|\stackrel{}{S}|`$ and $`\widehat{n}`$ is perpendicular to it. The deformation length is given by $`\delta =a/|\stackrel{}{S}|`$ where $`a`$ is the overlap area between the polygons. $`\stackrel{}{\mathrm{}}`$ is the branch vector, that connects the center of mass of the polygon to the point of application of the contact force, that is supposed to be the center of mass of the overlap area.

The normal elastic force is taken proportional to the deformation length as $`f_n^e=k_n\delta `$; the tangential force is calculated from the simplified Coulomb friction law, with a single friction coefficient $`\mu _s=\mu _d=\mu `$. Here $`\mu _s`$ is the static and $`\mu _d`$ the dynamic friction coefficient. This tangential force is implemented by an elastic spring $`f_t^e=k_t\xi `$, where $`\xi `$ grows linearly with the tangential displacement of the contact, whenever $`|f_t^e|<\mu f_n^e`$. We used the straightforward calculation of $`\xi `$ proposed by Brendel :

$`\xi (t)={\displaystyle _0^t}v_n(t^{})\mathrm{\Theta }(f_n^e(t^{})\mu |f_t^e(t^{})|)𝑑t^{}`$ (1)
where $`\mathrm{\Theta }`$ is the Heaviside function and $`\stackrel{}{v}`$ is the relative velocity at the contact, that depends on the linear velocity $`\stackrel{}{v}_i`$ and angular velocity $`\stackrel{}{\omega }_i`$ of the particles in contact according to:

$$\stackrel{}{v}=\stackrel{}{v}_i\stackrel{}{v}_j\stackrel{}{\omega }_i\times \stackrel{}{\mathrm{}}_i+\stackrel{}{\omega }_j\times \stackrel{}{\mathrm{}}_j$$
(2)
### B Boundary Forces

Let us now discuss how to apply the stress on the sample. One way to do that, would be to apply a perpendicular force on each edge of the polygons belonging to the external contour of the sample. Actually this does not work, because this force will act on all the fjords of the boundary. It produces an uncontrollable growth of cracks that with time ends up destroying the sample. Thus, results necessary to introduce a flexible membrane in order to restrict the boundary points that are subjected to the external stress.

The algorithm to identify the boundary is rather simple. The lowest vertex $`p`$ from all the polygons of the sample is chosen as the first point of the boundary list $`b_1`$. In Fig. 2 $`P`$ is the polygon that contains $`p`$, and $`qPQ`$ is the first intersection point between the polygons $`P`$ and $`Q`$ in counterclockwise orientation with respect to $`p`$. Starting from $`p`$, the vertices of $`P`$ in counterclockwise orientation are included in the boundary list until $`q`$ is reached. Next, $`q`$ is included in the boundary list. Then, the vertices of $`Q`$ between $`q`$ and the next intersection point $`rQR`$ in the counterclockwise orientation are included into the list. The same procedure is applied until one reaches the lowest vertex $`p`$ again. This is a very fast algorithm, because it only makes use of the contact points between the polygons, which are previously calculated to obtain the contact force.

The set of points that are in contact with the membrane are selected using a recursive algorithm. It is initialized with the vertices of the smallest convex polygon that encloses the boundary. (see Fig. 3). The lowest point of the boundary is selected as the first vertex of the polygon $`m_1=b_1`$. The second one $`m_2`$ is the boundary point $`b_i`$ that minimizes the angle $`\mathrm{}(\stackrel{}{b_1b_i})`$ with respect to the horizontal. The third one $`m_3`$ is the boundary point $`b_i`$ such that the angle $`\mathrm{}(\stackrel{}{m_2b_i},\stackrel{}{m_1m_2})`$ is minimal. The algorithm is recursively applied until the lowest vertex $`m_1`$ is reached again.

The points of the boundary are iteratively included in the list $`m_i`$ using the bending criterion proposed by Åstrøm : For each pair of consecutive vertices of the membrane $`m_i=b_i`$ and $`m_{i+1}=b_j`$ we choose that point from the subset $`\{b_k\}_{ikj}`$ that maximizes the bending angle $`\theta _b=\mathrm{}(\stackrel{}{b_kb_i},\stackrel{}{b_kb_j})`$. This point is included into the list, whenever $`\theta _b\theta _{th}`$. Here $`\theta _{th}`$ is a threshold angle for bending. This algorithm is repeatedly applied until there are not more points satisfying such bending condition.

The final result gives a set of segments $`\{\stackrel{}{m_im}_{i+1}\}`$ lying on the boundary of the sample. In order to apply the boundary forces, those segments are divided into two groups: A-type segments are those that coincide with an edge of a boundary polygon; B-type segments connect the vertices of two different boundary polygons.

0n each segment of the membrane $`\stackrel{}{m_im}_{i+1}`$ a force $`f_i=\sigma _iN_i`$ is applied, where $`\sigma _i`$ is the local stress and $`N_i`$ is the $`90^0`$ counterclockwise rotation of $`\stackrel{}{m_im}_{i+1}`$. This force is transmitted to the polygons in contact with it: if the segment is A-type, this force is applied in its midpoint; If the segment is B-type, half of the force is applied at each one of the vertices connected by this segment.

### C Molecular dynamic simulation

Before we implement the numerical solution of Newton’s equations it is convenient to make a dimensional analysis of the parameters. In such way we can keep the scale invariance of the model and reduce the parameter to a minimum of adimensional constants. All the polygons are supposed to have the same density. The mass $`m_i`$ of each polygon is measured in units of the mean mass $`m_0`$ of the Voronoi tessellation. The time is measured in fractions of the total loading time $`t_0`$. The evolution of the position $`\stackrel{}{x}_i`$ and the orientation $`\phi _i`$ of the $`ith`$ polygon is governed by the equations of motion:

$`\lambda ^2m_i\ddot{\stackrel{}{x}}_i+{\displaystyle \underset{c}{}}\stackrel{}{f_i^c}+{\displaystyle \underset{c_b}{}}{\displaystyle \frac{\sigma _i^b}{k_n}}\stackrel{}{f}_i^b`$ $`=`$ $`0`$ (3)
$`\lambda ^2I_i\ddot{\phi }_i+{\displaystyle \underset{c}{}}\stackrel{}{\mathrm{}}_i^c\times \stackrel{}{f_i^c}+{\displaystyle \underset{c_b}{}}{\displaystyle \frac{\sigma _i^b}{k_n}}\stackrel{}{\mathrm{}}_i^b\times \stackrel{}{f}_i^b`$ $`=`$ $`0`$ (4)
The sums go over all those particles and boundary segments that are in contact with the $`ith`$ polygon. The interparticle contact forces $`\stackrel{}{f_i^c}`$ and boundary forces $`\stackrel{}{f_i^b}`$ are given by:

$`\stackrel{}{f_i^c}`$ $`=`$ $`(\delta _i^c+\lambda \gamma mv_n^c)\widehat{n}_i^c+\zeta (\xi _i^c\lambda \gamma mv_t^c)\widehat{t}_i^c`$ (5)
$`\stackrel{}{f_i^b}`$ $`=`$ $`N_i^b\lambda \gamma m_i\stackrel{}{v_i}`$ (6)
Here $`\delta _i^c`$ and $`\xi _i^c`$ are the deformation length and the tangential displacement of the contact, which were defined in Sect. II A; $`\sigma _i^b`$ is the stress applied on the boundary segment $`T_i^b`$, defined in Sect. II B. Artificial viscous terms must be included in Eq. (6) to keep the stability of the numerical solution and reduce the acoustic waves generated during the loading process. $`\stackrel{}{v}^c`$ is the relative velocity at the contact (Eq. (2)) and $`m=(1/m_i+1/m_j)^1`$ is the effective mass of the two polygons in contact.

There are four microscopic parameters in the model: the viscosity $`\gamma `$, the ratio $`\lambda =t_s/t_o`$ between the characteristic period of oscillation $`t_s=\sqrt{k_n/m_0}`$ and the loading time $`t_0`$, the friction coefficient $`\mu `$, and the ratio $`\zeta =k_t/k_n`$ between the tangential $`k_t`$ and normal $`k_n`$ stiffness of the interparticle contacts.

The viscosity factor $`\gamma `$ is related to the normal restitution coefficient . It was taken large enough to have a high dissipation, but not too large to keep the numerical stability of the method. The ratio $`\lambda `$ was chosen small enough in order to to avoid rate-dependence in the strain response, as corresponds to the quasi-static approximation. Technically, that is done by looking for the value of $`\lambda `$ such that a reduction of it by half makes a change of the strain response less than $`5\%`$.

The two parameters $`\zeta `$ and $`\mu `$ determine the constitutive response of the system. For example, the micro-mechanical analysis of the strain response shows that the Young’s modulus and Poisson’s ratio depend on $`\zeta `$ . In the other hand, $`\mu `$ can be directly related to the friction angle of the material . Although the study of the dependence of the constitutive response on those parameters is an important point, such quantities have been kept fixed in this work.

The boundary conditions yield more dimensional parameters. The initial hight $`H_0`$ and width $`W_0`$ of the sample, and the characteristic length $`\mathrm{}_0`$ of the polygons define two geometrical parameters, which are the shape ratio $`W_0/H_0`$ and the granularity $`\mathrm{}_0/H_0`$ of the sample.

In order to keep overlaps much smaller than the characteristic area of the polygons. The ratio $`\sigma _i/k_n`$ between the stress applied on the membrane and the stiffness of the contacts is restricted to small values. This was implemented by fixing the contact stiffness to a value closed to the experimental granular stiffness $`k_n=160MPa`$. Then the stress is chosen in such way that it does not excess $`1\%`$ of this value.

## III Stress-strain Calculation

### A Theorical background

The macroscopic state of the system is characterized by the stress tensor and the void radio $`e`$. The area fraction of voids in the sample defines the void ratio. Initially $`e_0=0`$ due to the Voronoi tessellation used. The stress controlled test was restricted to stress states without-off diagonal components. The diagonal component, the axial $`\delta _1`$ and lateral $`\delta _3`$ stress, define the stress vector:

$$\stackrel{~}{\sigma }=\left[\begin{array}{c}p\\ q\end{array}\right]=\frac{1}{2}\left[\begin{array}{c}\delta _1+\delta _3\\ \delta _1\delta _3\end{array}\right]$$
(7)
where $`p`$ and $`q`$ are the pressure and the shear stress. The domain of admissible stresses is bounded by the failure surface. When the system reaches this surface it becomes unstable and fails.

Before failure, the constitutive behavior can be obtained performing small changes in the stress and evaluating the resultant deformation. An infinitesimal change of the stress vector $`d\stackrel{~}{\sigma }`$ produces an infinitesimal deformation of the sample, which is given by a change of height $`dH`$ and width $`dW`$. This defines the axial strain $`dϵ_1=dH/H`$ and lateral strain $`dϵ_3=dW/W`$ increments. The volumetric strain $`dϵ_v`$ and the shear strain $`dϵ_\gamma `$ increments define the incremental strain vector:

$$d\stackrel{~}{ϵ}=\left[\begin{array}{c}dϵ_v\\ dϵ_\gamma \end{array}\right]=\left[\begin{array}{c}dϵ_1+dϵ_3\\ dϵ_1dϵ_3\end{array}\right]$$
(8)
Each state of the sample is related to a single point in the stress space, and the quasi-static evolution of the system is represented by the movement of this point in the stress space. The constitutive relation is formulated taking the incremental strain as a function of the incremental stress and the stress state.

$$d\stackrel{~}{ϵ}=(d\stackrel{~}{\sigma },\stackrel{~}{\sigma })$$
(9)
If there is no rate - dependence in the constitutive equation, $`(d\stackrel{~}{\sigma })`$ is an homogeneous function of degree one. In this case, the application of the Euler identity  shows that Eq. (9) can be reduced to.

$$d\stackrel{~}{ϵ}=M(\widehat{\theta },\stackrel{~}{\sigma })d\stackrel{~}{\sigma }$$
(10)
Where $`\widehat{\theta }`$ is the unitary vector defining a specific direction in the stress space:

$$\widehat{\theta }=\frac{d\stackrel{~}{\sigma }}{|d\stackrel{~}{\sigma }|}\left[\begin{array}{c}\mathrm{cos}\theta \\ \mathrm{sin}\theta \end{array}\right],|d\stackrel{~}{\sigma }|=\sqrt{dp^2+dq^2}$$
(11)
The constitutive relation results from the calculation of $`d\stackrel{~}{ϵ}(\theta )`$, where each value of $`\theta `$ is related to a particular mode of loading. Some special modes are listed in the table:

The relation (10) has been proposed by Darve  , and it contains all the possible constitutive equations. In order to interpret our particular results, it is convenient to make some approximations: First, if the load increments are taken small enough, the tensor $`M(\theta )`$ can be supposed to be lineal in each stress direction. Then, we assume that the strain can be separated in an elastic (recoverable) and a plastic (unrecoverable) component:

$$d\stackrel{~}{ϵ}=d\stackrel{~}{ϵ}^e+d\stackrel{~}{ϵ}^p$$
(12)
$$d\stackrel{~}{ϵ}^e=D(\stackrel{~}{\sigma })d\stackrel{~}{\sigma }$$
(13)
$$d\stackrel{~}{ϵ}^p=J(\theta ,\stackrel{~}{\sigma })d\stackrel{~}{\sigma }$$
(14)
Here, $`D^1`$ defines the stiffness tensor, and $`J=MD`$ the flow rule of plasticity, which results from the calculation of $`d\stackrel{~}{ϵ}^e(\theta )`$ and $`d\stackrel{~}{ϵ}^p(\theta )`$.

### B The method

A numerical method has been implemented in order to find the elastic $`d\stackrel{~}{ϵ}^e`$ and plastic $`d\stackrel{~}{ϵ}^p`$ components of the strain as function of the stress state $`\stackrel{~}{\sigma }`$ and the stress direction $`\theta `$. Fig. 5 shows the three steps of the procedure:

1) The sample is driven to the stress state $`\stackrel{~}{\sigma }`$. First, it is isotropically compressed until it reaches the stress value $`\delta _1=\delta _3=pq`$. Next, it is subjected to axial loading, in order to increase the axial stress $`\delta _1`$ until $`p+q`$ (see Fig. 4). When the stress state $`\stackrel{~}{\sigma }=[pq]^T`$ is reached, ($`A^T`$ being the transpose of $`A`$) the sample is allowed to relax.

2) Loading the sample from $`\stackrel{~}{\sigma }`$ to $`\stackrel{~}{\sigma }+d\stackrel{~}{\sigma }`$ the strain increment $`d\stackrel{~}{ϵ}`$ is obtained. This procedure is implemented choosing different stress directions according to Eq. (11). Here the stress modulus is fixed to $`|d\stackrel{~}{\sigma }|=10^4p`$.

3) The sample is unloaded until the original stress state $`\stackrel{~}{\sigma }`$ is reached. Then one finds a remaining strain $`d\stackrel{~}{ϵ}^p`$, that corresponds to the plastic component of the incremental strain. Since the stress increments are taken small enough, the unloaded stress-strain path is practically elastic. Thus, the difference $`d\stackrel{~}{ϵ}^e=d\stackrel{~}{ϵ}d\stackrel{~}{ϵ}^p`$ represents the elastic component of the strain.

One could be concerned about the dependence of the strain response on the way how the stress state is reached. We found that there is not remarkable dependence of the strain response on the stress path, whenever the stress components are quasi - static and monotonically increased. Otherwise, a strong reduction in the plastic component of the strain is observed. In fact, when the plastic response is calculated after the sample is unloaded, the plasticity results smaller than that one calculated after a monotonic load. Furthermore, there is no plastic component in the strain response when elastic waves are previously generated in the sample. Those memory effects suggest that the plastic component of the strain depend on the history of the deformation, and is kept unchanged only if the sample is subjected to quasi-static and monotonic load.

Fig. 6 shows the load-unload paths and the corresponding strain response. They were taken from a stress-state with $`q=0.5p`$. The end of the load paths in the stress space map into a strain envelope response $`d\stackrel{~}{ϵ}(\theta )`$ in the strain space. Likewise, the end of the unload paths map into a plastic envelope response $`d\stackrel{~}{ϵ}^p(\theta )`$. The yield direction $`\varphi `$ can be found from this response, as the direction in the stress space where the plastic response is maximal. The flow rule can be obtained taking the direction $`\psi `$ of the maximal plastic response in the strain space. These angles do not agree, that reveals the necessity to analyze this behavior in the framework of the non-associated theory of plasticity (see Sect. IV C).

## IV Constitutive Relation

Fig.7 summarizes the global elasto-plastic behavior. The elastic response, calculated from Eq. (12), has a centered ellipse as envelope response. This can be related to the micro-contact structure using a local linear relation in each point of the stress space (see Sect.IV B). The solid line represents the failure surface, which separates the stable states of the unstable ones (see Sect. IV A). The plastic envelope response is almost on a straight line. The modulus and the orientation of this envelope depend on the stress state through a certain number of material parameters, which are given in Sect. IV C. All the quantities obtained in this section have been calculated from the average over five different samples of $`10\times 10`$ particles each one.

### A Failure surface

The failure line was calculated looking for the values of stress for which the system becomes unstable: for each pressure $`p`$, there is a critical shear stress $`q_c(p)`$, below which the sample reaches a stable state with an exponential decay of its kinetic energy. For shear stress values above the critical one, the sample develops an instability and fails. Fig. 8 shows the interface between these two stress states, that can be accurately fitted by the power law:

$$\frac{q_c}{p_0}=\mu ^{}\left(\frac{p}{p_0}\right)^\beta $$
(15)
Here $`p_0=1.0MPa`$ is the reference pressure, and $`\mu ^{}=0.78\pm 0.03`$ is the Mohr-Coulomb friction coefficient . The power law dependence on the pressure, with exponent $`\beta =0.92\pm 0.02`$ implies a significant deviation from the Mohr-Coulomb theory. Moreover, the empirical criteria of failure for most rocks  shows a power law dependence of the form of Eq. (15). It seems that additional features beyond the Mohr-Coulomb analysis are taking place when the sample fails, that will be discussed in sect. IV C.

### B stiffness

Hooke’s law of elasticity states that the stiffness tensor of isotropic materials can be written in terms of two material parameters: the Young’s modulus $`E`$ and the Poisson’s ratio $`\nu `$. However, the isotropy is not fulfilled when the stress state is far from the hydrostatic axis. Indeed, numerical simulations  and photo-elastic experiments  on granular materials show that the loading induces a significant departure from isotropy in the contact network.

The anisotropy of the granular sample can be characterized by the distribution of the micro-contact normal vectors $`\widehat{n}_i^c`$ (see Fig. 1). Our numerical simulations show that the structural changes of micro-contacts are principally due to the opening of contacts whose normal vectors are nearly aligned around the direction perpendicular to the load. Let us call $`N(\phi )\mathrm{\Delta }\phi `$ the number of contacts per particle, oriented between the angles $`\phi `$ and $`\phi +\mathrm{\Delta }\phi `$, measured with respect to the direction in which the sample is loaded. The lowest order of anisotropy can be described by:

$$N(\phi )=\frac{1}{2\pi }\left[N+(N_0N)\mathrm{cos}(2\phi )\right]$$
(16)
Here $`N`$ is the average coordination number of the polygons, whose initial value $`N_0=6.0`$ reduces as the load is increased. Fig. 9 shows this reduction. A critical line is found around $`q=0.12p`$, below which there are no structural changes in the contact network. Above this limit an induced anisotropy arises due to opened contacts whose amount follows a power law dependence.

In order to describe the effect of the anisotropy in the elastic response we proceed as follows: first, an additional parameter $`\alpha `$ is included in Hooke’s law

$$\left[\begin{array}{c}dϵ_1^e\\ dϵ_3^e\end{array}\right]=\frac{1}{E}\left[\begin{array}{cc}1\alpha & \nu \\ \nu & 1+\alpha \end{array}\right]\left[\begin{array}{c}d\sigma _1\\ d\sigma _3\end{array}\right]$$
(17)
Then, these three parameters are supposed to be depended on the internal damage parameter $`d`$:

$$d=\frac{N_0N}{N_0}$$
(18)
The tensor $`D`$ defined in Eq. (13) is calculated from Eq. (17) using the definition of the stress and strain vectors given in Eqs. (7) and (8). One obtains:

$$D=\frac{2}{E}\left[\begin{array}{cc}1\nu & \alpha \\ \alpha & 1+\nu \end{array}\right]$$
(19)
The diagonal components of this tensor are respectively the inverse of the bulk modulus and of the shear modulus. The non-diagonal component results from the anisotropy of the sample, and it couples the compression mode with the shearing deformation. These three variables are calculated from the elastic response $`d\stackrel{~}{ϵ}^e(\theta )`$ by the introduction of the following function:

$$R(\theta )=\frac{d\stackrel{~}{\sigma }^Td\stackrel{~}{ϵ}^e}{|d\stackrel{~}{\sigma }|^2}$$
(20)
by substitution of Eqs. (13) and (11) into Eq. (20), one sees that $`R`$ is the quadratic form of $`D`$:

$$R(\theta )=\widehat{\theta }^TD\widehat{\theta }=\frac{2}{E}\left[1\nu \mathrm{cos}(2\theta )\alpha \mathrm{sin}(2\theta )\right]$$
(21)
Using this equation, the components of $`D`$ can be evaluated as the Fourier coefficients of $`R`$:

$`{\displaystyle \frac{1}{E}}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _0^{2\pi }}R(\theta )𝑑\theta `$ (22)
$`\nu `$ $`=`$ $`{\displaystyle \frac{E}{2\pi }}{\displaystyle _0^{2\pi }}R(\theta )\mathrm{cos}(2\theta )𝑑\theta `$ (23)
$`\alpha `$ $`=`$ $`{\displaystyle \frac{E}{2\pi }}{\displaystyle _0^{2\pi }}R(\theta )\mathrm{sin}(2\theta )𝑑\theta `$ (24)
Figs. 10, 11 and 12 show the results of the calculation of the Young’s modulus $`E`$, the Poisson’s ratio $`\nu `$ and the anisotropy factor $`\alpha `$. Below the limit of isotropy, Hooke’s law can be applied: $`EE_0`$, $`\nu \nu _0`$ and $`\alpha 0`$. On the other hand, above the limit of isotropy a reduction of the Young’s modulus is found, along with an increase of the Poisson’s ratio and the anisotropy factor. The functional dependence of those parameters with the internal damage parameter $`d`$ are evaluated developing their Taylor’s series around $`d=0`$:

$`E(d)`$ $`=`$ $`E(0)+E^{}(0)d+O\left(d^2\right)`$ (25)
$`\alpha (d)`$ $`=`$ $`\alpha (0)+\alpha ^{}(0)d+O\left(d^2\right)`$ (26)
$`\nu (d)`$ $`=`$ $`\nu (0)+\nu ^{}(0)d+\nu ^{\prime \prime }(0)d^2+O\left(d^3\right)`$ (27)
The coefficients of this expansion are calculated from the best fitting of those expansions. Figs. 10 and 12 show that the linear approximation is good enough to reproduce the Young’s modulus and the anisotropy factor. The fit of the Poisson’s ratio, however, requires the inclusion of a quadratic approximation, implying that it has a non-linear dependence on the damage parameter (Fig. 11).

### C Plastic Flow

The formulation of the non-associated theory of plasticity requires the evaluation of three material functions: the yield direction $`\varphi `$, the flow direction $`\psi `$, and the plastic modulus $`h`$. These quantities can be calculated from the plastic response $`d\stackrel{~}{ϵ}^p(\theta )`$, as follows:

The yield direction is given by the incremental stress direction $`\varphi `$ with maximal plastic response:

$$|d\stackrel{~}{ϵ}^p(\varphi )|=\underset{\theta }{\mathrm{max}}|d\stackrel{~}{ϵ}^p(\theta )|$$
(28)
The flow direction is defined from the orientation of the plastic response at its maximum value:

$$\psi =\mathrm{arctan}(\frac{dϵ_\gamma ^p}{dϵ_v^p})|_{\theta =\varphi }$$
(29)
The plastic modulus is obtained form the modulus of the maximal plastic response:

$$\frac{1}{h}=\frac{|d\stackrel{~}{ϵ}^p(\varphi )|}{|d\stackrel{~}{\sigma }|}$$
(30)
Reciprocally, the plastic response can be expressed in terms of these quantities. Let us define the unitary vectors $`\widehat{\psi }`$ and $`\widehat{\psi }^{}`$. The first one is oriented in the direction of $`\psi `$ and the second one is the $`90^o`$ rotation of $`\widehat{\psi }`$. The plastic strain is written in this basis as:

$`d\stackrel{~}{ϵ}^p(\theta )`$ $`=`$ $`[(d\stackrel{~}{ϵ}^p)^T\widehat{\psi }]\widehat{\psi }+[(d\stackrel{~}{ϵ}^p)^T\widehat{\psi }^{}]\widehat{\psi }^{}`$ (31)
$``$ $`{\displaystyle \frac{1}{h}}\left[f(\theta )\widehat{\psi }+g(\theta )\widehat{\psi }^{}\right]`$ (32)
The plastic profiles $`f(\theta )`$ and $`g(\theta )`$ are shown in Fig. 13. The first one is approximately the same for all the stress states, and can be well fitted to a cosine function, centered on the yield direction $`\varphi `$ and truncated to zero for the negative values. The last profile depends on the stress value, and is difficult to evaluate, because it is of the same order as the statistical fluctuations. However, the contribution of $`g`$ to the total strain response is negligible. In order to simplify the description of the plastic response, the following approximation is made:

$$g(\theta )f(\theta )[[cos(\theta \varphi )]]=[[\widehat{\varphi }^T\widehat{\theta }]]$$
(33)
where $`[[]]`$ defines the function

$$[[x]]=\{\begin{array}{cc}\hfill 0:& x0\hfill \\ \hfill x:& x>0\hfill \end{array}$$
(34)
Now, The flow rule results from the substitution of Eqs.(32) and (33) into Eq. (14):

$$J(\theta )d\stackrel{~}{\sigma }=\frac{[[\widehat{\varphi }^Td\stackrel{~}{\sigma }]]}{h}\widehat{\psi }$$
(35)
The yield direction and the flow direction have been calculated for different stress states. The results are shown in Fig. 14. Both angles are quite different, which is a clear deviation from Drucker’s normality postulate . Indeed, many experimental results on soil deformation  have confirmed that these angles are completely different. Thus Drucker’s postulate is not fulfilled in the deformation of granular materials, and the main reason for that is the rearrangement of contacts on small deformations which are not taken into account in this theory. On the other hand, all the sliding, opening, and other micro-mechanical rearrangements can be well handled in the discrete element formulation, which results more adequate to describe the soil deformation.

The material constants are evaluated from the dependence of the plastic quantities on the stress: the yield direction and the flow direction can be roughly approximated by straight lines:

$`\varphi `$ $`=`$ $`\varphi _0+\varphi _0^{}{\displaystyle \frac{q}{p}}`$ (36)
$`\psi `$ $`=`$ $`\psi _0+\psi _0^{}{\displaystyle \frac{q}{p}}`$ (37)
The four material parameters $`\varphi _0=46^o\pm 0.75^o`$, $`\varphi _0^{}=88.3^o\pm 0.6^o`$, $`\psi _0=78.9^o\pm 0.2^o`$ and $`\psi _0^{}=59.1^o\pm 0.4^o`$ are obtained from the linear fit of the data. On the other hand, Fig. 15 shows that the plastic modulus depends on the stress through a power law relation:

$$h=h_0\left[1\frac{q}{\mu ^{}p_0}(\frac{p_0}{p})^\vartheta \right]^\eta $$
(38)
There are four additional material parameters: The plastic modulus $`h_0=14.5\pm 0.05`$ at $`q=0`$, the Mohr-Coulomb friction coefficient $`\mu ^{}`$ (see Eq. (15)), and the exponents $`\eta =2.7\pm 0.04`$ and $`\vartheta =0.981\pm 0.002`$ .

The plastic limit surface is given by the stress states where the plastic deformation becomes infinite. According to the flow rule ( Eq. (35)), it is found looking for the stress values where Eq. (38) vanishes:

$$\frac{q_p}{p_0}=\mu ^{}\left(\frac{p}{p_0}\right)^\vartheta $$
(39)
It is important to point out that the failure surface –given in Eq. (15)– does not correspond to the plastic limit surface. Actually, This matter has already been discussed in the framework of the Hill’s condition of instability  the bifurcation analysis , which predict that the instability should be reached strictly inside the plastic limit surface.

## V CONCLUDING REMARKS

The elasto-plastic response of a Voronoi tessellated sample of polygons has been calculated in the case of monotonic and quasi-static loading. It can be written in a simple form as:

$$d\stackrel{~}{ϵ}=D(d)d\stackrel{~}{\sigma }+\frac{[[\widehat{\varphi }^Td\stackrel{~}{\sigma }]]}{h}\widehat{\psi }$$
(40)
The plastic response reflects the non-associated features of realistic soils. Here the yield direction and flow direction are linearly related to the ratio $`q/p`$, and the plastic modulus obeys a power law relation with a weak pressure dependence.

The classical parameters of elasticity – the Young’s modulus and the Poisson’s ratio – are not material constants, because they depend on the internal damage parameter. Therefore, Eq. (40) is not complete, and it is necessary to include a constitutive equation that relates the internal damage to the external load. By focusing on the details of the dynamic of the micro-contacts, significant progress may be made in the macroscopic description of the deformation.

The elasto-plastic response leads to the identification of three different regimes which are shown in Fig. 16. Zone I corresponds to the isotropic regime, characterized by small plastic deformations and a linear elastic regime. In the zone II open contacts are detected, which must be taken into account in the calculation of the non linear elasticity. Zone III corresponds to unstable states so that the stress-strain relation can not be calculated here. The extrapolation of the strain response in this region shows that the plastic strain must have a finite value just before the instability is reached.

The above observation leads to the open question of the nature of the failure . Numerical simulations on strain controlled tests show that strain localization is the most typical mode of failure. The fact that it appears before the sample reaches the plastic limit surface suggests that the appearance of the instability is not completely determined by the macroscopic state.

The role of the microstructure on the strain localization has been intensely studied in the last years. Future work is the creation of samples with different granular textures –for example, changing the void ratio distributions and the polydispersity of the grains–. Then we can deal with the question of how a change in the microstructure affects the elasto-plastic response and the strain localization.

## Acknowledgments

We thank F. Darve, P. Vermeer, F. Kun, J. Åstrøm and S. Luding for helpful discussions and acknowledge the support of the Deutsche Forschungsgemeinschaft within the research group Modellierung kohäsiver Reibungsmaterialen.


format: text

text:
# Slow oscillations of magnetoresistance in quasi-two-dimensional metals

## Abstract

Slow oscillations of the interlayer magnetoresistance observed in the layered organic metal $`\beta `$-(BEDT-TTF)2IBr2 are shown to originate from the slight warping of its Fermi surface rather than from independent small cyclotron orbits. Unlike the usual Shubnikov-de Haas effect, these oscillations are not affected by the temperature smearing of the Fermi distribution and can therefore become dominant at high enough temperatures. We suggest that the slow oscillations are a general feature of clean quasi-two-dimensional metals and discuss possible applications of the phenomenon.

Quantum oscillations of magnetization (de Haas-van Alphen, dHvA, effect) and magnetoresistance (Shubnikov-de Haas, SdH, effect) in strong magnetic fields have been successfully used for many years in order to investigate Fermi surfaces and other electronic properties of metals Shoe . In three-dimensional (3D) metals the oscillations are determined by only a small fraction of conduction electrons near extremal cyclotron orbits on the Fermi surface (FS). Therefore they are very weak and can perfectly be described in terms of the so-called Lifshitz-Kosevich (LK) formalism Shoe . Recently, however, materials of lower dimensionality such as two-dimensional electron gas systems, layered metal oxides and organic conductors have been of high interest. In these materials the relative number of electrons contributing to the oscillations is much higher than in 3D metals and the main assumptions of the standard theory, i.e. a weak oscillation amplitude and constant chemical potential, are often no longer fulfilled.

Some crystalline organic metals, being very clean and having extremely large anisotropy of electronic properties, are excellent objects for studying specific features of the magnetic quantum oscillations in the low-dimensional limit Jero ; Harr ; Wosn1 ; Balt ; PhSh . Various deviations of the oscillatory magnetization from the LK theory observed in these compounds can be fairly well explained by a semi-phenomenological description of the quasi-two-dimensional (Q2D) dHvA effect developed in a number of theoretical works Harr ; Gv1 ; dHvADiff ; Pavel2 . The situation with the SdH effect is less satisfactory. Despite some progress in understanding several anomalies of the oscillatory magnetoresistance Harr ; PhSh ; ShubDiff ; Av  a consistent theory is still lacking. Moreover, there are still open qualitative questions. One of these questions concerns the origin of slow oscillations of magnetoresistance which were first found in the Q2D organic metal $`\beta `$-(BEDT-TTF)2IBr2 Kar1 ; Kar2 . Similar oscillations have been observed in other layered organic conductors, e.g. $`\beta `$-(BEDT-TTF)2I3 Kar3 ; Wosn2 , $`\kappa `$-(BEDT-TTF)2Cu2(CN)3 Ohmi , and $`\kappa `$-(BEDT-TSF)2C(CN)3 togo . Since the behavior of these oscillations strongly resembles that of the SdH effect, they have been supposed to originate from additional very small pockets of the FS. However, the band structure calculations which are basically believed to correctly reproduce the FS topology of organic metals (see e.g. OMRev  for a review) do not explain such small groups of carriers in any of the cited compounds.

In this Letter, we report on detailed studies of the oscillating interlayer magnetoresistance and magnetization of $`\beta `$-(BEDT-TTF)2IBr2 at various orientations of magnetic field. Our results provide an unequivocal evidence that the slow oscillations of the magnetoresistance do not reveal any new carriers but are ultimately caused by a weak warping of the single cylindrical FS in this Q2D metal. We propose a theoretical explanation of the phenomenon which appears to be in good agreement with the experiment.

The experiment was performed on a high-quality $`[R(290\text{K})/R(2\text{K})3000`$\] single crystal with the dimensions $`0.6\times 0.3\times 0.12\text{mm}^3`$. The sample was mounted into a measuring cell in a 3He-cryostat allowing simultaneous registration of the resistance and magnetic torque Weis  at different orientations of the magnetic field produced by a 14 T superconducting magnet. The field orientation was defined by the angle $`\theta `$ between the field direction and the normal to the highly conducting $`ab`$-plane.

The FS of $`\beta `$-(BEDT-TTF)2IBr2 is a slightly warped cylinder with the axis along the least conducting direction Kar2 ; Wosn3 ; Kar4 . The general behavior of the oscillating interlayer magnetoresistance and magnetization is shown in Fig. 1 and fully agrees with the previous reports Kar2 ; Wosn3 . Cyclotron orbits on the cylindrical FS give rise to rapid SdH (Fig. 1a) and dHvA (Fig. 1b) oscillations with the frequency $`F=3930`$ T that corresponds to the FS cross-sectional area $`S0.53S_{\text{B}Z}`$ ($`S_{\text{B}Z}`$ is the Brillouin zone area). The amplitudes of both the SdH and dHvA oscillations are clearly modulated due to a slight warping of the FS. A comparison of the beat and fundamental frequencies yields the evaluation of the warping: $`\mathrm{\Delta }S/S=2F_{\text{beat}}/F10^2`$.

In addition to the rapid SdH oscillations, the magnetoresistance exhibits prominent slow oscillations with the frequency $`F_{\text{slow}}42`$ T. Due to a large cyclotron mass, $`m4.2m_e`$ ($`m_e`$ is the free electron mass), the fundamental oscillations rapidly diminish with increasing temperature and can barely be resolved at the highest field at $`T=1.4`$ K. By contrast, the amplitude of the slow oscillations remains almost the same as at 0.6 K. Noteworthy, no trace of slow oscillations has been found in our magnetization measurements.

A clue for understanding the origin of the slow oscillations is the dependence of their frequency on the tilt angle $`\theta `$ displayed in Fig. 2a comment0 . Unlike the rapid oscillations having the $`1/\mathrm{cos}\theta `$-dependence typical of Q2D metals, the slow oscillations show a strong non-monotonic change of their frequency with $`\theta `$. Such a behavior immediately reminds that of the beat frequency $`F_{\text{beat}}(\theta )`$. At certain magnetic field orientations, repeated periodically with $`\mathrm{tan}\theta `$, the areas of cyclotron orbits in a Q2D metal become independent of the orbit positions in $`k`$-space Yam . This obviously results in periodic drops of $`F_{\text{beat}}`$ to zero Wosn3 . At the same angles the semiclassical part of the magnetoresistance turns out to sharply increase AMRO , giving rise to the so-called angle-dependent magnetoresistance oscillations (AMRO) Kar2 ; Kar4 . Fig. 2b shows the angular dependence of the background magnetoresistance of our sample at $`B=14`$ T, revealing two prominent AMRO peaks at $`\theta 33^{}`$ and $`20^{}`$. It clearly correlates with the angular dependence of $`F_{\text{slow}}`$: the latter rapidly decreases as the magnetoresistance approaches the AMRO peaks.

It is of course very tempting to directly compare $`F_{\text{slow}}`$ and $`F_{\text{beat}}`$. Unfortunately, the rapid oscillations were generally observed in a relatively small field interval in our experiment, so that the beat period could not be reliably measured in a sufficiently wide angular range. Nevertheless, estimations made at a few angles (see e.g. the data in Fig. 1) reveal the relationship $`F_{\text{slow}}(\theta )2F_{\text{beat}}(\theta )`$ within the experimental error bar. Thus, one can conclude that the slow oscillations and beats of the rapid oscillations have the same physical origin, i.e. both are directly related to the warping of the cylindrical FS.

In order to clarify the mechanism responsible for the slow oscillations, we first note that the interlayer conductivity contains several factors which, in general, oscillate in magnetic field with the frequency $`F`$ determined by the cross-sectional area of the FS cylinder. The amplitudes of the oscillations are modulated, due to the warping of the cylinder, with the frequency $`F_{\text{beat}}=(2t_{}/ϵ_F)FF`$ ($`t_{}`$ is the interlayer transfer integral and $`ϵ_F`$ the Fermi energy). The product of two oscillating quantities with modulated amplitudes $`\stackrel{~}{\alpha }`$ and $`\stackrel{~}{\beta }`$ yields a slowly oscillating term, e.g. $`(1+\stackrel{~}{\alpha }\mathrm{cos}x)(1+\stackrel{~}{\beta }\mathrm{cos}x)=1+(\stackrel{~}{\alpha }+\stackrel{~}{\beta })\mathrm{cos}x+(\stackrel{~}{\alpha }\stackrel{~}{\beta }/2)\mathrm{cos}2x+\stackrel{~}{\alpha }\stackrel{~}{\beta }/2`$. Here, the last term describes slow oscillations.

To be more explicit, we consider the interlayer conductivity as determined from the Boltzmann transport equation Mah :

$$\sigma _{zz}=e^2𝑑ϵ\left(n_F^{}\left(ϵ\right)\right)I(ϵ)\tau (ϵ),$$
(1)
where $`n_F^{}\left(ϵ\right)=1/\{4T\mathrm{cosh}^2\left[(ϵ\mu (B))/2T\right]\}`$ is the derivative of the Fermi distribution function, $`I(ϵ)\left|v_z(ϵ)\right|^2`$ the square of the electron interlayer velocity $`v_z`$ summed over all states at the energy $`ϵ,`$ and $`\tau (ϵ)`$ is the momentum relaxation time at the energy $`ϵ`$. In Born approximation, $`\tau (ϵ)`$ is inversely proportional to the density of states (DoS): $`\tau (ϵ)\rho ^1(ϵ)`$ and oscillates in a magnetic field Gv1 :

$$\tau (ϵ)\left[1+2\underset{p=1}{\overset{\mathrm{}}{}}(1)^p\mathrm{cos}\left(\frac{2\pi pϵ}{\mathrm{}\omega _c}\right)J_0\left(\frac{4\pi pt_{}}{\mathrm{}\omega _c}\right)R_D\right]^1,$$
(2)
where $`\omega _c=eB/m`$ is the cyclotron frequency, and $`R_D=\mathrm{exp}(2\pi ^2pk_BT_D/\mathrm{}\omega _c)`$ is the usual scattering damping factor Shoe  ($`T_D`$ is the Dingle temperature). If the cyclotron energy is comparable to the interlayer transfer integral, the oscillations of the quantity $`I(ϵ)`$ become also important. They are given by PhSh

$$I(ϵ)1+\frac{\mathrm{}\omega _c}{\pi t_{}}\underset{p=1}{\overset{\mathrm{}}{}}\frac{\left(1\right)^p}{p}\mathrm{cos}\left(\frac{2\pi pϵ}{\mathrm{}\omega _c}\right)J_1\left(\frac{4\pi pt_{}}{\mathrm{}\omega _c}\right)R_D.$$
(3)
The beats of the oscillations of $`\tau `$ and $`I`$ are given by the 0-th and 1-st order Bessel functions $`J_0`$ and $`J_1`$ in Eqs. (2) and (3), respectively. Their product eventually gives rise to a slowly oscillating term in the conductivity. Since the FS warping is large enough in our case, $`4\pi t_{}>\mathrm{}\omega _c`$, one can approximate the Bessel functions as $`J_0(z)=\sqrt{2/\pi z}\mathrm{cos}(z\pi /4)`$ and $`J_1(z)=\sqrt{2/\pi z}\mathrm{sin}(z\pi /4)`$. Further, taking into account the weak amplitude of the SdH effect (see Fig. 1) and strong harmonic damping (the second harmonic never exceeded 1% of the fundamental one in our experiment), we neglect oscillations of the chemical potential $`\mu `$ and restrict our consideration to the lowest order in the damping factors. Then, substituting Eqs. (2) and (3) into Eq. (1) and performing the integration over energy at finite temperature, we obtain:

$`\sigma _{zz}`$ $`=`$ $`\sigma _0\{1+2\sqrt{{\displaystyle \frac{\mathrm{}\omega _c}{2\pi ^2t}}}\sqrt{1+a^2}\times `$
$`\times \mathrm{cos}\left({\displaystyle \frac{2\pi \mu }{\mathrm{}\omega _c}}\right)\mathrm{cos}\left({\displaystyle \frac{4\pi t_{}}{\mathrm{}\omega _c}}{\displaystyle \frac{\pi }{4}}+\varphi \right)R_DR_T+`$
$`+{\displaystyle \frac{\mathrm{}\omega _c}{2\pi ^2t_{}}}[1+\sqrt{1+a^2}\mathrm{cos}\left(2[{\displaystyle \frac{4\pi t_{}}{\mathrm{}\omega _c}}{\displaystyle \frac{\pi }{4}}+{\displaystyle \frac{\varphi }{2}}]\right)]R_{D}^{}{}_{}{}^{2}\}`$
where $`\varphi =\mathrm{arctan}(a)`$ and $`a=\mathrm{}\omega _c/2\pi t_{}`$. $`R_T=\left(2\pi ^2k_BT/\mathrm{}\omega _c\right)/\mathrm{sinh}\left(2\pi ^2k_BT/\mathrm{}\omega _c\right)`$ is the temperature damping factor which is equal to that in the LK theory and comes from the integration with the Fermi distribution function. The scattering factor $`R_D^{}`$ has the same form as $`R_D`$ but includes a different Dingle temperature $`T_D^{}`$ instead of $`T_D`$ as will be discussed below. The coefficient $`\sigma _0`$ in Eq. (Slow oscillations of magnetoresistance in quasi-two-dimensional metals) can be estimated as $`\sigma _0=e^2N_{LL}2t^2d^2/\pi \mathrm{}^2\omega _ck_BT_D`$ where $`N_{LL}`$ is the Landau level degeneracy, $`N_{LL}/\mathrm{}\omega _c=m^{}/2\pi \mathrm{}^2`$ is the DoS at the Fermi level.

The second term in the curly brackets of Eq. (Slow oscillations of magnetoresistance in quasi-two-dimensional metals) describes the fundamental SdH oscillations modulated with the frequency $`F_{\text{beat}}=2t_{}m/e\mathrm{}`$. It is the last term in Eq. (Slow oscillations of magnetoresistance in quasi-two-dimensional metals) that gives the slow oscillations of $`\sigma _{zz}`$ (hence of $`R_{}1/\sigma _{zz}`$) with the frequency equal to the double beat frequency, in agreement with the experiment. The oscillations of $`F_{\text{beat}}`$ as a function of the angle $`\theta `$, in accordance with the AMRO effect Wosn3 ; Yam  lead to identical oscillations of $`F_{\text{slow}}`$ that explains the angular dependence plotted in Fig. 2. It is interesting that, unlike many anomalies which were studied earlier and associated with a strong enhancement of the oscillation amplitude and harmonic contents, the present phenomenon exists even in the limit of a weak amplitude and constant chemical potential, i.e. when no substantial deviations from the standard LK description were expected.

From Eq. (Slow oscillations of magnetoresistance in quasi-two-dimensional metals) it becomes clear why the slow oscillations are virtually independent of temperature (Fig. 1). Indeed, they are determined by the interlayer transfer integral (i.e. warping of the FS) but not by the electron energy. Therefore, they are not affected by the temperature smearing of the Fermi distribution function and do not contain the factor $`R_T`$. This is why the slow oscillations, though originating from the cyclotron motion on the same large orbits  on the FS cylinder as the fundamental SdH oscillations, persist up to relatively high temperatures at which the latter are completely suppressed.

Another notable feature of the slow oscillations is that their Dingle temperature $`T_D^{}`$ is different from $`T_D`$ entering the Dingle factor of the usual SdH oscillations. The usual Dingle temperature includes all mechanisms of the smearing of DoS oscillations. These are not only microscopical scattering events but also macroscopic spatial inhomogeneities of the sample. These inhomogeneities lead to macroscopic spatial variations of the electron energy in Eq. (1). Their effect is equivalent to the local shift of the chemical potential $`\mu `$. The total signal is an average over the entire sample and such macroscopic inhomogeneities lead to the damping of the magnetic quantum oscillations similar to that caused by temperature. Since the slow oscillations do not depend on $`\mu `$, they are not affected by this type of smearing and the corresponding Dingle temperature $`T_D^{}`$ is determined by only short-range scatterers. One can therefore estimate relative contributions from macroscopic inhomogeneities and from local defects to the scattering rate by comparing $`T_D`$ and $`T_D^{}`$. In Fig. 3 we plot the normalized amplitudes of the fundamental and slow oscillations as functions of inverse magnetic field. For the fundamental oscillations an angle $`\theta =19.3^{}`$ has been chosen. This angle corresponds to the AMRO peak, thus complications in the determination of the amplitude due to the beats could be avoided. Fits of the data according to Eq. (Slow oscillations of magnetoresistance in quasi-two-dimensional metals) (solid lines in Fig. 3) yield the Dingle temperatures $`T_D=(0.8\pm 0.02)`$ K and $`T_D^{}=(0.15\pm 0.02)`$ K. It is unlikely that the big difference between these values is only caused by the slightly different orientations. Therefore, we conclude that inhomogeneities play an important role in damping the SdH oscillations in the present sample.

Although the scattering is dominated by crystal imperfections at low temperature, electron-electron and electron-phonon interactions must also be taken into account in both $`R_D`$ and $`R_D^{}`$. This leads to a finite temperature dependence of the amplitude of the slow oscillations. A correct evaluation of this effect based on exact calculations of the imaginary part of the electron self-energy would allow to obtain additional information about the many-body interactions from the slow oscillations.

Finally, we note that the interference between the oscillating relaxation time and interlayer Fermi velocity considered above is not the only possible source of the slow oscillations. There are other oscillating quantities which do not enter Eq. (1) directly but are contained in the relaxation time condmat  and can, in principle, also lead to a similar effect. Taking into account these additional factors will lead to a change in the magnitude and phase of the slow oscillations, leaving, however, their dependence on temperature and magnetic field essentially the same.

In conclusion, we have shown that the slow oscillations in the Q2D metal $`\beta `$-(BEDT-TTF)2IBr2 originate from the warping of its cylindrical FS. We propose a model explaining this phenomenon as a general feature of clean Q2D metals developing when the cyclotron energy becomes comparable with the interlayer transfer energy. In particular, slow oscillations observed in some other layered organic metals Kar3 ; Wosn2 ; Ohmi ; togo  may have the same origin. Since the slow oscillations are not affected by the temperature smearing of the Fermi distribution function, one can use them for estimating the warping of the FS even at relatively high temperatures at which the fundamental SdH oscillations are suppressed. On the other hand, the temperature and field dependencies of the slow oscillations can give a valuable information on scattering processes.

The authors thank A.M. Dyugaev and I. Vagner for useful discussions. The work was supported by the EU ICN contract HPRI-CT-1999-40013, and grants DFG-RFBR No. 436 RUS 113/592 and RFBR No. 00-02-17729a.

Figure captions

Fig. 1. Interlayer resistance (a) and oscillating part of magnetic torque (b) of $`\beta `$-(BEDT-TTF)2IBr2 versus magnetic field at $`\theta 15^{}`$. The curves at different temperatures are offset for clarity.

Fig. 2. Angular dependencies of the frequency of the slow oscillations (a) and background resistance at $`B=14`$ T, $`T=0.44`$ K (b). Lines are guides for the eye.

Fig. 3. Amplitudes of the fast (a) and slow (b) oscillations normalized to the background resistance plotted versus inverse magnetic field at $`T=0.44`$ K. The lines are fits to Eq. (Slow oscillations of magnetoresistance in quasi-two-dimensional metals).


text:
# 1 Variations on Ω_{ϕ𝑖} lead to different physical situations given by 𝑤_{ϕ𝑜} and Ω_{ϕ𝑜}. We have taken Ω_{ϕ𝑖}={0.5,0.3,.2} short-dashed, long-dashed and solid lines, respectively. The vertical line marks the time at Ω_{ϕ𝑜}=0.7 with ℎ₀=0.67. Notice that 𝑤_{ϕ𝑜} increases with decreasing Ω_{ϕ𝑖}.

astro-ph/0203094

Acceptable Inverse Power Law Potential Quintessence with $`n=18/7`$

A. de la Macorraa111e-mail: macorra@fisica.unam.mx and G. Germanb222e-mail: gabriel@fis.unam.mx

| a Instituto de Física, UNAM |
| --- |
| Apdo. Postal 20-364, 01000 México D.F., México |
| b Centro de Ciencias Físicas, Universidad Nacional Autónoma de México, |
| Apartado Postal 48-3, 62251 Cuernavaca, Morelos, México |

ABSTRACT

We present a particle physics quintessence model which agrees well with existing cosmological data, including the position of the acoustic peaks. This model has an inverse power law potential (IPL) with $`n=18/72.57`$ and it gives $`w_{eff}=0.75`$, an acoustic scale $`l_A=307`$ and a density contrast $`\sigma _8=0.95`$.

Models with $`n>1`$ have been said to be disfavored by the analysis of the acoustic peaks. However, the results are not correct. The main reason is that the tracker approximation has been used in deriving the IPL constrains and for $`n<5`$ the scalar field has not reached its tracker value by present day.

The model can be derived from particle physics, using Affleck-Dine-Seiberg ”ADS” superpotential, for a non-abelian gauge group with $`N_c=8,N_f=1`$. The advantage of having $`N_f=1`$ is that there is only one degree of freedom below the condensation scale given by the condensate (quintessence) $`\varphi ^2=<Q\stackrel{~}{Q}>`$ field. The condensation scale is at $`1GeV`$ a very interesting scale since it connects the quintessence ”Q” with the standard model ”SM” scale. The similarity in energy scales between $`Q`$ and $`SM`$ scale gives an ”explanation” to the coincidence problem. The fact that only recently the universe is accelerating is a natural consequence of the Q scale and the evolution of $`\varphi `$.

The Maxima and Boomerang  observations on the cosmic radiation microwave background (”CMBR”) and the superonovae project SN1a  have lead to conclude that the universe is flat and it is expanding with an accelerating velocity. These conclusions show that the universe is now dominated by an energy density $`\mathrm{\Omega }_{\varphi o}=0.7\pm 0.1`$ (the subscript ”o” refers to present day quantities) with negative pressure. The SN1a data requires an equation of state $`w_{\varphi o}<2/3`$  while recent analysis on the CMBR peaks constrains the models to have $`w_{eff}=0.82_{.11}^{+.14}`$ , where $`w_{eff}`$ is an average equation of state. This energy is generically called the cosmological constant. Structure formation also favors a non-vanishing cosmological constant consistent with SN1a and CMBR observations . An interesting parameterization of this energy density is in terms of a scalar field with gravitationally interaction only called quintessence . The evolution of scalar field has been widely studied and some general approach con be found in . The evolution of the scalar field $`\varphi `$ depends on the functional form of its potential $`V(\varphi )`$ and a late time accelerating universe constrains the form of the potential .

One of the simplest and most interesting quintessence potentials are the inverse power law (IPL) . In some special cases they can be derived from non-abelian gauge theories  and we can also have consistent models with a gauge coupling unified with the standard model (SM) couplings .

From the CMBR analysis it has been inferred that IPL with $`n<1`$ are disfavored ,. However, in most cases one assumes a constant $`w`$ given by the tracker value $`w_{tr}=2/(2+n)`$ and for IPL models with $`n<5`$ the tracker solution  is not a good approximation to the numerical (exact) solution since the field has not reached its tracker value by present time. This fact implies that we cannot use the tracker solution (i.e. the constant $`w_{tr}`$) for the evolution of $`\varphi `$ in determining the acoustic peaks. It is no surprise that the values of the acoustic peaks differ greatly if we use the tracker solution approximation or we evolve the quintessence field from its initial conditions. Models with $`1<n<2.5`$ (for $`\mathrm{\Omega }_{\varphi i}0.25`$) are therefore still phenomenologically viable.

Tracker fields have the advantage that the value of $`w_{\varphi o}`$ does not depend on the initial condition. In fact one can have more then 100 orders of magnitude on the initial conditions $`\mathrm{\Omega }_{\varphi i}`$  and the value of $`w_{\varphi o}`$ will not change. In our case this is no longer so since the scalar field has not reached its tracker value and $`w_{\varphi o}`$ depends slightly on $`\mathrm{\Omega }_{\varphi i}`$, so there is a dependence on $`\mathrm{\Omega }_{\varphi i}`$ but there is no fine tuning required on the initial conditions (we do not need to adjust the initial condition more than one significant figure). Furthermore, in our model, derived from non-abelian gauge theories, we can determine the initial conditions in terms of the number of degrees of freedom of the system. So we do not need 100 orders of magnitude independence of initial conditions since they are well motivated and given within one order of magnitude at most. Of course, any model which requires a fine tuning on the initial conditions would not be theoretically acceptable but this is not our case. In both cases, tracker and our model, one still has the coincidence problem since the scale $`\mathrm{\Lambda }_c`$ has to be tuned so that $`\mathrm{\Omega }_\varphi 0.7`$ with $`h_o0.7`$ today.

There are two constrains on the equation of state parameter $`w_\varphi `$, one coming from direct observations SN1a which sets un upper limit, $`w_{\varphi o}<2/3`$  and the other is indirect and comes from numerical analysis of the CMBR data and gives a smaller value $`w_{eff}=0.82_{.11}^{+.14}`$ . Notice that the CMBR data gives a more negative $`w_\varphi `$ than the SN1a one but it is an average equation of state (from last scattering to present day) while the SN1a result gives an $`w_\varphi `$ in recent times. In IPL with $`n<5`$, where the quintessence field has not reached its tracker value yet, one has always $`w_{\varphi o}`$ larger than $`w_{eff}`$ in good agreement with SN1a and CMBR data. For IPL models it was shown  that $`w_{\varphi o}`$ depends on $`n`$ and the initial condition $`\mathrm{\Omega }_{\varphi i}`$. If we want $`w_{\varphi o}<2/3`$ assuming an $`\mathrm{\Omega }_{\varphi i}0.25`$ IPL models require an $`n`$ to be less 2.74  assuming no contribution from radiation at present time. If we include radiation with $`\mathrm{\Omega }_{ro}=4.17\times 10^5h_o^2`$ then the value of $`n`$ will decrease slightly, e.g for $`\mathrm{\Omega }_{\varphi i}=0.25`$ we have $`w_{\varphi o}2/3`$ for $`n2.5`$. Larger values of $`\mathrm{\Omega }_{\varphi i}`$ allow larger values of $`n`$, however we would not expect to have $`\mathrm{\Omega }_{\varphi i}`$ much larger since a ”reasonable” amount of energy must go into the standard model of elementary particles ”SM”. These results set an upper value of $`n`$ but there is still room for models with $`1<n<2.5`$ and if we take $`\mathrm{\Omega }_\varphi 0.3`$ then the value of $`n`$ can be as large as $`n2.66`$.

The CMBR constrain on $`w_\varphi `$ can be studied from the position of the third acoustic peak. The position of the third CMBR peak has been found to be not very sensitive to the different cosmological parameters and it is a good quantity to obtain the acoustic scale $`l_A`$ . The acoustic scale $`l_A`$, which sets the scale of the peaks, derived from the third acoustic peak is

$$l_A=316\pm 8$$
(1)
were we have taken $`l_3=845_{25}^{+12}`$  (see below for the definition of $`l_A`$).

Here we will present a model with the largest value of $`n`$ that is still in agreement with the observational cosmological data  and that can be nicely derived from particle physics. Since, for a non-abelian gauge group (see below) we have $`n=2+4N_f/(N_cN_f)=2+4/(z1)`$, where $`N_f`$ is the number of flavors of the gauge group $`SU(N_c)`$ and $`z=N_c/N_f`$, we can see that $`n`$ decreases with increasing $`z`$. The requirement on $`n`$ to be smaller than $`2.66(2.5)`$ for $`\mathrm{\Omega }_{\varphi i}0.3(0.25)`$, in order to give the correct phenomenology and the observed values of the acoustic scale and present day $`w_{\varphi o}`$, implies that $`z>9(7)`$ or $`zN_f<N_c`$. As we will see later the number of condensates of the gauge group $`SU(N_c)`$ is given in terms of $`N_f`$. Therefore, the model with the least number of condensates has $`N_f=1`$ and the smallest gauge group would be $`N_c=8`$. Furthermore, in string compactification for the heterotic string, the gauge group has at most rank 8 (as $`SU(8)`$). The value of $`n`$ for $`SU(8)`$ with $`N_f=1`$ is $`n=18/7=2.57`$. This model represents the limiting acceptable model.

The acoustic scale gives $`l_A=307`$ and $`w_{\varphi o}=0.68`$, which gives a good prediction of the acoustic peaks within observational limits. This acoustic scale should be compared with the tracker solution $`l_{Atr}=281`$ and $`w_{tr}=2/(2+n)=0.44`$ and the cosmological constant $`l_{ACte}=315`$ (i.e. $`w_{Cte}1`$) for the same initial conditions. We see that the tracking solution is not a good approximation since $`w_{\varphi o}`$ differs by more than $`38\%`$ and the acoustic scale $`l_A`$ by $`9\%`$ from the numerical solution of the scalar field, discrepancy large enough to rule out the model.

Since our model has $`n<5`$ the quintessence field has not reached its tracker value yet and the coincidence problem is not solved. However, there is a clear connection between the model condensation scale $`\mathrm{\Lambda }_c=1GeV`$ and the standard model scale. So, we could think of the ”solution” to the coincidence problem as the following: The scale of quintessence should not be given by today’s energy density but by the condensation scale $`\mathrm{\Lambda }_c`$. The natural value of this scale is that of the standard model. The subsequent evolution of the quintessence field is determined by the solution of the Friedmann-Robertson-Walker equations and the fact that only recently the universe is accelerating is a natural consequence of the quintessence dynamics starting at $`\mathrm{\Lambda }_c`$.

The model with $`n=18/7`$ can be easily obtained from a $`SU(8)`$ non-abelian gauge group with $`N_f=1`$ number of (chiral + antichiral) fields in the fundamental representation and with a condensation scale $`\mathrm{\Lambda }_c=1GeV`$, quite an interesting scale. Above the condensation scale the gauge coupling constant is small and the elementary fields are massless. At the condensation scale there is a phase transition, the gauge coupling constant becomes strong, and binds the elementary fields together forming meson fields. In this model there is only one degree of freedom, $`\varphi `$, in the confined phase and the Affeck-Dine-Seiberg ”ADS” superpotential  obtained is therefore exact.

At the beginning we have particles of the standard model (SM) and the quintessence model (Q). All fields, SM and $`Q`$ model, are massless and redshift as radiation until we reach the condensation scale $`\mathrm{\Lambda }_c`$ of Q group. Below this scale the fields of the quintessence gauge group will dynamically condense and we use ADS potential to study its cosmological evolution. The ADS potential is non-perturbative and exact (it receives no quantum corrections)  and it is given for a non-abelian $`SU(N_c)`$ gauge group with $`N_f`$ (chiral + antichiral) massless matter fields by

$$W=\left(N_cN_f\right)\left(\frac{\mathrm{\Lambda }_c^{b_o}}{det<Q\stackrel{~}{Q}>}\right)^{1/\left(N_cN_f\right)}$$
(2)
where $`b_o=3N_cN_f`$ is the one-loop beta function coefficient. The scalar potential in global supersymmetry is $`V=|W_\varphi |^2`$, with $`W_\varphi =W/\varphi `$, giving

$$V=c^2\mathrm{\Lambda }_c^{4+n}\varphi ^n$$
(3)
where we have taken $`det<Q\stackrel{~}{Q}>=\mathrm{\Pi }_{j=1}^{N_f}\varphi _j^2,c=2N_f,n=2+4\frac{N_f}{N_cN_f}`$ and $`\mathrm{\Lambda }_c`$ is the condensation scale of the gauge group $`SU(N_c)`$. Our model has $`N_f=1,N_c=8`$ and $`n=18/7`$. There are no baryons since $`N_f<N_c`$ and there is only one degree of freedom below $`\mathrm{\Lambda }_c`$ which is the condensate $`\varphi =<Q\stackrel{~}{Q}>`$.

The cosmological evolution of $`\varphi `$ with an arbitrary potential $`V(\varphi )`$ can be determined from a system of differential equations describing a spatially flat Friedmann–Robertson–Walker universe in the presence of a barotropic fluid energy density $`\rho _\gamma `$ that can be either radiation or matter, are

$`\dot{H}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\rho _\gamma +p_\gamma +\dot{\varphi }^2\right),`$
$`\dot{\rho }`$ $`=`$ $`3H\left(\rho +p\right),`$ (4)
$`\ddot{\varphi }`$ $`=`$ $`3H\dot{\varphi }{\displaystyle \frac{dV\left(\varphi \right)}{d\varphi }},`$
where $`H`$ is the Hubble parameter ($`H=100hkmMpc^1s^1`$), $`\dot{f}=df/dt`$, $`\rho `$ ($`p`$) is the total energy density (pressure) and we are setting the reduced Planck mass $`m_p^2=1/8\pi G1`$.

Solving eqs.(S0.Ex1) we have that the energy density of the $`Q`$ group $`\mathrm{\Omega }_\varphi `$ drops quickly, independently of its initial conditions, and it is close to zero for a long period of time, which includes nucleosynthesis (NS) if $`\mathrm{\Lambda }_c`$ is larger than the NS energy $`\mathrm{\Lambda }_{NS}`$ (or temperature $`T_{NS}=0.110MeV`$), and becomes relevant only until very recently . On the other hand, if $`\mathrm{\Lambda }_c<\mathrm{\Lambda }_{NS}`$ then the NS bounds on relativistic degrees of freedom must be imposed on the models. Finally, the energy density of $`Q`$ grows and it dominates at present time the total energy density with the $`\mathrm{\Omega }_{\varphi o}0.7`$ and a negative pressure $`w_{\varphi o}<2/3`$ leading to an accelerating universe .

The value of the condensation scale in terms of $`H_o`$ is

$$\mathrm{\Lambda }_c=\left(\frac{3y_o^2\varphi ^nH_o^2}{4N_f^2}\right)^{\frac{1}{4+n}}$$
(5)
The approximated value can be obtained since one expects, in general, to have $`y_o^2\varphi _o^n1`$ for a model with $`\mathrm{\Omega }_{\varphi o}=0.7`$ and $`w_{\varphi o}<2/3`$. The magnitude order of the condensation scale is therefore $`\mathrm{\Lambda }_c=H_o^{2/(4+n)}`$.

In our model, the cosmological evolution requires a condensation scale $`\mathrm{\Lambda }_c=1GeV`$ in order to give $`\mathrm{\Omega }_{\varphi o}=0.7\pm 0.1`$ with a Hubble parameter $`h_o=.65\pm .7`$ at present time. Since $`\mathrm{\Lambda }_c\mathrm{\Lambda }_{NS}`$ the energy density at NS is $`\mathrm{\Omega }_\varphi (NS)1`$ and there is no constrain from nucleosynthesis on the model.

In order to set the initial conditions for $`\varphi `$ we will assume that all relativistic degree of freedom (MSSM and Q) had the same fraction of energy density at high energies, when all fields were massless. The initial conditions could be set at the unification scale or at the reheating temperature $`T_{rh}`$. If $`T_{rh}`$ is larger then the supersymmetric masses (i.e. $`T_{rh}>10^3GeV`$), which is a natural assumption, then all degrees of freedom (MSSM and Q) would be relativistic at $`T_{rh}`$ and each degree of freedom would have the same energy density (assuming the standard reheating process which is gauge blind). Therefore, the initial energy density conditions would be exactly the same at the unification scale or reheating temperature.

The MSSM has $`g_{smi}=228.75`$ while the $`Q`$ group has $`g_{Qi}=(1+7/8)(2(N_c^21)+2N_fN_c)=266.25`$ degrees of freedom, where $`g_a=\mathrm{\Sigma }_aBosons+7/8\mathrm{\Sigma }_aFermions`$. Taking into account that some fields become massive at lower energies, we can determine the energy density at an arbitrary energy scale $`\mathrm{\Lambda }`$ and it is given by

$$\mathrm{\Omega }_Q\left(\mathrm{\Lambda }\right)=\frac{g_{Qf}\left(g_{smf}g_{Qi}/g_{smi}g_{Qf}\right)^{4/3}}{g_{smf}+g_{Qf}\left(g_{smf}g_{Qi}/g_{smi}g_{Qf}\right)^{4/3}}$$
(6)
where $`g_{smi},g_{smf},g_{Qi},g_{Qf}`$ are the initial (i.e. at high energy scale) and final (i.e. at $`\mathrm{\Lambda }`$) standard model and $`Q`$ model relativistic degrees of freedom, respectively. Taking $`\mathrm{\Lambda }=\mathrm{\Lambda }_c=1GeV`$ the MSSM has $`g_{smf}=10.75`$ and if the $`Q`$ group is still supersymmetric at $`\mathrm{\Lambda }_c`$ it has $`g_{Qf}=g_{Qi}`$ and $`\mathrm{\Omega }_{\varphi i}(\mathrm{\Lambda }_c)=0.29`$. If $`Q`$ is no longer supersymmetric $`g_{Qf}=g_{Qi}/2`$ and $`\mathrm{\Omega }_{\varphi i}(\mathrm{\Lambda }_c)=0.34`$. So a reasonable choice for the initial conditions is $`\mathrm{\Omega }_{\varphi i}(\mathrm{\Lambda }_c)=0.3`$.

We show in fig.(1) the evolution of $`\mathrm{\Omega }_\varphi `$ and $`w_\varphi `$ as a function of $`N=Log(a)`$, with $`a`$ the scale factor, for initial conditions $`\mathrm{\Omega }_{\varphi i}=0.5,0.3,0.2`$ short-dashed, long-dashed and solid lines, respectively. We see that $`w_{\varphi o}`$ decreases for larger initial condition $`\mathrm{\Omega }_{\varphi i}`$. For $`\mathrm{\Omega }_{\varphi i}=0.9,0.5,0.3,.2`$ one finds $`w_{\varphi o}=0.95,0.8,0.68,0.61`$ and an acoustic scale $`l_A=314,313,307,303`$, respectively. It is no surprise that for $`\mathrm{\Omega }_{\varphi i}=0.2`$ the value of $`w_{\varphi o}=0.61`$ lies outside the observed range $`w_{\varphi o}<2/3`$. This is because the model we are working with (i.e. n=18/7=2.57) gives almost the limiting value of $`w_{\varphi o}=2/3`$ for $`\mathrm{\Omega }_{\varphi i}=0.3`$ ( $`w_{\varphi o}`$ increases for smaller $`\mathrm{\Omega }_{\varphi i}`$) and that was the reason for using this model. However, as long as we take $`\mathrm{\Omega }_{\varphi i}0.27`$ the model satisfies the cosmological constrains. For scalar fields that have reached its tracker value by present day (i.e. $`n>5`$) the initial value of $`\mathrm{\Omega }_{\varphi i}`$ is not constrained  since it can vary for more than 100 orders of magnitude and the value of $`w_{\varphi o}`$ will be the same $`w_{\varphi o}=2/(2+n)`$ , however for $`n>5`$ one has $`w_{\varphi o}0.28`$ which is too large.

In previous works  we have studied quintessence models that are unified with the standard model gauge groups, i.e. the gauge coupling constant of all gauge groups is the same at the unification scale. In the model we are working here, $`N_f=1,N_c=8`$, the renormalization group equation given by $`\mathrm{\Lambda }_{RG}=\mathrm{\Lambda }_{Gut}exp[16\pi ^2/2b_og_{gut}^2]=9\times 10^{12}GeV`$ with $`\mathrm{\Lambda }_{gut}=10^{16}GeV,g_{gut}^2=4\pi /25.7`$ the unification scale and coupling , respectively, and $`b_o=3N_cN_f=23`$ the one-loop beta function coefficient. It is clear that $`\mathrm{\Lambda }_{RG}\mathrm{\Lambda }_c`$ and so the model cannot be unified with the SM groups. If we insist in gauge coupling unification we would need, on top of the original 1 chiral + 1 antichiral fields, 37 extra chiral fields. If all extra fields are chiral then they would not contribute to the ADS superpotential in eq.(2). Of course, we think that such a model is not natural but it is possible (4-D string models can have different number of chiral and antichiral fields).

The acoustic scale, that sets the scale of acoustic peaks, for a flat universe is given by

$$l_A=\pi \frac{\tau _o\tau _{ls}}{\overline{c}_s\tau _{ls}}$$
(7)
where $`\tau _o`$ and $`\tau _{ls}`$ are the conformal time today and at last scattering ($`\tau =𝑑ta^1(t),a(t)`$ the scale factor) and $`\overline{c}_s\tau _{ls}^1_0^{\tau _{ls}}𝑑\tau c_s`$ is the average sound speed before last scattering ($`c_s^2=3+(9/4)w_b(t)/w_r(t)`$ with $`w_b=\mathrm{\Omega }_bh^2,w_r=\mathrm{\Omega }_rh^2`$ the fraction of baryon and radiation energy density, respectively). The acoustic m-th peak $`l_m`$ is then given in terms of $`l_A`$ and a peak and model dependent phase shift $`\phi _m`$, $`l_m=l_A(m\phi _m)`$. It has been observed in  that the third peak is quite insensitive to different cosmological parameters that enter in determining $`\phi _30.341`$ and so the position of the third peak is a good quantity to extract the acoustic scale $`l_A`$. The data from Boomerang and Maxima set the first three acoustic peaks, (the first through $`l_{3/2}`$) and the acoustic scale at

$$l_1=213_{13}^{+10},l_2=541_{32}^{+20},l_3=845_{25}^{+12},l_{3/2}=416_{32}^{+22},l_A=316_8^{+8}.$$
(8)
We have solved eqs.(S0.Ex1) numerically with initial conditions $`\mathrm{\Omega }_{\varphi i}(\mathrm{\Lambda }_c)=0.3`$ at $`\mathrm{\Lambda }_c=1GeV`$ imposing $`h_o=0.65`$, $`\mathrm{\Omega }_{\varphi o}=0.75`$ and we obtain

$$l_A=307,w_{\varphi o}=0.68,w_{eff}=0.75,$$
(9)
where $`w_{eff}𝑑a\mathrm{\Omega }(a)w_\varphi (a)/𝑑a\mathrm{\Omega }(a)`$. The energy density al last scattering (LS) is negligible ($`\mathrm{\Omega }_\varphi (LS)=10^9`$) and the average sound speed at LS is $`\overline{c}_s=0.52`$. The result is not highly sensitive to the initial conditions and a change in $`\mathrm{\Omega }_{\varphi i}`$ of $`50\%`$ will still be ok .

We see from eq.(9) that the acoustic scale $`l_A`$ is within the observational range given by eq.(1). The prediction of the first three acoustic peaks and first through ($`l_{3/2}`$) is

$$l_1=223,l_2=542,l_3=829,l_{3/2}=414$$
(10)
for a baryon density $`w_b=\mathrm{\Omega }_bh_o^2=0.02`$ and $`n_s=1`$ the index of power spectrum of primordial density fluctuations. The peak values in eq.(10) are consistent with the observational data in eq.(8) and we have used the phase shifts given in .

The value of $`l_A`$ is sensitive to $`\mathrm{\Omega }_{\varphi o}`$ and even more to $`h_o`$. For increasing $`h_o`$ (with $`\mathrm{\Omega }_{\varphi o}`$ fixed) we find a decreasing $`l_A`$, e.g. for $`h_o=0.6,\mathrm{\hspace{0.17em}0.65},\mathrm{\hspace{0.17em}0.7},\mathrm{\hspace{0.17em}0.8}`$ one has $`l_A=309,\mathrm{\hspace{0.17em}307},\mathrm{\hspace{0.17em}288},\mathrm{\hspace{0.17em}272}`$ respectively, while for an increasing energy density one has an increasing $`l_A`$, e.g. $`\mathrm{\Omega }_{\varphi o}=0.6,\mathrm{\hspace{0.17em}0.7},\mathrm{\hspace{0.17em}0.75},\mathrm{\hspace{0.17em}0.8}`$ one gets $`l_A=291,\mathrm{\hspace{0.17em}300},\mathrm{\hspace{0.17em}307},\mathrm{\hspace{0.17em}317}`$ with fixed $`h_o=0.65`$. Since $`\overline{c}_s`$ depends on $`w_b`$ and an increase in $`w_b`$ makes $`c_s`$ smaller and we have therefore a slight increase in the acoustic scale, e.g. for $`w_b=.019,\mathrm{\hspace{0.17em}0.020},\mathrm{\hspace{0.17em}0.022},\mathrm{\hspace{0.17em}0.026}`$ one gets a value of $`l_A=305,\mathrm{\hspace{0.17em}307},\mathrm{\hspace{0.17em}309},\mathrm{\hspace{0.17em}314}`$ with $`\mathrm{\Omega }_{\varphi o}=0.75,h_o=0.65`$ fixed. A change in $`n_s`$ does not affect $`l_A`$ but it changes the acoustic peaks through the phase shifts $`\phi _m`$ slightly.

Another relevant cosmological quantity is the density contrast on scales of $`8h^1Mpc`$, $`\sigma _8`$, which is constraint by the galaxies cluster abundance. For a flat universe the empirical fit of different authors (which converge within one $`\sigma `$) are: Eke et al. have $`\sigma _8=(0.52\pm 0.08)\mathrm{\Omega }_m^{0.52+0.13\mathrm{\Omega }_m}`$, Viana et al. report their best fit at $`\sigma _8=0.56\mathrm{\Omega }_m^{0.47}(1\pm 0.3)`$ , while Steinhardt et al.  have $`\sigma _8=[(0.50.1\mathrm{\Theta })\pm 0.1]\mathrm{\Omega }_m^\gamma ,\mathrm{\Theta }=(n1)+(h_o0.65),\gamma =0.210.22w+0.25\mathrm{\Theta }`$. For $`n=1,h_0=0.65`$ one has $`\sigma _8=1.02\pm .15,\mathrm{\hspace{0.17em}1.07}\pm 0.03,\mathrm{\hspace{0.17em}0.98}\pm 0.2`$ for Eke, Viana and Steinhardt respectively. Recent analysis give slightly lower values of $`\sigma _8=(0.46_{0.07}^{+0.05})\mathrm{\Omega }_m^{0.52}`$ for  (see also ) and depends quite strongly on $`\mathrm{\Omega }_m`$ (for smaller $`\mathrm{\Omega }_m`$ one has a larger $`\sigma _8`$). The central value for $`\mathrm{\Omega }_m=0.25`$ is $`\sigma _8=0.945`$ which agrees quite well with our the value obtained in our model $`\sigma _8=0.95`$.

To conclude, we have shown that inverse power law potentials with $`n>1`$ are not disfavored with the existing cosmological data and we have shown an explicit example with $`n2.57`$. In particular, the values of $`w_{\varphi o}`$, the acoustic scale and peaks and the density contrast $`\sigma _8`$ lie within the observed data. The model has been derived from particle physics, using ADS superpotential, from a non-abelian gauge group with $`N_c=8,N_f=1`$. Since $`N_f=1`$ there is only one degree of freedom below the condensation scale, i.e. the condensate or quintessence field $`\varphi `$. The condensation scale is at $`1GeV`$ a very interesting scale which connects quintessence with the standard model.

This work was supported in part by CONACYT project 32415-E and DGAPA, UNAM project IN-110200.