format: text

text:
# Shocks in Vertically Oscillated Granular Layers

## I Introduction

### I.1 Background

A successful continuum description of granular media could exploit the powerful methods of standard fluid dynamics to describe a variety of granular flow phenomena, but the validity of a continuum description of granular materials has been questioned and remains an open problem Kadanoff (1999). Several investigators have proposed equations of motion for rapid granular flow using continuum fields – number density $`n`$, velocity $`𝐮`$, and granular temperature $`T`$ ($`\frac{3}{2}T`$ is the average kinetic energy due to random particle motion) Haff (1983); Campbell (1990); Dufty (2001). In one approach, the interaction between grains is modeled with binary, inelastic hard sphere collision operators in kinetic theory to derive granular continuum equations to Euler Goldshtein and Shapiro (1995), Navier-Stokes Jenkins and Richman (1985), and Burnett Sela and Goldhirsch (1998) order. However, attempts to use such continuum equations to predict properties of granular media have mainly been restricted to idealized systems that are difficult to produce in experiment Haff (1983); Goldshtein et al. (1996); Jenkins (1998); Kamenetsky et al. (2000) and to steady state or asymptotic time limits Haff (1983); Jenkins (1998); Goldshtein et al. (1996); Rericha et al. (2002). The viability of these continuum equations for experimentally realizable, three-dimensional (3D), time-dependent granular systems has scarcely been examined. In this paper, we use 3D inelastic hard sphere molecular dynamics (MD) simulations as well as 3D simulations of continuum equations to Navier-Stokes order to investigate the strongly time-dependent properties of normal shocks in a vertically oscillated granular layer.

### I.2 Shock waves

If the Mach number $`Ma`$ (the ratio of the local mean fluid speed to the local speed of sound) is greater than unity at the point where a fluid encounters an obstacle, a compression wavefront is formed near the object and steepens to form a shock. In an ideal fluid with no viscosity or heat conduction, the wavefront steepens until the fields develop a mathematical discontinuity between the compressed region and the undisturbed region. In a viscous, conducting fluid, the fields vary continuously, but the front steepens until the shock thickness (the width of the region in which the fields change from the undisturbed values to the compressed values) is on the order of a mean free path.

When a supersonic fluid impinges perpendicularly onto a flat, impenetrable plate, a shock forms normal to the plate. The normal component of the mean velocity (with respect to the plate) is zero at the plate, but the velocity of the undisturbed fluid far from the plate is still supersonic, flowing towards the plate. The Mach number decreases abruptly from greater than unity to near zero as the plate is approached. The density, temperature, and pressure also change rapidly in the region of the shock front. This normal shock propagates back upstream in a direction opposite that of the upstream flow velocity.

This paper concerns normal shocks in vertically oscillating granular layers. A distinguishing feature of granular materials is that granular flows reach supersonic speeds under common laboratory conditions Haff (1983); Rericha et al. (2002); Goldshtein et al. (1995). Hence shock formation, which is achieved only under extreme conditions in ordinary gases, is commonplace in granular materials. Shock propagation has been extensively studied in rarefied gases Courant and Friedrichs (1948); Zel’dovich and Raizer (1966), but the similarities and differences between shocks in ordinary fluids and in granular media is the subject of ongoing research Rericha et al. (2002); Potapov and Campbell (1996); Goldshtein et al. (1995). A recent laboratory and molecular dynamics study examined time-independent behavior of oblique shocks in granular flow between two closely spaced plates, and the results were compared with 2D simulations of Navier-Stokes order continuum equations Rericha et al. (2002).

### I.3 Model system

The system consists of a layer of grains on an impenetrable plate that oscillates sinusoidally in the direction of gravity with amplitude $`A`$ and frequency $`f`$. The dimensionless peak acceleration is $`\mathrm{\Gamma }=4\pi ^2f^2A/g`$, where $`g`$ is the acceleration due to gravity. This simple oscillating system has been found to yield standing wave pattern formation Melo et al. (1994), convection Knight et al. (1996), clustering Falcon et al. (1999), steady-state flow fields far from the plate Brey et al. (2001), and shocks Goldshtein et al. (1995). We examine shock formation and propagation in layers that are approximately 9 particles deep as poured except in Sec. III E, where layer depth is varied. Throughout this paper, the phrase “the layer” refers to the dense region in which the volume fraction of particles is greater than $`4\%`$ of the random close-packed volume fraction, $`\nu _{max}=0.65`$.

When the maximum acceleration of the plate is greater than $`g`$, the layer leaves the plate during each cycle. When the layer is off the plate, it is cooled by inelastic collisions between particles, while the particles are simultaneously accelerated towards the plate by gravity; this leads to a large mean velocity compared to the sound speed (See Sec. III C). We show in Sec. III that when the layer contacts the plate later in the cycle, the Mach number is much greater than unity with respect to the plate, i.e., the flow is supersonic.

Shocks similar to those presented here were found for a wide range of $`\mathrm{\Gamma }`$ and $`f`$. Here we report results obtained from a particular set of parameter values: $`\mathrm{\Gamma }=3`$ and $`f=0.095\sqrt{g/\sigma },`$ where $`\sigma `$ is the diameter of a particle. This would correspond to a frequency of approximately $`30`$ Hz for particles with a diameter of $`0.1`$ mm. The particles have a coefficient of restitution $`e=0.9,`$ except in Sec. III D, where $`e`$ is varied to investigate the effect of inelasticity on shock propagation. The number of particles in the container per unit area of the bottom plate is $`10/\sigma ^2`$, except in Sec. III E, where the number of particles is varied to investigate shocks in layers of different depths. That is, the layer would be $`10\sigma `$ deep if the particles were arranged in a simple cubic lattice. In actual packings seen experimentally, $`10`$ particles/$`\sigma ^2`$ corresponds to an average depth of approximately $`H=9\sigma `$ as poured Bizon et al. (1998). For consistency with previous investigations, we use $`H`$ (the depth of the layer as poured) to express the depth of the layer throughout this paper.

The values of $`\mathrm{\Gamma }`$ and $`f`$ in this study were examined previously in experiments and MD studies of spatial pattern formation in shallow granular layers Bizon et al. (1998). In the present simulations, pattern formation was intentionally suppressed by considering a container smaller than one full wavelength in either horizontal direction (with periodic boundary conditions at the side walls); the absence of patterns permits better visualization of the dependence of the fields on height $`z`$.

We present 3D continuum and MD simulations of the oscillated granular system, and analyze the behavior of the fields $`n,`$ $`T,`$ and $`Ma`$. Section II describes the methods used to simulate and analyze the oscillatory time-dependent flow. Section III reports our results, including an examination of varying inelasticity and layer depth, and Sec. IV presents our conclusions.

## II Methods

### II.1 Molecular dynamics simulation

We use an inelastic hard sphere molecular dynamics simulation that was validated in previous studies Bizon et al. (1998); Moon et al. (2002) by comparison to experimentally observed standing wave patterns for varying control parameters. The collision model assumes instantaneous binary collisions. Linear momentum is conserved, but energy is dissipated through inelastic collisions. As in the continuum model (Sec. II B), the dissipation is assumed to be characterized by a single parameter, the normal coefficient of restitution $`e`$:

$$e=v_n^{}/v_n,$$
(1)
where $`v_n=(𝐯_1𝐯_2)(𝐫_1𝐫_2)/|𝐫_1𝐫_2|`$, and $`𝐯`$’s and $`𝐫`$’s are velocity and position vectors for a pair of colliding particles before the collision (no superscript) and after the collision (with a superscript $``$). Particles are monodisperse, and $`|𝐫_1𝐫_2|`$ and $`|𝐫_1^{}𝐫_2^{}|`$ are identically $`\sigma `$. In previous studies, the normal coefficient of restitution depended on $`v_n`$ Bizon et al. (1998). We investigate the effects of such a model in Sec. III D. However, to be consistent with the continuum model, we use a velocity-independent coefficient of restitution for most of this study. We also neglect the surface friction between particles as well as between the particles and the plate.

The particles are constrained between the bottom plate which oscillates sinusoidally between height $`z=0`$ and $`z=2A`$, and a ceiling which is fixed at $`z=200\sigma `$. In the two horizontal directions, the width is $`20\sigma `$, and the boundary condition is periodic. Except in Sec. III E, 3937 particles were used; the number of particles was chosen to correspond to $`H=9\sigma `$ and the total mass of the layer matches with that in the continuum simulations. In Sec. III E, 1969 particles were used to correspond to $`H=4.5\sigma `$, and 5906 particles were used to correspond to $`H=13.5\sigma `$. The container has the same $`e`$ as the particles, and the mass of the container is assumed to be infinitely large compared to that of the granular layer.

Fields are calculated by dividing the container height into small bins of size $`\sigma `$. The number density is the number of particles found in a bin divided by the volume of the bin. The number density $`n`$ from MD simulations is averaged over the same phase angle of the oscillatory state for ten cycles of the plate. The velocity $`𝐯`$ of each particle is used to calculate the mean velocity in each bin, $`𝐮=𝐯`$, where the brackets represent the average over all the particles found in a bin at the same phase angle of the oscillatory state for ten cycles of the plate. Granular temperature is then defined as $`T=\frac{1}{3}|𝐯𝐮|^2`$. Throughout this paper, we use units such that the particles have mass unity.

### II.2 Continuum simulation

We numerically solve the continuum equations of Navier-Stokes order proposed by Jenkins and Richman Jenkins and Richman (1985). These equations were derived for a dense gas composed of frictionless (smooth), inelastic hard spheres by applying Grad’s 13-moment method to the inelastic Enskog-Boltzmann equation. This model yields hydrodynamic equations for number density (or equivalently, volume fraction $`\nu =\frac{\pi }{6}n\sigma ^3`$), momentum, and granular temperature:

$$\frac{n}{t}+(n𝐮)=0,$$
(2)
$$n\left(\frac{𝐮}{t}+𝐮𝐮\right)=\underset{¯}{𝐏}ng\widehat{𝐳},$$
(3)
$$\frac{3}{2}n\left(\frac{T}{t}+𝐮T\right)=𝐪+\underset{¯}{𝐏}:\underset{¯}{𝐄}\gamma ,$$
(4)
where the components of the symmetrized velocity gradient tensor $`\underset{¯}{𝐄}`$ are given by: $`E_{ij}=\frac{1}{2}\left(_ju_i+_iu_j\right).`$ The components of the stress tensor $`\underset{¯}{𝐏}`$ are given by the constitutive relation:

$$P_{ij}=\left[p+(\lambda \frac{2}{3}\mu )E_{kk}\right]\delta _{ij}+2\mu E_{ij},$$
(5)
and the heat flux is calculated from Fourier’s law:

$$𝐪=\kappa T.$$
(6)
To calculate the pressure, we use the equation of state and radial distribution function at contact proposed by Goldshtein et al. Goldshtein et al. (1996) to include both dense gas and inelastic effects:

$$p=nT\left[1+2(1+e)G(\nu )\right],$$
(7)
where

$$G(\nu )=\nu g_0(\nu ),$$
(8)
and the radial distribution function at contact, $`g_0`$, is:

$$g_0(\nu )=\left[1\left(\frac{\nu }{\nu _{max}}\right)^{\frac{4}{3}\nu _{max}}\right]^1,$$
(9)
where $`\nu _{max}=0.65`$ is the 3D random close-packed volume fraction.

Equations (2-4) differ from those for a compressible, dense gas of elastic particles by the energy loss term

$$\gamma =\frac{12}{\sqrt{\pi }}(1e^2)\frac{nT^{3/2}}{\sigma }G(\nu ),$$
(10)
which arises from the inelasticity of collisions between particles. The bulk viscosity is given by

$$\lambda =\frac{8}{3\sqrt{\pi }}n\sigma T^{1/2}G(\nu ),$$
(11)
the shear viscosity by

$$\mu =\frac{\sqrt{\pi }}{6}n\sigma T^{1/2}\left[\frac{5}{16}\frac{1}{G(\nu )}+1+\frac{4}{5}\left(1+\frac{12}{\pi }\right)G(\nu )\right],$$
(12)
and the thermal conductivity by

$$\kappa =\frac{15\sqrt{\pi }}{16}n\sigma T^{1/2}\left[\frac{5}{24}\frac{1}{G(\nu )}+1+\frac{6}{5}\left(1+\frac{32}{9\pi }\right)G(\nu )\right].$$
(13)
In marked contrast to incompressible isothermal liquid layers, $`n`$ and $`T`$ for an oscillating granular layer vary by orders of magnitude in both space and time during each cycle of oscillation; hence the viscosities and thermal conductivity are very far from constant and must be calculated at each timestep.

The equations are numerically integrated using a second order finite difference scheme on a uniform grid in 3D and first order adaptive time stepping. For densities near close-packed, terms of order $`1/G(\nu )`$ are negligible in Eqs. (12) and (13). For reasons of numerical stability, these terms were neglected and artificial dissipation Roache (1976) was introduced in the low-density region above and below the layer. In the bulk of the layer, where the shock forms and propagates, the density near the shock is high, and the artificial dissipation has little effect on the solution.

As in the MD simulation, the granular fluid is contained between two impenetrable horizontal plates at the top and bottom of the container. We initially conducted the simulation in a cell with a ceiling that oscillates with the bottom plate at a height $`200\sigma `$ above the bottom plate. For $`H=9\sigma `$, where most of these simulations were conducted, an average of less than one particle per oscillation reached a height of $`z80\sigma `$ in the MD simulation. Changing the height of the cell from $`200\sigma `$ to $`80\sigma `$ in the continuum simulation resulted in a significant speedup in computation time, but did not change the shock profile. The fractional root mean square difference over one cycle between the shock location produced in the tall cell ($`200\sigma `$ high) and the shorter cell ($`80\sigma `$ high) was less than $`1\%`$. Therefore, for computational efficiency, the continuum simulations were conducted with a ceiling located at a height $`80\sigma `$ above the bottom plate, oscillating with the bottom plate. In MD simulations of shallow layers, $`H5`$, many particles reach heights of $`z>80\sigma `$, so we use a cell height of $`200\sigma `$ in continuum simulations for $`H5`$ (see Sec. III E).

As in the MD simulations, we use periodic boundary conditions in the horizontal directions. Initial simulations were conducted with a container $`20\sigma `$ wide in the horizontal directions, as in the MD simulation. For reasons of computational efficiency, we reduced the horizontal width of the cell to $`11\sigma `$, which produced no significant changes in the shock behavior. The width in the horizontal directions is $`11\sigma `$ for continuum simulations throughout this paper.

The continuum equations are solved in the reference frame of the container, with the container oscillations accounted for by a sinusoidal external forcing term in the momentum equation, Eq. (3). We use impenetrable boundary conditions at the plates: $`u_z=0`$ in the frame of the plate. Boundary conditions for horizontal velocities and temperature are set to match those found in the MD simulation. For most of the oscillation, their vertical derivatives at the plates were found to be negligible in the MD simulation, so for simplicity the boundary conditions used for the continuum simulation were: $`u_x/z=u_y/z=0,`$ and $`T/z=0,`$ although these derivatives are not identically zero in the MD simulation for the entire oscillation cycle of the plate.

### II.3 Shock tracking and shock width

With each cycle of the plate, a normal shock wave is formed and propagates through the layer. This shock wave separates a compressed region near the plate from the undisturbed region which is still falling towards the plate. Both the undisturbed region and the compressed region vary in space and time. In Sec. III we find that when the shock is in the layer, granular temperature increases with height in the compressed region, drops quickly in the shock, and then increases again in the undisturbed region (Fig. 1). The shock is identified in simulations as the region in which the temperature falls from its maximum at the top of the compressed region to the local minimum at the bottom of the undisturbed region. The shock width $`d`$ is defined as the distance between the locations of this maximum and this minimum. For purposes of finding the propagation speed of the shock, the “shock location” is defined as the steepest region of temperature decrease, which is located near the center of the shock region (shaded regions in Fig. 1).

## III Results

### III.1 Shock profile during formation and propagation

Both simulations start with a slightly randomized flat layer near the plate which is then oscillated. A periodic state is reached after a transient state of several cycles of the plate. Although most of the particles move together as a high-density layer, some particles are always found above and below the high-density region. We define “the layer” as the region in which the volume fraction is greater than $`4\%`$ of the random close-packed volume fraction $`\nu _{max}`$. Below this volume fraction, there are too few particles for accurate averaging in the MD simulation.

During each oscillation, the shock forms and propagates at the same phase angle with respect to the plate’s oscillation, so the dynamics of the system are fully characterized by the time interval between time $`ft=0`$ and one cycle later, $`ft=1.`$ We now investigate the formation and propagation of the shock by examining the behaviors found in MD and continuum simulations for various times $`ft`$ during the cycle (Fig. 1).

#### III.1.1 $`ft=0`$: The layer falls towards the plate

At $`ft=0`$ the container is at its minimum height. The particles, thrown off the plate in the previous cycle, are now falling towards the plate. The temperature is nearly zero for most of the material, increasing with height only in the low density region above the layer (Fig. 1). The effects of the artificial dissipation which is used in the continuum simulation are most pronounced at this time in the cycle, while the layer is dilated and above the plate. The shock formed in the previous cycle has propagated out through the top of the layer and is not present in this picture. The difference between the MD and continuum simulation at large heights has little effect on the shock behavior at later times in the cycle.

#### III.1.2 $`ft=0.15`$: Shock forms as the layer strikes plate

The layer now begins to strike the rising plate. The layer compresses and particle velocities are randomized by collisions, increasing $`\nu `$ and $`T`$ drastically near the plate. This results in the formation of a shock region where $`\nu `$ and $`T`$ change rapidly from the compressed, high temperature region near the plate to the dilute, low temperature undisturbed region which is still falling towards the plate. At $`ft=0.15`$ the shock is just beginning to form and is broad and not fully developed (see Sec. III C). During this time and throughout the formation and propagation of the shock, the shock profile from the continuum model shows remarkable agreement with the profile found in the MD simulation.

#### III.1.3 $`ft=0.22`$: Shock propagates through the layer

The layer continues to compress on the plate, resulting in the propagation of the shock up through the layer. Here the shock is fully developed and has propagated through much of the layer, but there is still a significant number of particles falling towards the plate. The shock has steepened since $`ft=0.15`$, with $`\nu `$ and $`T`$ changing from the upstream to the downstream values within a distance of about $`2\sigma .`$ As will be discussed in Sec. III C, the width of the shock continues to change throughout the cycle as the shock propagates through parts of the layer with different densities. A discontinuity in the derivatives of $`\nu `$ and $`T`$ appears at the leading edge of the shock, showing the boundary of the undisturbed region. Collisions cause the layer to gradually cool behind the shock, creating a lower temperature near the plate.

#### III.1.4 $`ft=0.40`$: The layer begins to leave the plate

The shock has propagated through the bulk of the layer and now enters the very dilute region above the layer. In this region the temperature becomes uncertain in the MD simulation because there are too few particles to use in averaging. Simultaneously, the plate is approaching its maximum height, and the layer begins to leave the plate as the downward acceleration of the plate becomes larger than $`g`$. The layer continues to cool, setting the stage for the next oscillation.

### III.2 Shock dynamics

The formation and propagation of the shock is shown in Fig. 2 for three values of the restitution coefficient $`e`$. The effect of varying inelasticity will be discussed in Sec. III D. In each case, when the bottom of the falling layer hits the plate, the center of mass of the layer is still falling towards the plate, and the layer is compressed. At this time, a shock is formed near the plate. The shock propagates up through the center of mass and out into the low density region above the layer. Then the layer leaves the plate, dilates and cools during its free flight, and once again falls towards the plate. As the shock propagates through the layer, it goes through regions of differing $`\nu `$ and $`T`$, and the shock velocity shows some slight variation with height, as seen in the curvature of the shock location (Fig. 2).

### III.3 Mach number and shock width

A prerequisite for shock formation is that the local Mach number of the flow be greater than unity with respect to the object causing the disturbance. We calculate the speed of sound from MD and continuum simulations using a relation derived from the continuum equation of state Eq. (7) Savage (1988):

$$c=\sqrt{T\chi \left(1+\frac{2}{3}\chi +\frac{\nu }{\chi }\frac{\chi }{\nu }\right)},$$
(14)
where $`\chi =1+2(1+e)G(\nu ).`$

As the layer falls towards the plate, it is accelerated by gravity and becomes supersonic in the reference frame of the plate. When the layer hits the plate, it is stopped by the impenetrable boundary, and the Mach number falls to nearly zero with respect to the plate (Fig. 3). The undisturbed region is still falling towards the plate at supersonic speeds, and the shock is formed at the boundary between the compressed (low $`Ma`$) region and the undisturbed (high $`Ma`$) region. The Mach number in the undisturbed region is largest near the shock where the layer has been falling the furthest under gravity, and gradually decreases with increasing height. The Mach number is directly calculated from the values of $`\nu ,`$ $`𝐮,`$ and $`T,`$ so the apparent increased discrepancy between the two simulations as compared to the discrepancy in previous figures is due to propagation of errors in these quantities in calculating $`Ma.`$

While the shock is in the layer, the Mach number always reaches its maximum at the lower boundary of the undisturbed region. Thus the maximum value of $`Ma`$ at a given time represents the value at the leading edge of the shock. The flows in this system are hypersonic, with $`Ma>10`$ just above the plate immediately before the shock is produced (Fig. 4). As the shock moves through the layer, the maximum $`Ma`$ decreases since the fastest flow has been stopped by the shock, but the flow just above the shock remains supersonic until after the shock leaves the layer at time $`ft=0.35`$, and $`Ma`$ always changes from supersonic to subsonic as the shock is crossed moving towards the plate.

Since the density of the layer varies in space and time, the mean free path of the particles in the layer changes as a function of height and time. The mean free path $`\xi `$ may be estimated as Huang (1987):

$$\xi =\frac{\sqrt{\pi /8}}{n\left(\pi \sigma ^2\right)}.$$
(15)
From $`ft=0.12`$ to $`ft=0.17`$ the shock is forming but is not fully developed. The dilute bottom of the layer hits the plate first at $`ft=0.12`$, so that while the initial shock is forming, a higher density region is still falling toward the plate. As higher and higher density material hits the plate, the shock width $`d`$ broadens (Fig. 4) as it forms from $`ft=0.12`$ to $`ft=0.17`$. After $`ft=0.17,`$ the shock steepens until the fully formed shock is about one mean free path in width, where it remains until the shock leaves the layer at $`ft=0.35.`$ After the shock leaves the layer, $`d`$ remains approximately one mean free path, but the mean free path itself increases in the low density region above the layer, leading to a broadening of the shock front as the shock leaves the layer (compare $`ft=0.40`$ to $`ft=0.22`$ in Fig. 1).

### III.4 Effect of inelasticity

Previous studies have investigated shocks in elastic and inelastic materials in the asymptotic limit of infinite time Goldshtein et al. (1996) and in inelastic materials at finite times Kamenetsky et al. (2000) for a piston moving into a uniform granular gas at a constant speed in the absence of gravity. We wish to investigate the effect of inelasticity on shock dynamics at finite times in the granular shaker system. However, the oscillatory state of the layer is not appropriate for comparing a gas of elastic hard spheres to one of inelastic hard spheres. An elastic gas, which has no mechanism for removing energy, will never reach an oscillatory state since the container injects energy into the layer at each cycle of the plate. In order to investigate the effect of inelasticity on shocks in granular materials, we slightly modify the system to look at the propagation of a single shock from identical initial conditions rather than looking at oscillatory behavior.

We run simulations with different coefficients of restitution $`e`$, starting from identical initial conditions of a dilated layer falling towards the plate. These initial conditions are taken from the oscillatory state of the system for $`e=0.9`$ using the same parameters as those used in previous sections. At time $`ft=0.33`$, when the layer is near the peak of its flight above the plate, the coefficient of restitution of the particles is suddenly increased or decreased.

Changing $`e`$ significantly changes properties of the shock and the layer following the next collision with the plate (Fig. 2). Since the initial conditions are the same for each case and the layer is in free fall, very few differences result from the change in $`e`$ until the layer collides with the plate. When the layer hits the plate it is compacted and many collisions occur between particles. The lower the coefficient of restitution, the more energy the particles lose through these collisions, and the more compact the layer becomes. In the case of high $`e`$, the layer rebounds and dilates quickly after the collision, while in the case of lower $`e`$, the layer remains very compact and leaves the plate almost as a solid body. The center of mass of the high $`e`$ layer flies very high before beginning to fall back towards the plate, while the lower $`e`$ value leads to a more shallow flight, with the center of mass beginning to fall earlier in the cycle. The shock trajectory exhibits curvature similar to that of the center of mass, remaining mostly straight for high $`e,`$ but bending downward significantly for low $`e.`$

After the collision with the plate, we track the propagation of the shock through the layer with a range of new values of $`e.`$ Starting from the same initial conditions, the shock starts at the same time and height regardless of $`e.`$ However, the velocity of the shock depends on $`e`$, as Fig. 5 shows. Although the shock speed changes somewhat as it propagates through the layer (see Fig. 2) , the average shock speed $`v_{shock}`$ increases monotonically with $`e`$ (see Fig. 5). The special case of elastic particles appears to match with the limit $`e1`$; thus there is no qualitative difference between elastic and inelastic gases.

Throughout this paper, both MD and continuum simulations were based on a collision model in which $`e`$ was independent of impact velocity. A more realistic hard-sphere collision model could include a coefficient of restitution $`e=e\left(v_n\right)`$ that depends on the normal component of the relative colliding velocity of the particles, $`v_n`$, such that the collisions become elastic as $`v_n`$ approaches zero. To test the effect of such a velocity-dependent restituition coefficient, we ran MD simulations including $`e\left(v_n\right)`$ of the same form as that used in studies of pattern formation in vertically oscillated layers Bizon et al. (1998): $`e\left(v_n\right)=10.1\left(v_n/\sqrt{g\sigma }\right)^{3/4}`$ for $`v_n<\sqrt{g\sigma }`$, and $`e=0.9`$ otherwise. The velocity-dependence of the coefficient of restitution produces negligible changes in the shock profiles and trajectories, and less than a $`1\%`$ change in the average velocity of the shock (Fig. 6). For the range of restitution values used in this paper ($`0.8e1.0`$), very few particles collide with low enough relative velocity for the velocity-dependence to significantly affect the outcome. However, for velocity independent $`e0.8`$, the MD simulation was unsuccessful due to inelastic collapse, while velocity dependent $`e\left(v_n\right)`$ may be used to avoid inelastic collapse with more inelastic particles.

### III.5 Shocks in layers of different depths

Thus far, we have restricted our investigation of shocks in oscillated layers to layer depths of $`H=9\sigma `$. We now investigate how changing $`H`$ affects shock formation and propagation for particles with $`e=0.9`$.

For fluidized layers, increasing the layer depth increases the pressure compressing the layer on the plate during collision with the plate, resulting in tighter packing of the particles near the plate. For $`4\sigma H15\sigma `$, increasing the layer depth increases the volume fraction of the layer near the plate, so that the trajectory of the center of mass of a layer of depth $`H=13.5\sigma `$ is nearly identical to that for $`H=9\sigma `$, even though the former has 1.5 times as many particles in the layer (Fig. 7). Both trajectories are actually lower than the center of mass trajectory for $`H=4.5\sigma `$.

For $`H=4.5\sigma `$, the layer is almost entirely gaseous, so much so that many particles reach heights greater than $`z=80\sigma `$ in each oscillation of the plate in MD simulations with a ceiling at $`z=200\sigma `$. Therefore, we use a cell of height $`200\sigma `$ for continuum simulations of $`H=4.5\sigma `$ (see Sec. II B).

Changing layer depth also affects the formation and propagation of shocks in the layer. For more dilute layers, the shock forms earlier in the cycle, when the lower edge of the dilute layer collides with the plate. Thus, increasing layer depth causes the shock to form later in the cycle, at which point the shock then propagates faster through the higher densities found in deeper layers (Fig. 7).

For $`H=13.5\sigma `$, the pressure on the layer during collision with the plate is so high that in some small regions near the plate, particles begin to form ordered lattices and the volume fraction increases beyond the random close-packed volume fraction. For deep layers $`H15\sigma `$, a significant portion of the layer is compressed into a solid lattice. The equation of state Eq. (7) does not allow volume fractions $`\nu >\nu _{max}`$, so continuum simulations cannot reproduce densities greater than random close-packed found in $`H15\sigma `$. For $`H4\sigma `$, the particles do not move as a coherent layer, but rather act as a gas filling the container for the entire cycle of the plate.

## IV Conclusions

Shocks form in the vertically oscillated granular system with each collision of the layer with the plate. Properties of these shocks captured in both MD and continuum simulations show good agreement even though the system exhibits a strong time dependence and large spatial gradients. Results from the continuum simulation agree with the MD results even for ranges where the validity of the continuum equations is in question, including when $`e`$ is significantly less than one ($`e=0.8`$) and when the volume fraction approaches close-packed. We also found that deeper layers exhibit denser packing near the plate and higher shock speeds than shallow layers. In addition, adding a simple velocity-dependence to the coefficient of restitution did not significantly change the behavior of the system. Finally, we found that increasing the coefficient of restitution causes no singular behavior in the shock velocity as $`e`$ approaches unity.

These results demonstrate that continuum equations can describe even a complicated three-dimensional, time-dependent granular system. These results also suggest directions for further research. A more realistic model of particle collisions could include friction between particles. The simulations should be further tested by direct comparison to laboratory measurements on shocks. Finally, an investigation of the role of a shock wave in transporting energy from the plate to the layer could provide useful information about the mechanism by which oscillations of the plate provide the energy to drive flows and patterns.

###### Acknowledgements.

We thank Daniel I. Goldman, W. D. McCormick, and Erin C. Rericha for helpful discussions. This work was supported by the Engineering Research Program of the Office of Basic Energy Sciences of the Department of Energy (Grant DE-FG03-93ER14312) and a grant from the Texas Advanced Research Program (Grant ARP 003658-0055-2001).


text:
# 1 Introduction

## 1 Introduction

Supersymmetry (SUSY) extension is one of the most promising ways to provide a solution to the gauge hierarchy problem in the standard model . The minimal version of this extension is called the minimal supersymmetric standard model (MSSM). Interestingly, the experimental data support the unification of the three gauge couplings at the scale $`M_{GUT}2\times 10^{16}\text{GeV}`$ with the MSSM particle contents . At high energies, our world may be described by a SUSY grand unified theory (GUT) with a simple gauge group such as SU(5) or SO(10) into which all the gauge groups in the standard model are embedded and unified.

On the other hand, much information about quark and lepton mass matrices has been accumulated in these decades. Especially, there exist, at present, strong evidences of neutrino masses and mixings through the interpretation of the (active) neutrino oscillations as the solutions to the solar neutrino deficit  and the atmospheric neutrino anomaly . According to the results, the MSSM has to be extended so as to incorporate neutrino masses and mixings.

In these points of view, the supersymmetric SO(10) GUT model is one of the most well-motivated theories, since the right-handed neutrinos are naturally incorporated and unified into the same multiplet, 16, together with the ordinary matters in the MSSM. Furthermore, this model has an additional advantage that the smallness of the neutrino masses can be naturally explained through the see-saw mechanism  with large right-handed Majorana neutrino masses.

There are many possibilities for the Higgs sector that will make the model realistic: namely, it correctly breaks the SO(10) gauge symmetry into the standard model one, reproduces observed fermion masses, mixings and CP-phase, and realizes the doublet-triplet Higgs mass splitting (in terminology of SU(5) GUT) etc. In this paper, we concentrate our discussion on the fermion masses and mixings, and consider the so-called “minimal” SO(10) model where only one 10 and one $`\overline{\mathrm{𝟏𝟐𝟔}}`$ Higgs multiplets have Yukawa couplings with quarks and leptons. It is very interesting to investigate this model, since it can fix all the fermion mass matrices including, at present, ambiguous and undetermined ones.

The minimal model was first seriously considered by Babu and Mohapatra  (non-SUSY case) and Lee and Mohapatra  (SUSY case), where the neutrino oscillation parameters were predicted with inputs of six quark masses, three mixing angles in the Cabibbo-Kobayashi-Maskawa (CKM) matrix and three charged-lepton masses. Unfortunately, the predictions were revealed to contradict the current experimental data. For example, the predicted mixing angle, $`\mathrm{sin}^22\theta _{23}`$, in the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix was found to be too small to be consistent with the the atmospheric neutrino oscillation data . It seems to be inevitable that the minimal model should be extended by introducing new Higgs multiplets so as to incorporate the realistic neutrino oscillation parameters.

However, note that, in Refs.  and , only the CP invariant case was analyzed. It was recently pointed out  that, once CP-violating phases are taken into account in the minimal model, there exits the parameter region which cannot be excluded by the current neutrino oscillation data. This observation brings a new hope to obtain the minimal SO(10) model compatible with all the observed data of fermion mass matrices. Nevertheless, more elaborate studies are needed, since, in , the GUT relation among the charged fermion mass matrices was applied to the mass matrices at the electroweak scale. This treatment is incomplete because the GUT relation is valid only at the GUT scale. In fact, we can see that Yukawa couplings evolve according to the renormalization group equations (RGEs), and the GUT relation is lost after the running in the “long desert” between the GUT scale and the electroweak scale.

In this paper, we pursue this program in more correct and complete way. Our strategy is the following. First, we evaluate the data of charged fermion mass matrices at the GUT scale by extrapolating the data at the electroweak scale according to their Yukawa coupling REGs. Next, these GUT scale data are applied to the GUT relation that is generalized to include CP-violating phases. This leads us to the explicit form of neutrino mass matrix at the GUT scale via the GUT relation. Lastly, running it back to the electroweak scale according to the RGE for the effective dimension-five operator , we compare our result with the current neutrino oscillation data.

This paper is organized as follows: in the next section, we give our basic framework of the minimal SUSY SO(10) model with one 10 and one $`\overline{\mathrm{𝟏𝟐𝟔}}`$ Higgs multiplets, and discuss the GUT relation among the fermion mass matrices. In Sec. 3, RGEs we use in our analysis are introduced. In Sec. 4, numerical analysis is performed, and the results are presented. Other predictions from our resultant mass matrices are discussed in Sec. 5. The last section is devoted to summary and comments.

## 2 Minimal SO(10) model and fermion mass matrices

We consider the minimal SUSY SO(10) model with a pair of 10+$`\overline{\mathrm{𝟏𝟐𝟔}}`$ Higgs multiplets, only which have Yukawa couplings (superpotential) with the matter multiplets such as

$`W_Y=Y_{10}^{ij}\mathrm{𝟏𝟔}_iH_{10}\mathrm{𝟏𝟔}_j+Y_{126}^{ij}\mathrm{𝟏𝟔}_iH_{126}\mathrm{𝟏𝟔}_j,`$ (1)
where $`\mathrm{𝟏𝟔}_i`$ is the matter multiplet of the i-th generation, $`H_{10}`$ and $`H_{126}`$ are the Higgs multiplet of 10 and $`\overline{\mathrm{𝟏𝟐𝟔}}`$ representations under SO(10), respectively. Note that, by virtue of the gauge symmetry, the Yukawa couplings, $`Y_{10}`$ and $`Y_{126}`$, are complex symmetric $`3\times 3`$ matrices. In our model, these Yukawa couplings are assumed to be real as in the original model  to keep the number of free parameters minimum.

In the SUSY GUT scenario, the GUT gauge symmetry is broken at the GUT scale, $`M_{GUT}10^{16}`$ GeV, into the standard model one. There are some different breaking chains suitable for our aim. Here, suppose the intermediate pass through the Pati-Salam subgroup , $`G_{422}=SU(4)_c\times SU(2)_L\times SU(2)_R`$, for simplicity. Under this gauge symmetry, the above Higgs multiplets are decomposed as $`\mathrm{𝟏𝟎}(\mathrm{𝟔},\mathrm{𝟏},\mathrm{𝟏})+(\mathrm{𝟏},\mathrm{𝟐},\mathrm{𝟐})`$ and $`\overline{\mathrm{𝟏𝟐𝟔}}(\mathrm{𝟔},\mathrm{𝟏},\mathrm{𝟏})+(\mathrm{𝟏𝟎},\mathrm{𝟑},\mathrm{𝟏})+(\overline{\mathrm{𝟏𝟎}},\mathrm{𝟏},\mathrm{𝟑})+(\mathrm{𝟏𝟓},\mathrm{𝟐},\mathrm{𝟐})`$, while $`\mathrm{𝟏𝟔}(\mathrm{𝟒},\mathrm{𝟐},\mathrm{𝟏})+(\overline{\mathrm{𝟒}},\mathrm{𝟏},\mathrm{𝟐})`$ for the matter multiplet. Breaking down to the standard model gauge group, $`SU(4)_c\times SU(2)_RSU(3)_c\times U(1)_Y`$, is accomplished by vacuum expectation value (VEV) of the $`(\overline{\mathrm{𝟏𝟎}},\mathrm{𝟏},\mathrm{𝟑})`$ Higgs multiplet. Note that the Majorana masses for the right-handed neutrinos are also generated by this VEV through the Yukawa coupling $`Y_{126}`$ in Eq. (1).

After this symmetry breaking, we find two pair of Higgs doublets which are in the same representation as the pair in the MSSM. One pair comes from $`(\mathrm{𝟏},\mathrm{𝟐},\mathrm{𝟐})\mathrm{𝟏𝟎}`$ and the other comes from $`(\overline{\mathrm{𝟏𝟓}},\mathrm{𝟐},\mathrm{𝟐})\overline{\mathrm{𝟏𝟐𝟔}}`$. Using these two pairs of the Higgs doublets, the Yukawa couplings of Eq. (1) are rewritten as

$`W_Y`$ $`=`$ $`\overline{u}_i\left(Y_{10}^{ij}H_{10}^u+Y_{126}^{ij}H_{126}^u\right)q_j+\overline{d}_i\left(Y_{10}^{ij}H_{10}^d+Y_{126}^{ij}H_{126}^d\right)q_j`$ (2)
$`+`$ $`\overline{\nu }_i\left(Y_{10}^{ij}H_{10}^u3Y_{126}^{ij}H_{126}^u\right)\mathrm{}_j+\overline{e}_i\left(Y_{10}^{ij}H_{10}^d3Y_{126}^{ij}H_{126}^d\right)\mathrm{}_j`$
$`+`$ $`\overline{\nu }_i\left(Y_{126}^{ij}v_R\right)\overline{\nu }_j,`$
where $`\overline{u}`$, $`\overline{d}`$, $`\overline{\nu }`$ and $`\overline{e}`$ are the right-handed $`SU(2)_L`$ singlet quark and lepton superfields, $`q`$ and $`\mathrm{}`$ are the left-handed $`SU(2)_L`$ doublet quark and lepton superfields, $`H_{10}^{u,d}`$ and $`H_{126}^{u,d}`$ are up-type and down-type Higgs doublet superfields originated from $`H_{10}`$ and $`H_{126}`$, respectively, and the last term is the Majorana mass term of the right-handed neutrinos with VEV of the $`(\overline{\mathrm{𝟏𝟎}},\mathrm{𝟏},\mathrm{𝟑})`$ Higgs, $`v_R`$. Quarks and leptons acquire Dirac masses through VEVs of these Higgs doublets. Note that, in that case, the entry of the Clebsch-Gordan coefficient, $`3`$, in the lepton sector plays the crucial role so that the unwanted GUT relations, $`m_d=m_e`$ and $`m_s=m_\mu `$, are corrected  in the same manner discussed by Georgi and Jarlskog .

Here, remember that the gauge coupling unification succeeds with only the MSSM particle contents. This means that only one pair of Higgs doublets remains light and the others should be heavy ($`M_G`$). It is highly non-trivial to accomplish this task naturally, while keeping the crucial role of $`H_{126}^{u,d}`$ to the fermion masses. This problem is essentially the same problem as the doublet-triplet Higgs mass splitting problem, and is referred as the “doublet-doublet splitting problem”. In this paper, we assume that this splitting is realized.

Now let us see the low energy superpotential of our model. After realizing the doublet-doublet Higgs mass splitting, we obtain only one pair of light Higgs doublets ($`H_u`$ and $`H_d`$) in the MSSM, which are, in general, admixture of all Higgs doublets having the same quantum numbers in the original model such as

$`H_u`$ $`=`$ $`\stackrel{~}{\alpha }_uH_{10}^u+\stackrel{~}{\beta }_uH_{126}^u+\mathrm{}`$
$`H_d`$ $`=`$ $`\stackrel{~}{\alpha }_dH_{10}^d+\stackrel{~}{\beta }_dH_{126}^d+\mathrm{},`$ (3)
where “$`+\mathrm{}`$” denotes the mixing with other Higgs doublets in the model (about which we do not discuss), $`\stackrel{~}{\alpha }_{u,d}`$ and $`\stackrel{~}{\beta }_{u,d}`$ denote elements of the unitary matrix which rotate the flavor basis in the original model into the (SUSY) mass eigenstates. Omitting the heavy Higgs mass eigenstates, the low energy superpotential is described by only the light Higgs doublets $`H_u`$ and $`H_d`$ such that

$`W_Y`$ $`=`$ $`\overline{u}_i\left(\alpha ^uY_{10}^{ij}+\beta ^uY_{126}^{ij}\right)H_uq_j+\overline{d}_i\left(\alpha ^dY_{10}^{ij}+\beta ^dY_{126}^{ij}\right)H_dq_j`$ (4)
$`+`$ $`\overline{\nu }_i\left(\alpha ^uY_{10}^{ij}3\beta ^uY_{126}^{ij}\right)H_u\mathrm{}_j+\overline{e}_i\left(\alpha ^dY_{10}^{ij}3\beta ^dY_{126}^{ij}\right)H_d\mathrm{}_j`$
$`+`$ $`\overline{\nu }_i\left(Y_{126}^{ij}v_R\right)\overline{\nu }_j,`$
where the formulas of the inverse unitary transformation of Eq. (3), $`H_{10}^{u,d}=\alpha ^{u,d}H_{u,d}+\mathrm{}`$ and $`H_{126}^{u,d}=\beta ^{u,d}H_{u,d}+\mathrm{}`$, have been used. Note that the elements of the unitary matrix, $`\alpha ^{u,d}`$ and $`\beta ^{u,d}`$, are in general complex parameters, through which CP-violating phases are introduced into the fermion mass matrices.

Providing the Higgs VEVs, $`H_u=v\mathrm{sin}\beta `$ and $`H_d=v\mathrm{cos}\beta `$ with $`v=174\text{GeV}`$, the quark and lepton mass matrices can be read off as

$`M_u`$ $`=`$ $`c_{10}M_{10}+c_{126}M_{126}`$
$`M_d`$ $`=`$ $`M_{10}+M_{126}`$
$`M_D`$ $`=`$ $`c_{10}M_{10}3c_{126}M_{126}`$
$`M_e`$ $`=`$ $`M_{10}3M_{126}`$
$`M_R`$ $`=`$ $`c_RM_{126},`$ (5)
where the mass matrices are defined as $`M_{10}=Y_{10}\alpha ^dv\mathrm{cos}\beta `$ and $`M_{126}=Y_{126}\beta ^dv\mathrm{cos}\beta `$, respectively, $`c_{10}=(\alpha ^u/\alpha ^d)\mathrm{tan}\beta `$, $`c_{126}=(\beta ^u/\beta ^d)\mathrm{tan}\beta `$, and $`c_R=v_R/(\beta ^dv\mathrm{cos}\beta )`$. All the quark and lepton mass matrices are characterized by only two basic mass matrices, $`M_{10}`$ and $`M_{126}`$, and three coefficients. This means the strong predictability of the model to the fermion mass matrices in the MSSM.

Now let us count how many free parameters are included in the above mass matrices, not counting $`c_R`$. Remember that we began with real symmetric matrices $`Y_{10}`$ and $`Y_{126}`$. By flavor rotation of the GUT matter multiplets, we can begin with the basis where $`M_{10}`$ is real and diagonal, which includes three free parameters. There are seven parameters in $`M_{126}`$: six parameters originate from the symmetric real matrix $`Y_{126}`$, and one from the phase of $`\beta ^d`$. Adding the four parameters from complex $`c_{10}`$ and $`c_{126}`$, we find 14 free parameters in total, which are evaluated as six quark masses, three angles and one CP-phase in the CKM matrix and three charged-lepton masses, but one parameter is left free . Now two free parameters are newly introduced in addition to 12 free parameters analyzed in  and , by which our model is generalized to include CP-violating phases.

The parameter, $`c_R`$, appears in the right-handed Majorana neutrino mass matrix, by which masses of the right-handed Majorana neutrinos are determined after $`M_{126}`$ are fixed. Since $`v_RM_G`$, $`c_Rv_R/v10^{14}\text{GeV}`$ is naturally expected. However, we treat $`c_R`$ as a free parameter in the following analysis, since the Higgs multiplet breaking $`G_{422}`$, in general, mixes with other Higgs multiplets being in the same representation under $`G_{422}`$. Also there should exist Higgs multiplets which also have $`G_{422}`$ breaking VEV to cancel the D-term VEV and ensure SUSY. In addition, threshold corrections would introduce an uncertainty of $`v_R`$.

The mass matrix formulas in Eq.(5) leads to the GUT relation among the quark and lepton mass matrices,

$`M_e=c_d\left(M_d+\kappa M_u\right),`$ (6)
where

$`c_d`$ $`=`$ $`{\displaystyle \frac{3c_{10}+c_{126}}{c_{10}c_{126}}}`$
$`\kappa `$ $`=`$ $`{\displaystyle \frac{4}{3c_{10}+c_{126}}}.`$ (7)
Without loss of generality, we can begin with the basis where $`M_u`$ is real and diagonal, $`M_u=D_u`$. Since $`M_d`$ is the symmetric matrix, it can be described as $`M_d=V_{KM}^{}D_dV_{KM}^{}`$ by using the CKM matrix $`V_{KM}`$ and the real diagonal mass matrix $`D_d`$. Considering the basis-independent quantities, $`\text{tr}[M_e^{}M_e]`$, $`\text{tr}[(M_e^{}M_e)^2]`$ and $`\text{det}[M_e^{}M_e]`$, and eliminating $`|c_d|`$, we obtain two independent equations such that

$`\left({\displaystyle \frac{\text{tr}[\stackrel{~}{M_e}^{}\stackrel{~}{M_e}]}{m_e^2+m_\mu ^2+m_\tau ^2}}\right)^2`$ $`=`$ $`{\displaystyle \frac{\text{tr}[(\stackrel{~}{M_e}^{}\stackrel{~}{M_e})^2]}{m_e^4+m_\mu ^4+m_\tau ^2}},`$ (8)
$`\left({\displaystyle \frac{\text{tr}[\stackrel{~}{M_e}^{}\stackrel{~}{M_e}]}{m_e^2+m_\mu ^2+m_\tau ^2}}\right)^3`$ $`=`$ $`{\displaystyle \frac{\text{det}[\stackrel{~}{M_e}^{}\stackrel{~}{M_e}]}{m_e^2m_\mu ^2m_\tau ^2}},`$ (9)
where $`\stackrel{~}{M_e}=V_{KM}^{}D_dV_{KM}^{}+\kappa D_u`$. With the input data of six quark masses, three angles and one CP-phase in the CKM matrix and three charged-lepton masses, we can determine $`\kappa `$ and $`|c_d|`$, but one parameter, the phase of $`c_d`$, is left undetermined . Now the mass matrices, $`M_{10}`$ and $`M_{126}`$, are described by

$`M_{10}`$ $`=`$ $`{\displaystyle \frac{3+|c_d|e^{i\sigma }}{4}}V_{KM}^{}D_dV_{KM}^{}+{\displaystyle \frac{|c_d|e^{i\sigma }\kappa }{4}}D_u,`$ (10)
$`M_{126}`$ $`=`$ $`{\displaystyle \frac{1|c_d|e^{i\sigma }}{4}}V_{KM}^{}D_dV_{KM}^{}{\displaystyle \frac{|c_d|e^{i\sigma }\kappa }{4}}D_u,`$ (11)
as the functions of $`\sigma `$, the phase of $`c_d`$, with the solutions, $`|c_d|`$ and $`\kappa `$, determined by the GUT relation.

Note that the GUT relation of Eq. (6) is valid only at the GUT scale. Therefore, the input data needed to determine $`\kappa `$ and $`|c_d|`$ should be the data at the GUT scale. To evaluate the GUT scale data, RGE analysis of Yukawa couplings is necessary. This is the subject of the next section.

## 3 One-loop RGEs

The data we need at the GUT scale are six quark masses, three angles and one CP-phase in the CKM matrix and three charged-lepton masses. In general, it is not a straightforward task to obtain the data at GUT scale by extrapolating the data at the electroweak scale. This is because RGEs of quark and charged-lepton Yukawa couplings include unknown neutrino Yukawa couplings. We can obtain a reliable answer only if neutrino Yukawa couplings are negligible.

Fortunately, in models with the see-saw mechanism, the right-handed neutrinos are decoupled at the scale below the heavy Majorana masses, and the neutrino Yukawa couplings become irrelevant and disappear from RGEs of quark and charged-lepton Yukawa couplings. This fact makes our analysis straightforward, at least, below the heavy Majorana mass scale. In the following analysis, we use the RGEs independent of the neutrino Yukawa couplings at all scales between the electroweak and the GUT scales. Although this treatment fails at the scale above the right-handed Majorana neutrino masses, we expect that our results are reliable, since the heavy Majorana mass scale is close to the GUT scale. We will give a comment on this point in the last section.

We analyze one-loop RGEs for gauge couplings and quark and charged-lepton Yukawa couplings in the MSSM. For simplicity, the universal soft SUSY breaking masses are assumed, $`m_{SUSY}=M_Z91\text{GeV}`$. In this case, the GUT scale is found to be $`M_G2\times 10^{16}\text{GeV}`$. The RGEs used in our analysis are as follows . For gauge couplings,

$`{\displaystyle \frac{dg_i}{d(\mathrm{ln}\mu )}}={\displaystyle \frac{b_i}{16\pi ^2}}g_i^3,`$ (12)
where $`b_i=(3,1,33/5)`$ for $`i=3,2,1`$ corresponding to $`SU(3)_c`$, $`SU_L`$ and $`U(1)_Y`$ gauge couplings, respectively. We use the electroweak scale inputs such that $`\alpha _3(M_Z)=g_3^2/4\pi =0.118`$, $`\alpha _2(M_Z)=g_3^2/4\pi =0.0335`$ and $`\alpha _1(M_Z)=g_3^2/4\pi =0.0168`$. On the other hand, RGEs for the Yukawa couplings of up-type quarks ($`Y_u`$), down-type quarks ($`Y_d`$), and charged-leptons ($`Y_e`$) are found to be

$`{\displaystyle \frac{dS_F}{d(\mathrm{ln}\mu )}}={\displaystyle \frac{1}{16\pi ^2}}\left(\beta _FS_F+S_F\beta _F\right),`$ (13)
where $`S_F=Y_F^{}Y_F`$ with fermion index $`F=u,d,e`$, and

$`\beta _u`$ $`=`$ $`3S_u+S_d+\left\{\text{tr}\left[3S_u\right]\left({\displaystyle \frac{16}{3}}g_3^2+3g_2^2+{\displaystyle \frac{13}{15}}g_1^2\right)\right\}\mathrm{𝟏}_{3\times 3}`$
$`\beta _d`$ $`=`$ $`3S_d+S_u+\left\{\text{tr}\left[3S_d+S_e\right]\left({\displaystyle \frac{16}{3}}g_3^2+3g_2^2+{\displaystyle \frac{7}{15}}g_1^2\right)\right\}\mathrm{𝟏}_{3\times 3}`$
$`\beta _e`$ $`=`$ $`3S_e+\left\{\text{tr}\left[3S_d+S_e\right]\left(3g_2^2+{\displaystyle \frac{9}{5}}g_1^2\right)\right\}\mathrm{𝟏}_{3\times 3}.`$ (14)
To determine the input values of fermion Yukawa couplings at $`M_Z`$, we have to fix $`\mathrm{tan}\beta `$ by hand. For given $`\mathrm{tan}\beta `$, they are given by

$`S_u`$ $`=`$ $`\text{diag}({\displaystyle \frac{m_u^2}{v^2}},{\displaystyle \frac{m_c^2}{v^2}},{\displaystyle \frac{m_t^2}{v^2}})\left(1+\mathrm{tan}^2\beta \right),`$ (15)
$`S_d`$ $`=`$ $`V_{KM}\text{diag}({\displaystyle \frac{m_d^2}{v^2}},{\displaystyle \frac{m_s^2}{v^2}},{\displaystyle \frac{m_b^2}{v^2}})V_{KM}^{}\left(1+\mathrm{tan}^2\beta \right),`$ (16)
$`S_e`$ $`=`$ $`\text{diag}({\displaystyle \frac{m_e^2}{v^2}},{\displaystyle \frac{m_\mu ^2}{v^2}},{\displaystyle \frac{m_\tau ^2}{v^2}})\left(1+\mathrm{tan}^1\beta \right).`$ (17)
Here, we have worked out in the basis where $`S_u`$ and $`S_e`$ are diagonal. Given $`\mathrm{tan}\beta `$, it is easy to solve the above RGEs numerically. Although, as can be seen from the above RGEs, $`S_u`$ is no longer diagonal after running, we can extract the data we need at the GUT scale by the usual manner. For example, the CKM matrix at the GUT scale is given by $`V_{KM}=V_u^{}V_d`$, where $`V_u`$ and $`V_d`$ are the unitary matrices which diagonalize $`S_u`$ and $`S_d`$, respectively.

Solving the GUT relation as discussed in the previous section by using the extrapolated data, we can obtain the light Majorana neutrino mass matrix, $`M_\nu `$, as a function of $`\sigma `$ and $`c_R`$ through Eqs. (5), (10) and (11) and the see-saw mechanism, $`M_\nu =M_D^TM_R^1M_D`$. Since this is the mass matrix at the GUT scale, it should be extrapolated to the electroweak scale according to the RGE for the effective dimension-five operator  such that

$`{\displaystyle \frac{dM_\nu }{d(\mathrm{ln}\mu )}}={\displaystyle \frac{1}{16\pi ^2}}\left\{S_e^TM_\nu +M_\nu S_e+\left(6\text{tr}\left[S_u\right]6g_2^2{\displaystyle \frac{6}{5}}g_1^2\right)M_\nu \right\}.`$ (18)
After the running, our resultant neutrino mass matrix are compared with the neutrino oscillation data. In the next section, the effect of the RGE running will be found to be small and almost negligible in our model.

## 4 Numerical analysis and results

Now we are ready to perform all the numerical analysis. Fist of all, note that the solution of the GUT relation exists, only if we take appropriate input parameters. This is because the number of our free parameters (fourteen) are almost the same as the number of inputs (thirteen). Therefore, in the following analysis, we vary two input parameters, $`m_s`$ and CP-phase $`\delta `$ in the CKM matrix, within the experimental errors so as to find the solution. We take input fermion masses at $`M_Z`$ as follows (in GeV):

$`m_u=0.00233,m_c=0.677,m_t=176`$
$`m_d=0.00469,m_b=3.00`$
$`m_e=0.000487,m_\mu =0.103,m_\tau =1.75`$
Here, the experimental values extrapolated from low energies to $`M_Z`$ were used . For the CKM mixing angles in the “standard” parameterization, we input the center values measured by experiments as follows :

$`s_{12}={\displaystyle \frac{0.219+0.226}{2}},s_{23}={\displaystyle \frac{0.037+0.043}{2}},s_{13}={\displaystyle \frac{0.002+0.005}{2}}.`$
Since it is very difficult to search all the possible parameter region systematically, we present only reasonable results we found. In the following, we show our analysis in detail in the case $`\mathrm{tan}\beta =45`$. Here we take $`m_s=0.0729\text{(GeV)}`$ and $`\delta =86.4^{}`$. In this case, the CKM matrix at $`M_Z`$ is given by

$`V_{KM}(M_Z)=\left(\begin{array}{ccc}0.975& 0.222& 0.0002200.00349i\\ 0.2220.000136i& 0.9740.0000311i& 0.0400\\ 0.008690.00340i& 0.03900.000777i& 0.999\end{array}\right).`$ (22)
After RGE evolutions to the GUT scale, fermion masses (up to sign) and the CKM matrix are found as follows (in GeV) :

$`m_u=0.00103,m_c=0.301,m_t=134`$
$`m_d=0.00171,m_s=0.0265,m_b=1.56`$
$`m_e=0.000413,m_\mu =0.0872,m_\tau =1.69`$
and

$`V_{KM}(M_G)=\left(\begin{array}{ccc}0.975& 0.222& 0.0001750.00279i\\ 0.2220.000121i& 0.974+0.000129i& 0.0320\\ 0.006950.00272i& 0.03120.000626i& 0.999\end{array}\right)`$ (26)
in the standard parameterization.

These outputs at the GUT scale are used as input parameters in order to solve Eqs. (8) and (9). The contours of solutions of each equations are depicted in Fig. 1, where the signs of the input fermion masses have taken as $``$ for $`m_u`$, $`m_c`$, $`m_d`$, and $`m_s`$, and $`+`$ for $`m_t`$ and $`m_b`$. We find a solution

$`\kappa =0.0134+0.000786i,`$
$`|c_d|=6.38,`$ (27)
which leads to the reasonable results for the neutrino sector.<sup>1</sup><sup>1</sup>1 Although there are other solutions, we do not address them. They lead to small $`\mathrm{sin}^22\theta _{23}`$ inconsistent with the atmospheric neutrino oscillation data, and are out of our interest.

Now we can describe all the fermion mass matrices as functions of $`\sigma `$ and $`c_R`$ by using the above solution through the mass matrix formulas of Eqs. (5), (10) and (11). Since the PMNS mixing angles are independent of $`c_R`$, we can plot the angles as functions of $`\sigma `$. The results are depicted in Fig. 2, and we can see that the resultant mixing angles are very sensitive to $`\sigma `$. For $`\sigma =3.198[\text{rad}]`$, we obtain $`\mathrm{sin}^22\theta _{12}=0.722`$, $`\mathrm{sin}^22\theta _{23}=0.881`$ and $`\mathrm{sin}^22\theta _{13}=0.164`$ at the GUT scale. After running this results back to the electroweak scale according to RGE of Eq. (18), we find

$`\mathrm{sin}^22\theta _{12}=0.723,\mathrm{sin}^22\theta _{23}=0.895,\mathrm{sin}^22\theta _{13}=0.164.`$ (28)
Note that RGE running effects are almost negligible. The ratio $`\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{}^2`$ is also independent of $`c_R`$, where $`\mathrm{\Delta }m_{}^2`$ and $`\mathrm{\Delta }m_{}^2`$ are the oscillation parameters relevant for the solar and the atmospheric neutrino deficits, respectively. In Fig. 3, the ratio at the GUT scale is depicted as a function of $`\sigma `$. After the RGE running, we find, at the electroweak scale,

$`{\displaystyle \frac{\mathrm{\Delta }m_{}^2}{\mathrm{\Delta }m_{}^2}}=0.188.`$ (29)
RGE running effects are almost negligible also for the ratio. The neutrino mass matrix at the electroweak scale and the PMNS matrix which lead to the above results are as follows:

$`M_\nu =c_R^1\left(\begin{array}{ccc}14.7+48.0i& 46.996.5i& 446370i\\ 46.996.5i& 347+90.3i& 1030+702i\\ 446370i& 1030+702i& 497914i\end{array}\right)(\text{GeV}),`$ (33)
and

$`U_{PMNS}=\left(\begin{array}{ccc}0.168+0.838i& 0.467+0.0940i& 0.00508+0.207i\\ 0.0519+0.498i& 0.6510.0473i& 0.01890.569i\\ 0.0745+0.116i& 0.4500.381i& 0.431+0.669i\end{array}\right).`$ (37)
They have been calculated in the basis where the charged-lepton mass matrix is diagonal and all the matrix elements are real and positive. The neutrino mass eigenvalues are determined by fixing $`c_R`$ appropriately. For example, taking $`c_R=3.01\times 10^{13}`$ leads to $`\mathrm{\Delta }m_{}^2=3.76\times 10^4\text{eV}^2`$ and $`\mathrm{\Delta }m_{}^2=2.00\times 10^3\text{eV}^2`$. In this case, three mass eigenvalues of the right-handed neutrinos are found to be (in GeV) $`M_{R_1}=1.64\times 10^{11}`$, $`M_{R_2}=2.50\times 10^{12}`$ and $`M_{R_3}=8.22\times 10^{12}`$, respectively. These results are in agreement with the results of the recent neutrino oscillation analysis  with the large angle solution with Mikheev-Smirnov-Wolfenstein effect  for the solar neutrino deficit.

Our results are very sensitive to the input values of $`m_s`$ and $`\delta `$. As mentioned above, the solution of Eqs. (8) and (9) exists only if the input values of $`m_s`$ and $`\delta `$ are appropriately fixed. In Fig. 1, as $`m_s`$ is taken to be larger for fixed $`\delta `$, the contour of Eq.(8) moves to the left, while the circle moves to the right. Thus, the solution disappears eventually. On the other hand, as $`\delta `$ is taken to be larger for fixed $`m_s`$, both of the contours move to the left, but the circle moves faster. Eventually, the solution disappears. This disappearance of the solution occurs by several percent changes of the input parameters. Therefore, with fixed $`m_s`$ ($`\delta `$), $`\delta `$ ($`m_s`$) is almost determined in order for the solution to exist.

Next let us see a more important parameter dependence. When $`\delta `$ and $`m_s`$ are fixed appropriately and the solutions are found, we obtain the similar plot as in Fig. 2 for the mixing angles in the PMNS matrix. However, the hight of the peaks depends on the input $`\delta `$ values. As in the past works  , we mostly find that the resultant $`\mathrm{sin}^22\theta _{23}`$ is too small to be consistent with the experimental data. Only the special values of the inputs and $`\sigma `$ can give the results being in agreement with all the neutrino oscillation data. The reasonable results we found are as follows:

| $`\mathrm{tan}\beta `$ | $`m_s(M_Z)`$ | $`\delta `$ | $`\sigma `$ | $`\mathrm{sin}^22\theta _{12}`$ | $`\mathrm{sin}^22\theta _{23}`$ | $`\mathrm{sin}^22\theta _{13}`$ | $`\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{}^2`$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 40 | 0.0718 | $`93.6^{}`$ | 3.190 | 0.738 | 0.900 | 0.163 | 0.205 |
| 45 | 0.0729 | $`86.4^{}`$ | 3.198 | 0.723 | 0.895 | 0.164 | 0.188 |
| 50 | 0.0747 | $`77.4^{}`$ | 3.200 | 0.683 | 0.901 | 0.164 | 0.200 |
| 55 | 0.0800 | $`57.6^{}`$ | 3.201 | 0.638 | 0.878 | 0.152 | 0.198 |

The mass eigenvalues of the light Majorana neutrinos are determined by $`c_R`$. In other words, $`c_R`$ is a function of the oscillation parameter, and appropriate value of $`c_R`$ is fixed so as to be suitable for the neutrino oscillation data. These relations are listed in the following.

| $`\mathrm{tan}\beta `$ | $`c_R\times m_{\nu _3}\text{(TeV)}`$ | $`c_R\times \sqrt{\frac{\mathrm{\Delta }m_{}^2}{10^3\text{eV}^2}}`$ |
| --- | --- | --- |
| 40 | 2.18 | $`6.26\times 10^{13}`$ |
| 45 | 2.08 | $`6.02\times 10^{13}`$ |
| 50 | 1.81 | $`5.23\times 10^{13}`$ |
| 55 | 1.42 | $`4.04\times 10^{13}`$ |

Here, $`m_{\nu _3}`$ is the mass eigenvalue of the heaviest light Majorana neutrino. As mentioned above, our resultant neutrino oscillation parameters are sensitive to all the input parameters. In other words, if we use the neutrino oscillation data as the input parameters, the other input, for example, the CP-phase in the CKM matrix can be regarded as the prediction of our model. It is a very interesting observation that the CP-phases listed above are in the region consistent with experiments .

## 5 Other predictions

Now all the mass matrices have been completely determined in our model as discussed in the previous section. There are other interesting physical observables which can be calculated by using the concrete fermion mass matrices in our model.

The averaged neutrino mass relevant for the neutrino-less double beta decay  can be read off from $`(1,1)`$ component of the light Majorana neutrino mass matrix $`M_\nu `$. The CP-violation in the lepton sector is characterized by the Jarlskog parameter  defined as

$`J_{CP}=\text{Im}\left[U_{e2}U_{\mu 2}^{}U_{e3}^{}U_{\mu 3}\right],`$ (38)
where $`U_{fi}`$ is the PMNS matrix element.

It is well known that the SO(10) GUT model possesses a simple mechanism of baryogenesis through the out-of-equilibrium decay of the right-handed neutrinos, namely, the leptogenesis . The amount of the created baryon asymmetry is characterized by the CP-violating parameter $`ϵ`$ estimated as

$`ϵ{\displaystyle \underset{j=2,3}{}}{\displaystyle \frac{\text{Im}\left[(M_DM_D^{})_{1j}^2\right]}{v^2\mathrm{sin}^2\beta (M_DM_D^{})_{11}}}{\displaystyle \frac{M_{R_1}}{M_{R_j}}}.`$ (39)
Here, $`M_D`$ is the neutrino Dirac mass matrix in the basis where the right-handed Majorana mass matrix $`M_R`$ is real and diagonal, and $`M_{R_j}`$ is the mass eigenvalue of the right-handed Majorana neutrino of the i-th generation.

These quantities are evaluated by using the results presented in the previous section, and results are listed in the following.

| $`\mathrm{tan}\beta `$ | $`m_\nu _{ee}(\text{eV})`$ | $`J_{CP}`$ | $`ϵ`$ |
| --- | --- | --- | --- |
| 40 | 0.00122 | $`0.00110`$ | $`7.39\times 10^5`$ |
| 45 | 0.00118 | $`0.00429`$ | $`6.80\times 10^5`$ |
| 50 | 0.00119 | $`0.00631`$ | $`6.50\times 10^5`$ |
| 55 | 0.00117 | $`0.00612`$ | $`11.2\times 10^5`$ |

Here, $`c_R`$ was fixed so that $`\mathrm{\Delta }m_{}^2=2\times 10^3\text{eV}^2`$. Unfortunately, both of the averaged neutrino mass and the Jarlskog parameter may be too small to expect their evidences in future experiments. On the other hand, the CP-violating parameter is too large to be consistent with the observed baryon asymmetry. The problem of this too large baryon asymmetry can be easily avoided by assuming an inflationary universe whose reheating temperature is smaller than the right-handed neutrino masses. However, in this case, the leptogenesis scenario no longer works, and another scenario of baryogenesis such as the Affleck-Dine mechanism  may be applicable.

Sizable lepton-flavor violation (LFV) can be expected in SUSY models, and the LFV processes are one of the most important processes as the low-energy SUSY search . In detailed analysis, concrete fermion mass matrices are necessary. It is worth investigating the LFV processes in our model.

## 6 Summary and comments

We have discussed the minimal SUSY SO(10) model, where only one 10 and one $`\overline{\mathrm{𝟏𝟐𝟔}}`$ Higgs multiplets have Yukawa couplings with fermions. This model can determine all the fermion mass matrices with only a few free parameters. It is known that, in the absence of CP violation, this model cannot incorporate the realistic neutrino mass matrix consistent with the neutrino oscillation data. We examined more general case in which CP-violating phases in the fermion mass matrices were introduced. In this case, using the experimental data of six quark masses, three angles and one CP-phase in the CKM matrix and three charged-lepton masses, the light Majorana neutrino mass matrix can be determined as the function of only one free parameter $`\sigma `$ through the GUT relation and the see-saw mechanism. Here we do not count $`c_R`$, which plays only the role to fix the overall scale of the neutrino masses. In order to connect the mass matrix data between the electroweak and the GUT scales, we have analyzed the one-loop RGEs of the charged fermion Yukawa couplings and the effective dimension-five operator of the light Majorana neutrinos.

We found that there was the parameter region in which the predicted neutrino mass matrix can be consistent with the current neutrino oscillation data. This parameter region is severely constrained, and all the parameters, even the CP-phase in the CKM matrix, are almost fixed for given $`\mathrm{tan}\beta `$. In other words, we can regard the neutrino oscillation data as the inputs while the CP-phase as the output. Interestingly, our results consistent with the neutrino oscillation data are also consistent with experimental results of CP-violation in the quark sector.

In our RGE analysis, we have used RGEs introduced in Sec. 3 at all the energy scales between the electroweak scale and the GUT scale. This treatment fails at the energy larger than the right-handed Majorana neutrino mass scale, where the right-handed neutrinos are no longer decoupled. Above that scale, we have to take RGE of neutrino Yukawa coupling into account. However, as found in Sec. 4, the Majorana mass scale is about $`10^{13}\text{GeV}`$, and close to the GUT scale in the sense of the RGE (logarithmic) running. Thus, we expect that our results can be still reliable. In fact, we can find that changes of our input values at the GUT scale remain within several percents even if the neutrino Yukawa coupling RGE is taken into account, and that the results presented in Sec. 4 are almost unchanged.

In addition, we have neglected the SUSY threshold corrections in our RGE analysis. It is known that with large $`tan\beta `$ the down-type quark mass matrix is potentially affected by large SUSY threshold corrections and we can neglect them with a limited region of the soft SUSY breaking parameters . If the SUSY threshold corrections are large, it should be taken into account in our analysis, and will change the input values at the GUT scale and eventually the final result. It is an interesting question whether the result gives a better or worse fit, and it is worth analyzing the RGEs taking SUSY threshold corrections into account.

When the CHOOZ data  is concerned, our results seems to be in sever situation . However, more comprehensive analysis  gives rather loose constraint. In any case, our results lie in the allowed region at 99% C.L. in the global analysis , and we have concluded that they are in agreement with the current neutrino oscillation data.<sup>2</sup><sup>2</sup>2 We thank Osamu Yasuda for his analysis concluding that our results is 2.6 $`\sigma `$ away from the current best fit value.  This situation seems to be inevitable in our model, since two angles, $`\mathrm{sin}^22\theta _{23}`$ and $`\mathrm{sin}^22\theta _{13}`$, are strongly correlated and have the peaks at almost the same $`\sigma `$ value as can be seen in Fig. 2. This can be easily avoided, when we extend the minimal model and introduce new Higgs multiplets. An extended model has new free parameters in fermion mass matrices, and can be in much better agreement with the neutrino oscillation data, larger $`\mathrm{sin}^22\theta _{23}`$ but smaller $`\mathrm{sin}^22\theta _{13}`$. However, what we have found in this paper is that the minimal SO(10) model, the model with only minimal set of Higgs multiplets, is still viable without any extention.

We have concentrated our discussion only on the Yukawa sector with one $`\mathrm{𝟏𝟎}`$ and one $`\overline{\mathrm{𝟏𝟐𝟔}}`$ Higgs multiplets. In order to complete our model, it is necessary to construct a concrete Higgs sector which can realize all the assumptions in our model, the correct GUT symmetry breaking, doublet-doublet (triplet) Higgs mass splitting etc. Such a Higgs sector may include new Higgs multiplets which affect the fermion mass matrices. If it is the case, our model becomes less predictive, although the model can be in better agreement with the neutrino oscillation data.<sup>3</sup><sup>3</sup>3 In this point of view, the Higgs sector proposed in  is interesting, because it remains the Yukawa sector minimal. This issue is highly non-trivial, and a further work is needed.

Although, in this paper, the SUSY SO(10) model has been considered, non-SUSY SO(10) model with intermediate symmetry breaking scales is also worth investigating. In this case, the strategy is completely analogous to that we discussed in this paper, except that RGEs are replaced by the ones of non-SUSY case.

## Acknowledgments

We would like to thank Rabindra N. Mohapatra for useful discussions and comments. The author (N.O) would like to thank Siew-Phang Ng for his useful advice on our numerical analysis.


text:
# Scale-Free Networks are Ultrasmall

## ACKNOWLEDGMENTS

We would like to thank Daniel ben-Avraham and Tomer Kalisky for useful discussion.


text:
# The semiclassical limit of chaotic eigenfunctions

## Abstract

A generic chaotic eigenfunction has a non-universal contribution consisting of scars of short periodic orbits. This contribution, which can not be explained in terms of random universal waves, survives the semiclassical limit (when $`\mathrm{}`$ goes to zero). In this limit, the sum of scarred intensities is a simple function of $`\eta \sqrt{\pi /2}(f1)h_T^1(\lambda _i^2)^{1/2}`$, with $`f`$ the degrees of freedom, $`h_T`$ the topological entropy and $`\{\lambda _i\}`$ the set of positive Lyapunov exponents. Moreover, the fluctuations of this representation go to zero as $`1/|\mathrm{ln}\mathrm{}|`$. For this reasson, we will be able to provide a detailed description of a generic chaotic eigenfunction in the semiclassical limit.

Berry  and Voros  proposed a semiclassical description of chaotic eigenfunctions by considering the surface of constant energy as the unique classical invariant able to support them. This would be the case if the time required for the definition of individual eigenfunctions was infinite. However, the required Heisemberg time $`T_H`$ is finite for finite values of $`\mathrm{}`$, even though it goes to infinity in the semiclassical limit (when $`\mathrm{}`$ goes to zero).

Recently, we have derived a semiclassical theory which describes individual eigenfunctions in terms of short periodic orbits $`(\mathrm{𝐏𝐎𝐬})`$ . The number of $`\mathrm{𝐏𝐎𝐬}`$ is such that the sum of their periods is around $`T_H`$. Using these invariants, we will give a description of chaotic eigenfunctions with fluctuations going to zero with $`\mathrm{}`$.

It is worth to emphasize that the fluctuations of the Berry-Voros description go to infinity in the semiclassical limit (SL). But this statement demands explanation because it is common to find in the literature expresions like ” …this description is supported by Shnirelman theorem.”. For example, the peaks of the Husimi distribution of a chaotic eigenfunction have very strong fluctuations with respect to their classical ergodic measure; by assuming a random wave hypothesis they go logarithmically to infinity as $`\mathrm{}0`$ . However, each of the peaks occupy a volume (the semiclassical volume of a quantum state) going to zero with $`\mathrm{}`$. Then, a smoothing on a region independent of $`\mathrm{}`$, a classical smoothing, is sufficient for washing out fluctuations in the SL. In this context, Shinirelman theorem  may be expresed as follows: a classical smoothing is sufficient for the elimination of fluctuations (with respect to the classical ergodic measure) of generic chaotic eigenfunctions in the SL.

We are going to show that the highest peaks of the Husimi are not completely random. On the contrary, they live (within an error of order $`\sqrt{\mathrm{}}`$ in each direction) on short POs, and their phases (considering in this case the corresponding Bargmann function ) are connected semiclassically. That is, high fluctuations correspond with the scars of short POs introduced by Heller  (although the meaning of short is other; see Eq. (4)). In this respect, the following commonly used expression is definitely wrong: ”…scarred eigenfunctions do not contradict Shnirelman theorem because they are of null measure.”. Actually, scar phenomena are of a semiclassical nature (they manifest on scales going to zero with $`\mathrm{}`$). Then again, a classical smoothing washes out these very rich structures in the SL.

Bogomolny  and Berry  incorporated the idea of scars into the description of chaotic eigenfunctions. These authors concluded that chaotic eigenfunctions consist of a dominant universal contribution, decorated by scars of short POs; but this result was obtained after making an average over an energy interval that also eliminates fluctuations. In our study, we consider individual eigenfunctions, without either a classical smoothing nor an energy averaging. Then, we arrive to the conclusion that a generic chaotic eigenfunction is scarred by a set of short POs which characterizes the state, and the sum of their intensities survives the SL. This fact stresses that localization on short POs is the signature of chaos in this limit. Moreover, we will give the mean value and dispersion of the scarred intensities showing that relative fluctuations go to zero as $`1/|\mathrm{ln}\mathrm{}|`$, and providing a detailed description of a generic chaotic eigenfunction.

We have developed a systematic semiclassical construction of wave functions living in the neighbourhood of unstable POs. Reference  provides the construction of resonances without transverse excitations and ref.  applies the recipe to the Bunimovich stadium billiard. In ref.  we construct resonances of a PO with transverse excitations, all at the same Bohr-Sommerfeld (BS) quantized energy. Then, scar functions are the result of a minimization of the energy dispersion in this basis. Finally, ref.  gives the details for the construction of scar functions in the stadium billiard. Scar functions are the objects on which we are going to focus our description. We are thinking about structures living along the manifolds up to the first homoclinic point. Then, it is not necessary to include interference effects, and the description remains essentially simple.

Let $`\gamma `$ be an unstable $`\mathrm{𝐏𝐎}`$ of the system with period $`T_\gamma `$ and Lyapunov exponent $`\lambda _\gamma `$ per unit time (consider for the moment a conservative Hamiltonian system with two degrees of freedom). Let $`\varphi _\gamma `$ be the corresponding scar function with BS energy $`E_\gamma `$. Finally, let $`\phi _\mu `$ be the set of normalized eigenfunctions of the system with eigenenergies $`E_\mu `$. In ref.  we have shown that the set of intensities $`I_\mu =|\varphi _\gamma |\phi _\mu |^2`$ semiclassically satisfies, $`_\mu I_\mu =1`$, $`_\mu E_\mu I_\mu =E_\gamma `$, and

$$\sigma _\gamma \sqrt{\underset{\mu }{}(E_\mu E_\gamma )^2I_\mu }=\mathrm{}\lambda _\gamma \mathrm{\Gamma }/2\text{[12]}.$$
(1)
In the SL, the universal dispersion $`\mathrm{\Gamma }`$ goes to zero as $`2\pi /|\mathrm{ln}\mathrm{}|`$ (expresions of $`\mathrm{\Gamma }`$ for finite values of $`\mathrm{}`$ are given in ref. ). Then, the life time $`\tau _\gamma \mathrm{}/\sigma _\gamma `$ of $`\varphi _\gamma `$ diverges logarithmically in the same way as the Ehrenfest time does (but they are different times). This extremely low decay  is governed by a Gaussian law ($`|\varphi _\gamma (0)|\varphi _\gamma (t)|^2=e^{t^2/\tau _\gamma ^2}`$), and the smooth part of the intensities $`I_\mu `$ (the strength function) results a Gaussian function. Then, defining $`ϵ(EE_\gamma )\rho _E`$ (with $`\rho _E`$ the energy density), the mean value of the intensities at $`ϵ`$ is

$$\overline{I}(ϵ)=e^{ϵ^2/2n^2}/\sqrt{2\pi n^2},$$
(2)
where $`n\sigma _\gamma \rho _E`$ is the mean number of eigenenergies contained in one energy dispersion.

Another consequence of the Gaussian decay is that Eq. (1) is also valid for systems with $`f`$ degrees of freedom if we replace $`\lambda _\gamma `$ by $`(\lambda _i^2)_\gamma ^{1/2}`$, with $`\{\lambda _i\}_\gamma `$ the set of $`f1`$ positive Lyapunov exponents of $`\gamma `$ (for an exponential decay the change would be $`\lambda _\gamma `$ by the K-S entropy $`(\lambda _i)_\gamma `$).

Now, we are able to establish a criterium in order to decide when a $`𝐏.𝐎`$ is short. We propose a relation of the form $`T_\gamma <\beta \tau _\gamma `$, with $`\beta `$ a constant to be determined from the theory of short POs. This theory says that in the SL an eigenfunction is defined by POs of period lower than $`T_0=h_T^1\mathrm{ln}(T_Hh_T)`$ , with $`h_T`$ the topological entropy. Then, in the SL $`T_0`$ would be equal to $`\beta `$ times $`\tau `$ (the life time of a generic scar function), with

$$\tau =\pi ^1(\lambda _i^2)_{sys}^{1/2}|\mathrm{ln}\mathrm{}|.$$
(3)
$`\{\lambda _i\}_{sys}`$ is the set of positive Lyapunov exponents of the system. Finally, we say that a PO $`\gamma `$ is short when

$$T_\gamma /\tau _\gamma <\sqrt{2\pi }\eta \mathrm{"}short\mathrm{𝐏𝐎}condition\mathrm{"},$$
(4)
where the classical invariant

$$\eta \sqrt{\pi /2}(f1)h_T^1(\lambda _i^2)_{sys}^{1/2},$$
(5)
will play a central role below. Evidently, any PO is short at sufficiently high energies.

Of all the intensities of a given scar function we are actually interested in those with the highest values (they are certainly related to the phenomenon of localization on short POs). We are going to study the properties of these high intensities in terms of a statistical model. The main purpose is to estimate their mean values and dispersions, and to provide a range in the spectrum where they live.

The fluctuations of an intensity $`I`$ at $`ϵ`$, around its mean value $`\overline{I}(ϵ)`$ (see Eq. (2)), will be described as usual  by a chi-squared distribution with one degree of freedom for systems with time reversal symmetry (appendix D treats systems without time reversal). Then, the probability of finding an intensity lower than $`I`$ (the accumulated probability) is given, for $`aI/\overline{I}(ϵ)>1`$, by $`F_ϵ(I)=1\sqrt{2/\pi a}e^{a/2}(1a^1+3a^2+\mathrm{})`$. In order to obtain the distribution for an arbitrary intensity (independent of its position in the spectrum), we make an averaging; if $`N`$ is the number of intensities contributing to the scar function, $`F(I)(1/N)_{N/2}^{N/2}F_ϵ(I)𝑑ϵ`$. This expresion is computed in terms of Gaussian integrals after the following expansion: $`e^{a/2}=e^ye^{yϵ^2/2n^2}(1yϵ^4/8n^4+\mathrm{})`$, where $`y\sqrt{\pi /2}nI`$ (note that $`ϵ^2/n^2=𝒪(y^1)`$). Then,

$$F(I)=1\frac{\sqrt{2}n}{N}\frac{e^y}{y}\left(1\frac{9}{8y}+\frac{305}{128y^2}+\mathrm{}\right).$$
(6)
Now, if $`x_1`$ is the greatest of the intensities, $`x_2`$ the second one and so on, the probability density of $`x_j`$ is given by Eq. (16). A simple estimation of $`x_j`$ is derived from $`1F(x_j)=j/N`$. With this in mind, we define the random variable $`z_j`$ by the relation

$$1F(x_j)=e^{z_j}(j/N).$$
(7)
Combining Eqs. (6) and (7), $`x_j`$ is given by

$$y_j\sqrt{\pi /2}nx_j\alpha \mathrm{ln}(\alpha +9/8)+b+b^2/2,$$
(8)
with $`\alpha z_j+\mathrm{ln}(\sqrt{2}n/j)`$ and $`b\mathrm{ln}(\alpha +287/128)/(\alpha +17/8)`$. Equation (8) works for $`\alpha >1`$.

Equation (17) gives the mean value $`\overline{z_j}`$ and dispersion $`\sigma _{z_j}`$ of $`z_j`$. With these, the mean value of $`y_j`$ is obtained from Eq. (8) by setting $`\overline{\alpha }`$ ($`=\overline{z_j}+\mathrm{ln}(\sqrt{2}n/j)`$) in place of $`\alpha `$ and adding to the rhs the term $`c\sigma _{z_j}^2/2(\overline{\alpha }+9/8)^2`$ (because $`\overline{\mathrm{ln}(\alpha +9/8)}\mathrm{ln}(\overline{\alpha }+9/8)c`$). On the other hand, the dispersions of $`y_j`$ and $`z_j`$ are equal to the leading order. Then, we arrive to the first conclusion: relative fluctuations go to zero in the SL as follows,

$$\sigma _{x_j}/x_j1/[\sqrt{j}\mathrm{ln}(\sqrt{2}n/j)]=𝒪(1/|\mathrm{ln}\mathrm{}|).$$
(9)
The next question is to know where $`x_j`$ can be found. The probability $`p(ϵ)`$ of finding $`x_j`$ near $`ϵ`$ is proportional to the probability density $`dF_ϵ/dI`$ at $`x_j`$; that is, $`p(ϵ)e^{(y_j+1/2)ϵ^2/2n^2}(1y_jϵ^4/n^4+\mathrm{})`$. Then, $`x_j`$ is restricted to a range $`\mathrm{\Delta }E_j\sigma _\gamma /\sqrt{y_j+1/2}`$ around $`E_\gamma `$ which in units of $`2\pi \mathrm{}/T_\gamma `$ (the distance between consecutive BS quantized energies) goes to zero as $`|\mathrm{ln}\mathrm{}|^{3/2}`$. In conclusion: scars of $`\gamma `$ accumulate in the vicinity of $`E_\gamma `$ in the SL.

In the following, we are going to discuse a fundamental question. How many universal waves participate in the description of a chaotic eigenfunction?. For maps the question is simple; there are $`N_u`$ waves with $`N_u`$ equal to the area of the map divided by $`2\pi \mathrm{}`$. However, for Hamiltonian systems it is not clear in general which one is the right, if any, Poincaré surface of section.

At classical level, we can say that during an evolution for a time lower than the period $`T_{min}`$ of the shortest PO, it is impossible to extract non-universal information of the system. Using this time, Berry  distinguishes universal aspects of the spectral rigidity from the non-universal ones. However, $`T_{min}`$ is not very useful because it is simply an estimation. Evidently, we need a time representing a short universal evolution but such that its determination is based on a deep knowledge of the system. For these reasons we propose the inverse of the topological entropy for this time. Then, if an eigenfunction is defined after an evolution equal to $`T_H`$ and we can do the travel in steps no greater than $`h_T^1`$, where each step defines a universal wave, the number of required universal waves is

$$N_u=T_Hh_T.$$
(10)
There is a way of verifying the accuracy of Eq. (10). In billiards, Birkhoff coordinates define the right Poincaré surface of section (because all classical or quantal information is contained on the boundary). The area of the section is $`2\mathrm{}k`$, with $``$ the lenght of the boundary and $`k`$ the eigenwave number. Then, $`N_u=k/\pi `$, and assuming that also Eq. (10) is right, we arrive to the following expresion for the topological entropy $`h_T^{}`$ per bounce

$$h_T^{}=\overline{l}/\pi 𝒜,$$
(11)
with $`𝒜`$ the area of the billiard and $`\overline{l}`$ the mean length per bounce. We have verified Eq. (11) within an error of $`2\%`$ in the stadium billiard with radius unity and several different areas.

So far, we have described scar functions in the basis of eigenfunctions. Now according to the theory of short POs , eigenfunctions would be described equivalently in terms of scar functions. That is, imagine the spectrum of scar functions given by all BS energies of all short POs. This means for a given PO, all BS energies in those energy regions where $`\gamma `$ is short (following the criterium of Eq. (4)). Then, an eigenfunction $`\phi _\mu `$ with eigenenergy $`E_\mu `$ is represented in the basis of scar functions by the intensities $`I_\gamma =|\phi _\mu |\varphi _\gamma |^2`$, which are concentrated around $`E_\mu `$ in the spectrum of BS energies. The dispersion $`\sigma _\mu ((E_\gamma E_\mu )^2I_\gamma )^{1/2}`$ depends on the Lyapunov exponents of the orbits with BS energies in the vicinity of $`E_\mu `$. However, in the SL there is a uniformization and the energy dispersion results independent of the position in the spectrum and equal to $`\mathrm{}/\tau `$, with $`\tau `$ given by Eq. (3). The smooth part of the intensities $`I_\gamma `$ is given by Eq. (2), with

$$n=\sigma \rho _E=T_H/2\pi \tau ,$$
(12)
and $`ϵ=(EE_\mu )\rho _E`$ . If $`x_1`$ is the highest of the intensities, $`x_2`$ the second one and so on, the mean value of $`x_j`$ is given by Eq. (8) (and discusion thereafter), and the relative dispersion by Eq. (9). Of all short POs, those with possibility of having intensity $`x_j`$ satisfy $`|E_\mu E_\gamma |<\sigma /\sqrt{y_j+1/2}`$.

The following question is to decide a criterium for scarring. We will say that $`\phi _\mu `$ is scarred by $`\gamma `$ if $`I_\gamma `$ is greater than the greatest of the intensities provided by a random model of universal waves. The highest universal intensity is, to the leading order, $`I^{(u)}=2\mathrm{ln}(N_u)/N_u`$ . Then, using Eqs. (5), (10) and (12), the condition $`x_j>I^{(u)}`$ reduces to

$$y_j>\eta \mathrm{"}scarringcondition\mathrm{"}.$$
(13)
The number $`n_{scar}`$ of POs satisfying Eq. (13) is given by

$$\frac{n_{scar}}{n}e^\eta \sqrt{\frac{8\mathrm{ln}(1+1/\delta \eta )}{(1+4\eta /\delta )}},$$
(14)
with $`\delta =(9+\sqrt{73})/2`$. This formula (and the next too) interpolates the behaviors for large $`\eta `$ (obtained from Eqs. (6) and (7)) and $`\eta `$ going to zero (see appendix C). Finally, the sum of scarred intensities results

$$\underset{j=1}{\overset{n_{scar}}{}}x_j\frac{2}{\sqrt{\pi }}e^\eta \left(\frac{2}{\sqrt{\pi }}1\right)e^{2\eta }.$$
(15)
We emphasize that the unexpected results of Eqs. (14) and (15) depend decisevely on the use of scar functions in order to measure localization on short POs. On the contrary, by using wave packets in the transverse direction to the motion (the so called vacuum states in ref. ), there result $`n_{scar}/n=0`$ and $`x_j=0`$ in the SL. In this respect, the departure of universal behavior found in ref.  has a weight going to zero in the SL.

We stress that all formulae of statistical nature verify an impresive agreement with numerical simulations. Moreover, Eqs. (1) and (2) were extensively verified in the stadium billiard. Finally, we present an example where all the ideas developed in the article have been tested. Figure 1) shows the decomposition of a very high excited chaotic state (plotted in configuration space in ref. ). We found $`9`$ scarred short POs ($`n_{scar}8.5`$ from Eq. (14)) providing a contribution of the $`38\%`$ (Eq. (15) predicts $`33\%`$). Moreover, scarred intensities are in good agreement with predictions (see Fig. 1d)).

In conclusion, localization on short POs survives the SL and depends exclusively on $`\eta `$. Evidently, this localization will be strong in systems with few degrees of freedom. In particular for $`f=2`$, and assuming that $`\lambda h_T`$, this localization results a universal property.

This work was partially supported by SETCYP-ECOS A98E03.

Appendix A. Let $`I_1,I_2,\mathrm{},I_N`$ be a set of independent random variables, with common probability density $`f(I)`$, living in the range $`[0,\mathrm{})`$. Let $`x_1`$ be the greatest of the variables, $`x_2`$ the second one an so on. The joint probability density is given by $`p(x_1,x_2,\mathrm{},x_N)=N!f(x_1)f(x_2)\mathrm{}f(x_N)`$, for $`x_1x_2\mathrm{}x_N`$; otherwise $`p=0`$. Then, the probability density of $`x_j`$ results (after integration over the others variables)

$$p(x_j)=\frac{N!f(x_j)F(x_j)^{Nj}[1F(x_j)]^{j1}}{(Nj)!(j1)!},$$
(16)
with $`F(x)=_0^xf(x^{})𝑑x^{}`$ the accumulated probability.

Appendix B. We will derive the mean value and dispersion of the random variable $`z_j`$ defined in Eq. (7). Assuming that $`z_j=𝒪(1)`$, there results from Eq. (7) $`1F(x_j)=𝒪(j/N)`$, and then $`F(x_j)^{Nj}=\mathrm{exp}[(Nj)\mathrm{ln}F(x_j)]=\mathrm{exp}[N(1F(x_j))+𝒪(j^2/N)]`$. Moreover, $`N!/(Nj)!=N^j[1+𝒪(j^2/N)]`$. Using these approximations in Eq. (16), the probability density of $`z_j`$ is given by

$`p(z_j)j^je^{j(z_j+e^{z_j})}/(j1)!.`$

For instance, a numerical computation gives $`\overline{z_1}0.577`$ and $`\sigma _{z_1}1.28`$. On the other hand, for large values of $`j`$ we change to a new variable $`w=\sqrt{j}z_j`$, and expanding in powers of $`1/\sqrt{j}`$ there results

$`p(w){\displaystyle \frac{e^{w^2/2}}{\sqrt{2\pi }}}\left(1{\displaystyle \frac{1}{12j}}\right)\left(1+{\displaystyle \frac{a}{\sqrt{j}}}+{\displaystyle \frac{b}{j}}+{\displaystyle \frac{c}{j^{3/2}}}\right),`$

with $`a=w^3/6`$, $`b=w^4(w^23)/72`$ and $`c=w^5(5w^445w^2+54)/6480`$. Finally, using this density we have

$$\overline{z_j}\frac{1}{2j}+\frac{1}{12j^2}\mathrm{and}\sigma _{z_j}\frac{1}{\sqrt{j}}\left(1+\frac{1}{4j}\right),$$
(17)
in excellent agreement with numerical data.

Appendix C. For $`\eta `$ going to zero, $`n_{scar}`$ satisfies $`\overline{I}(n_{scar}/2)=I^{(u)}=\sqrt{2/\pi }\eta /n`$. Then, using Eq. (2) there results $`n_{scar}/n2^{3/2}\sqrt{\mathrm{ln}(1/\eta )}`$.

Appendix D. The corresponding formulae for systems without time reversal are the following: $`y_j\sqrt{2\pi }nx_j\alpha \mathrm{ln}(\alpha +3/4)/2+b+b^2`$, in place of Eq. (8), with $`\alpha z_j+\mathrm{ln}(\sqrt{2\pi }n/j)`$ and $`b=\mathrm{ln}(\alpha +31/16)/(4\alpha +5)`$. Equation (13) is the same. Equation (14) is replaced by

$`n_{scar}/ne^\eta \sqrt{2\pi \mathrm{ln}(1+1/\eta )},`$

and Eq. (15) by

$`{\displaystyle \underset{j=1}{\overset{n_{scar}}{}}}x_j{\displaystyle \frac{(1+\eta )}{\sqrt{3/4+\eta }}}e^\eta \left({\displaystyle \frac{2}{\sqrt{3}}}1\right)e^{2\eta }.`$