Mode resolved travel time statistics for elastic rays in three-dimensional billiards

0.58 MB·English

by  Adrian Ortega

Mode resolved travel time statistics for elastic rays in three-dimensional billiards
A. Ortega,1 K. Stringlo,2 and T. Gorin2
1Instituto de F´ısica, Universidad Nacional Auto´noma de M´exico, M´exico D. F. 01000, M´exico
2Departamento de F´ısica, Universidad de Guadalajara, Guadalajara 44840, Jal´ısco, M´exico
(Dated: January 16, 2012)
We consider the ray limit of propagating ultrasound waves in three-dimensional bodies made
from an homogeneous, isotropic, elastic material. Using a Monte Carlo approach, we simulate the
propagation and proliferation of elastic rays using realistic angle dependent reflection coefficients,
takingintoaccountmodeconversionandray-splitting. Forafewsimplegeometries, weanalysethe
long time equilibrium distribution focussing on the energy ratio between compressional and shear
waves. Finally, we study the travel time statistics, i.e. the distribution of the amount of time a
2 giventrajectoryspendsasacompressionalwave,ascomparedtothetotaltraveltime. Theseresults
1 areintimatelyrelatedtorecentelastodynamicsexperimentsonCodawaveinterferometrybyLobkis
0 and Weaver[Phys. Rev. E 78, 066212 (2008)].
2
PACSnumbers: 05.45.Mt,43.35.+d,43.40.+s
n
Keywords: Codawaveinterferometry,elasticrays,raysplitting,three-dimensionalbilliard
a
J

3
I. INTRODUCTION Thismakesitverydifficulttosimulateray-dynamicsover
1
longtimes. Oneofourmainachievementsconsiststhere-

fore in the development of an efficient algorithm which
] Elastic rays are the fundamental building block of the
D applies a Monte Carlo approach to deal with these diffi-
theoryofCodawaveinterferometry[1]whichoverrecent
C years has been developed into a well established method culties. Thebodiesemployedaresimilartotheonesused
by Lobkis and Weaver but not identical. Still our sim-
. for the analysis of seismological data [2], among oth-
n ulations allow to verify certain assumptions made about
ers[3]. Tothebestofourknowledge,onlylittleisknown
i
l aboutelasticraysasthe particlelimitoftheelasticwave the wave field and to analyse the effects of eventual vio-
n
lations.
[ equation,eventhoughitispreciselythatlimitthetheory
in Ref. [1] relies upon. The main difficulty with elastic Together with the introduction, this paper is divided
2 into 6 sections. The followingSec. II defines ourconcept
rays is due to the phenomenon of “ray-splitting” which
v ofanelasticray. Sec.IIIdescribestherandommodecon-
causes the dynamics to become effectively statistical in
8
versionmodelintroducedinRef.[8]. Ourmainresultsare
nature [4].
7
describedinSec. IV,andtheir relationto experiments is
5 There is a close analogyin elastic raysandCoda wave
discussedinSec.V.WeconcludethepaperwithSec.VI.
2 interferometry on the one hand and the orbits of clas-
0. sically chaotic quantum systems and semiclassical the-
1 ory in the diagonal approximation [5–7] on the other.
1 This analogy is crucial also for measurable quantities II. ELASTIC RAYS
1 such as the distortion (introduced in Ref. [8]) and scat-
:
v tering fidelity [9–12]. In Ref. [8] Lobkis and Weaver Rays are a well known concept from geometric optics,
i started a series of experiments with reverberant ultra- where they are used to approximate the propagation of
X
sound in three-dimensional Aluminum samples. These electromagneticwaves in the limit where the wavelength
r experiments could be described in terms of the Coda is small compared to the typical dimension L of the sys-
a
wave interferometry (based on the ray picture of a dif- tem. A single ray stands for a transversally concen-
fuse wave field), but also in terms of scattering fidelity tratedelectromagnetic wavepacketwhich movesthrough
and random matrix theory [11, 13]. So far the focus an optical system. In the propagation direction, the
has been on the form of the fidelity decay [13], but not wavepacket may also be localized, but that need not be
on the overall decay time of the fidelity signal. In or- so. Without localization one arrives at a stationary sit-
der to explain their results Lobkis and Weaver assumed uation, where wave energy continuously flows along the
that after a certain transient time, the elastic wave field ray path.
settles on a equilibrium state, where the energy is dis- In the present case, we assume that the wavepacket
tributed homogeneouslyandisotropicallyoverthe whole is finite, even in the longitudinal direction. It therefore
Aluminum sample, with relative energy shares in shear- contains a finite amount of energy Q. Since we neglect
andcompressionalwavescorrespondingtothe equiparti- any form of energy dissipation, that amount of energy is
tion ratio [8, 14]. conserved for all times. However, due to mode conver-
In this paper, we present numerical simulations of the sion and ray-splitting, the amount of energy Q is shared
propagation of elastic rays in finite three-dimensional among an increasing number of branches of the elastic
bodies. Due to ray-splitting an elastic ray spreads into ray.
anexponentiallyincreasingnumberofdifferentbranches. Elastic wavesin a homogeneousandisotropic medium
2
come in two different modes: P-waves (compressional
1.5
waves) where the medium undergoes harmonic displace-
ments in the longitudinal (propagation) direction, and 1
S-waves, where the medium undergoes harmonic dis-
0.5
placements in the transversal direction (shear waves).
Both wave forms have different propagation speeds, c 0
d x
y
and cs respectively. For the Aluminum samples studied R -0.5
in [8, 11, 13]:
-1
c =635cm/ms, c =315cm/ms. (1) -1.5
d s
-2
Asweshallseebelow,bothwaveformsconvertaccording
0 0.2 0.4 0.6 0.8 1 1.2 1.4
tospecificrulesintoone-anotherwhentherayisreflected
at the free surfaces of the body. In the classical ray pic- ϑ
ture this is handled by providing the elastic ray with an
additionalinternaldegreeoffreedom,aswillbeexplained FIG. 1. The different reflection coefficients as a function of
in detail below, in Sec. IIB. the angle of incidence ϑ. RPP (red line), RSP (green line),
One should be well aware of the fact that the ray pic- RSS (blue line), RPS (pink line). Note, Ryx means incident
ture, when applied to propagating elastic waves in a fi- x-waveand reflected y-wave.
nite solid, breaks down very soon. The transversal size
of a localized wavepacket increases rapidly in time such
Following [22], when a plane elastic wave hits a free
that after a few reflections, it reaches the extention of
surface, we first define a reflection plane. This is the
the whole body. In quantum mechanics this time scale
plane spanned by the surface normal and the propaga-
would be called the Ehrenfest time [17]. However, this
tion direction of the ray. An incident P-wave splits into
doesnotmeanthattheraypicturebecomesuselesswhen
an out-going P-wave and an out-going S-wave with po-
longer times are involved. In the field of quantum chaos
larization vector lying in the plane of reflection. The
it is well established that one can construct semiclassi-
reflection amplitudes for both waves are given by:
cal approximations, on nothing else than these classical
trajectories[18]. Theseapproximationsmayremainvalid 2 2
sin(2ϑ) sin(2Θ)−κ cos (2Θ) sinϑ
formuchlongertimes–inthecaseofchaotictwodegree- RPP = sin(2ϑ) sin(2Θ)+κ2 cos2(2Θ) , κ= sinΘ
of-freedomsystemsuptotimesoftheorderoftheHeisen-
(2)
berg time [19]. Ultimately, the present work might help
to pave the way towards a similar semiclassical theory 2κ sin(2ϑ) cos(2Θ)
R = , (3)
for elastodynamic systems. SP sin(2ϑ) sin(2Θ)+κ2 cos2(2Θ)
Forthenumericalsimulations,weusegeometries(rect-
whereκ=c /c >1suchthatthe exitangleΘisalways
angularblock,tetrahedron)withlongstraightedges,be- d s
smaller than the entrance angle ϑ. These angles are al-
cause similar bodies have been used in the experiments
waysmeasuredwithrespecttothenormalofthesurface.
by Lobkis and Weaver and also because of their simplic-
The reflection amplitudes R and R determine the
ity. These bodies however have the disadvantage that PP SP
amplitudes of the two out-going waves.
diffraction on the edges and corners may have consider-
The case of an incident S-wave is more complicated.
able effects [20, 21]. The inclusion of diffractive orbits
Beforeconsideringthereflectionitself,wehavetodecom-
will be left to a future work.
pose the wave into one component with polarization in
the reflection plane (SV-wave) and another with polar-
ization perpendicular to it (SH-wave). For the SV-wave
A. Reflection coefficients
we then have similar reflection coefficients as for the P-
wave:
Ingeneral,raysareofahybridnature. Fromamacro-
sin(2ϑ) sin(2Θ′)−κ2 cos2(2ϑ) sinΘ′
scopic perspective, the ray has negligible width and is R = , κ=
regarded as a one-dimensional object, while, from a mi- SS sin(2ϑ) sin(2Θ′)+κ2 cos2(2ϑ) sinϑ
croscopic perspective, the ray is seen as a plane wave (4)
whichallowstouserelativelysimplerulestocalculateits −κ sin(4ϑ)
R = , (5)
behavior under reflections. Since the system is assumed PS sin(2ϑ) sin(2Θ′)+κ2 cos2(2ϑ)
to be macroscopicinsize,the shape ofthe reflectingsur-
face does not really matter. It is always approximated where now Θ′ > ϑ such that these equations only apply
locallyasaplanesurface. Forourpurposesitistherefore as long as κsinϑ < 1. This introduces the critical angle
sufficienttoconsiderthereflectionlawsofplanewavesat (of incidence) ϑcr =arcsin(κ−1). However, for SV-waves
plane surfaces. For simplicity, we shall also assume that incident at larger angles, the reflected P-wave becomes
all reflecting surfaces are free surfaces. a surface wave with negligible contribution to the wave
3
field, while the reflection amplitude for the SV-wave be- of the reflected ray is obtained solely from the SV com-
comes ponent,wherasthe S-wavebranchisgivenasa superpo-
sition of the reflected SH-wave and the S-wave branch
cos2(2ϑ)−2iβ sinϑ sin(2ϑ) of the SV-component. These relations explain why it is
R =− , (6)
SS cos2(2ϑ)+2iβ sinϑ sin(2ϑ) necessary to keep track of the polarization direction of
the S-mode rays.
whereβ = sin2ϑ−κ−2. Theabsolutevaluesquaredof For the classical description of the secondary rays, we
R is onepin that case, whichmeans that allthe energy need neither the wavepacket-widths nor the local inten-
SS
is transferred to the reflecting S-wave. sities,however,wedoneedthetotalenergyshareofeach
SH-waves, i.e. shear waves with polarization direction branch. In principle, the calculation of the total en-
perpendicular to the reflection plane, cannot convert to ergyshareisratherinvolvedsinceitrequirestheintegra-
P-waves. They are reflected according to the standard tion of the energy flux across a surface perpendicular to
law where the reflection angle is equal to the angle of the propagation direction, taking into account the vari-
incidence. ation of the wavepacket intensity. Fortunately, for the
reflected branch without mode conversion, neither the
wavepacket’s intensity profile nor its propagation speed
change. As a consecuence, the total energy contained in
B. Classical ray limit
that branch is simply proportional to the intensity cal-
culated from the corresponding reflection coefficients:
In order to define the classical ray limit, we consider
a localized wavepacket of total energy Q. At any mo- • An incident P-wave carrying the energy Q will
ment in time this wavepacket has a well defined posi- transfertheenergyR2 Qtothe reflectedP-mode
PP
tion ~r and propagation direction ~e . With respect to branch.
v
its momentary position ~r, the wavepacket has a certain
• An incident S-wave of the same energy whose po-
extension along the propagation direction and perpen-
larization direction makes an angle θ with the re-
diculartoit. Typically,wewouldimagineacigar-shaped
flectionplanemustbesplitintotheSV-component
wavepacket,orientedalongthe propagationdirection. In
of energy Q = Q cos2θ and the SH-component
the classical limit, we ignore the transversal extension V 2
of energy Q = Q sin θ. The total energy of the
of the wavepacketand obtain thereby a one dimensional H
reflected S-mode branch then is R2 Q +Q .
object. For the studies to follow, the wavepacket’s ex- SS V H
tension along the propagation direction does not mat- Energy conservation now implies that the respective
ter. Sinceelasticwavespropagateintwodifferentmodes, mode-convertedbranchmustcarrytheremainingenergy.
theclassicaldescriptionmustbeextendedbysomeinter- That this is indeed the case, is shownexemplarily for an
nal variables: We need one binary variable to specify incident P-mode ray in the appendix.
the mode of the wavepacket which may be longitudinal Returningtothepropagationoftheelasticray,wenote
(P) or transversal (S). Moreover, if the wavepacket is that subsequent reflections lead to further subdivisions
in transversal mode, we need to record the polarization into an exponentially increasing number of branches.
direction by a unit vector~e , perpendicular to the prop- Topologically,anelasticraymaybeconsideredasatree,
p
agation direction. wheretheinitialenergyistransportedfromthetrunkto-
The classicaldescriptionof elastic wavesruns into dif- wardstheouter(higherorder)branches. Evenformodest
ficulties whenitcomestoreflections. The problemisray traveltimes, it is soonimpossible to keeptrackofallthe
splitting [4], whichwillbetreatedasfollows: We assume branches of the elastic ray.
that in the vicinity of the wavepacket center, an elastic We therefore employ a Monte Carlo method to sam-
ray may be described as a plane wave. For that plane ple implicitely only over those branches which carry the
wave, the reflection laws of the previous section IIA ap- largest amount of energy. For that purpose we trans-
ply. Asexplainedthere,thereflectedwavewillingeneral late energy share into probability. Thus, we replace the
besplitintoanS-modeandaP-modebranch,propagat- deterministic energy distribution model by a statistical
ing into different directions. Accordingly, the reflected model, where we define an ensemble of rays which all
ray splits in two, where the local intensities and polar- start in the same initial state corresponding to a wave
izationdirectionsinthewavepacketcentermaybecalcu- packet with energy Q0. At each reflection, we choose at
latedfromthereflectioncoefficientsofthecorresponding random,whethertherayfollowsonebranchortheother.
plane waves. Finally,purely geometricconsiderationsal- Theprobabilitywithwhichoneoptionortheotherischo-
low to determine the widths of those branches. In that sen, are simply given by the relative energy shares cal-
way, we obtain the complete information about the two culatedfromthecorrespondingreflectioncoefficients. At
reflected branches of the incident ray. Note that an inci- theend,theprobabilitytotravelthroughacertainhigher
dentS-moderaymustbeconsideredasalinearcombina- order branch is given by the product of probabilities for
tionofaSV-wave(polarizationdirectioninthereflection the choicesmadealongthe historyofthe givenray. This
plane)andaSH-wave(polarizationdirectionperpendic- probability agrees with the energy ratio between the to-
ular to the reflection plane). Hence, the P-wave branch tal energy of the wavepacket in that particular branch
4
and the energy of the initial wavepacket. Thus, in a nu- The authors obtain this from an analogy to room acous-
merical simulation, the energy share of a higher order tics[15]andadetailedcalculationofthemodeconversion
branchcanbeestimatedfromthenumberofmembersin probability during individual reflections [16]. The result
the ensemble which terminate in that given branch. is:
Our ensemble of rays may as well start in different
S
initial states also chosen at random. That is the case, β =0.59c . (9)
d
4V
when we intend to calculate the evolution of a statisti-
cal ensemble of initial conditions. In what follows, we In order to explain the experiments performed on the
will consider the following three different types of initial various Aluminum samples, one requires the statistics of
conditions: tP, the accumulated amount of time a given ray of du-
rationt spends in P-mode. As an estimate for the mean
(i) Deterministic initial conditions, where we com- ht i and its variance var(t ) the authors find
P P
pletely specify one particularray,starting atacer-
tain point, in a specific direction, either in S- or t 2t R2
ht i= , var(t )= . (10)
P-mode, and in case with a specific polarization. P R+1 P β (1+R)3
(ii) Surface initial conditions, where we start P-waves The first equation is easily explained on the basis of er-
at a particular point on the surface of the body, godicity. A more involved calculation is required for the
while the direction of the ray is chosen at random second.
within the half sphere pointing into the body. For Fig. 2 we performed numerical simulations for the
randommodeconversionmodeldescribedabove. Weuse
(iii) Homogeneousinitialconditions,wherewestartrays random numbers with exponential probability densities
at random positions inside the body, with random togeneratesequencesofalternatingP-andS-modetime
directions (isotropic on the whole unit sphere) in segments forming a random realisation of an elastic ray.
S- or P-mode with probabilities chosen according With the total time t being fixed we can than measure
to the equipartition ratio in Eq. (7) and in case the average P-mode occupation time ht i as well as its
P
random polarization direction. variance. These quantities areshownin Fig. 2 as a func-
tion of t, for three different cases: When the trajectory
starts with a P-mode segment, with a S-mode segment,
III. RANDOM MODE CONVERSION MODEL or by choosing P-mode or S-mode at random acording
to the equipartition ratio. We find that in all cases, ht i
P
In Ref. [8] the authors introduce a simple model, ig- and var(tP) quickly convergeto the theoreticalvalues as
noring any geometric effects, where elastic rays undergo t becomes sufficiently large.
random and statistically independent mode conversions
according to the conversion rates α (S-to-P) and β (P-
IV. NUMERICAL RESULTS
to-S). Here,themeanfreetime betweentwoconversions
is given by α−1 (for the S-mode segments) and β−1 (for
theP-modesegments),respectively. Asanadditionalin- In our simulations we use bodies of two different
gredient, the authors invoke the equipartition principle, shapes,arectangularblockofdimensions9cm×13cm×
which states that energy should be distributed equally 7.6cm (giving rise to integrable or possibly pseudo-
among the different degrees of freedom. Taking into ac- integrable dynamics) and a regular tetrahedron with
countthedifferentwavespeedsforP-andS-modewaves, edges of length 10cm, where the dynamics is probably
the ratio of the energy share between S- and P-mode ergodic. In addition, we study each of the two bod-
waves is [14] ies with and without an internal sphere. That sphere
is suposed to provide another free surface for the waves
c 3 (rays)movinginsidethebody,whichrendersthedynam-
d
R=2 . (7) ics chaotic. The inner sphere for the rectangle is chosen
(cid:18)c (cid:19)
s
to be relatively large in order to reduce bouncing ball
This number also determines the ratio between the con- orbits. The bodies (including the inner spheres) are de-
version rates, since the equilibrium condition of the cor- picted in Fig. 3. According to the random mode conver-
responding rate equation demands that on average, the sionmodel,thevolumeandthesurfaceareaofthebodies
numberofelasticraysinP-modeandS-modearerelated are important parameters. These are given in table I.
by As explained earlier, the simulation of the classical
rays is done by launching a large number of rays (up
hN i β to 4×106) with different initial conditions. Those are
S
= =R. (8)
hN i α chosen to be either of type (ii) (Surface) or of type (iii)
P
(Homogeneous). The formermight be consideredas cor-
One last piece of information is necessary in order to responding to the experimental situation in Ref. [8, 13],
be able to estimate the values of the conversion rates. but if at all this may be true only qualitatively.
5
Body Volume [cm3] Surface [cm2]
0.11 Rectangular block 889.2 568.4
...with innersphere 709.6 722.3
Regular tetrahedron 117.8 173.2
0.09 ...with innersphere 110.6 191.3
t
/
i
P
ht 0.07 TABLE I.
0.05
1.2
t 0.0008 1.15
/
) 0.0006
P
(t 0.0004
ar 1.1
v 0.0002 β
/
m
0 u
n
0 0.5 1 1.5 2 2.5 β 1.05
t/[ms]
1
FIG. 2. Mean and Variance of the P-mode travel time for
theLobkis’andWeavers’srandomconversionmodel. TheP- 0.95
to-S-conversion rate was chosen as β = 162. The different 0 0.5 1 1.5 2 2.5 3 3.5 4
colors refer todifferent initial conditions: starting in S-mode t/[ms]
(red), starting in P-mode (blue), starting randomly in S- or
P-mode according to the equipartition ratio (green). In the
FIG. 4. The conversion rate βnum divided by the travel time
case of the equipartitioned initial condition, the solid green
t, determined from the mean free P-mode occupation time,
line shows the value of the theoretical expectation, Eq. (10).
as a function of t, for homogeneous initial conditions, for all
Othersolid lines (red and blue) show simple rational best fit
fourdifferentgeometries. Tetrahedronwithsphere(redline),
functions to guide theeye.
withoutsphere(greenline),Rectanglewithsphere(blueline),
and without sphere (pinkline).
ray happened to travel in P-mode. The average over
those periods over all rays is just the mean free P-mode
travel time, and therefore equal to β−1. In Fig. 4 the so
determinedconversionrateβnum is comparedto the the-
oretical estimate (9) for all four geometries considered.
Wefindthatassoonastislargeenough(fortherandom
mode conversion model to become valid, there need to
occur sufficiently many reflections), βnum is quite close
to the theoretical estimate. In fact, as can be observed
FIG. 3. On left a regular tetrahedron with length of 10 cm
in Fig. 4, βnum never deviates more than 3% from the
with inner sphere radius 1.2 cm on right a rectangular block
theoretical estimate. Nevertheless, there are systematic
withlengthof9cmwidthof13cmandheightof7.6cmwith
innersphere radius 3.5 cm. differences for the different geometries. While βnum ends
up about 2% above the theoretical estimate in the case
ofthe rectangularblock withinner sphere, βnum ends up
about 3% below, in the other three cases. We believe
A. Conversion rates and equipartition ratio
that this behavior is still acceptable in view of the fact
that Eq. (9) relies on rather rough estimates.
Toverifytheaccuracyoftherandomconversionmodel
outlined above, we start by analysing the conversion
rates. Forthatpurposeweperformsimulationswherethe
1. Equipartition ratio
initial conditions of the rays are chosen according to the
expected equilibrium state. Thus, we start rays at ran-
dompositionsinsidethebody,withrandomdirections,in IntheFigs.5and6,westudytheenergysharebetween
S- or P-mode according to the theoretical equipartition S- and P-mode rays. Since our simulations assume an
ratio (7), and if in S-mode we choose the polarization equal amount of energy associated with each ray, that
direction also at random. For each ray of pre-defined energysharecanbe computedfromthe relativefrequen-
duration t, we then record all periods during which the cies with which we find a particular realization of the
6
50 20
40
15
o o
ati 30 ati
r r
y y
nerg nerg 10
Pe 20 Pe
/ /
S S
5
10
0 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
t/[ms] t/[ms]
FIG. 5. S-mode vs. P-mode energy ratio for samples of dif- FIG. 6. S-mode vs. P-mode energy ratio for the rectangular
ferentgeometries. RaysarestartedonthesurfaceinP-mode blockwithandwithoutinternalsphere. Theinitialconditions
withrandomdirections(Surfacetypeinitialconditions). The were chosen homogeneous and isotropic (blue and pink line)
color coding for the different geometries is the same as in andstartingfromafixedpointonthesurface(alwaysstarting
Fig. 4. The black horizontal line shows the equipartition ra- in P-mode) (red and green line). The black horizontal line
tio R, according toEq. (7). shows the valueof theequipartition ratio R, Eq. (7).
inner sphere, this is not the case, and the equilibrium
rays in either one of the two possible modes. In Figs. 5
values a very different.
and 6, we plot the ratio of these frequencies vs. t, the
Extendingthesesimulationsuptotimesoftheorderof
travel time.
5ms, we have confirmed that the S/P-energy ratio con-
InFig.5thesimulationalwaysstartswithinitialcondi-
verges to finite values in the limit of large times, except
tionsonthesurface[initialconditionsoftype(ii)],where
5 for the rectangular block without inner sphere. Surpris-
we applied the Monte Carlo sampling with 4×10 ran-
ingly, we have found that these values do not depend on
domrealizations,foreachofthefourdifferentgeometries.
the type of initial conditions applied (Surface or Homo-
Theoretically, we expect that the S/P energy ratio con-
geneous). Concretely, we found the following values:
verges to R = as given in Eq. (7). Since in the present
case,theraysalwaysstartinP-mode,theenergyratioat
16.538 : Tetrahedron without sphere
small times is close to zero. We observe that two tetra-
16.308 : Tetrahedron with sphere
hedrons show a similar behavior, which is very distinct R= , (11)
 14.685 : Rectangle with sphere
from that of the two rectangular blocks. They approach
9.5/3.2 : Rectangle without sphere
the theoretical value for R via damped irregular oscilla- 
tions which are initially very strong. Comparingthe two while the theoretical value is R = 16.384. We can see
tetrahedrons, the presence of the inner sphere, tends to that the rectangular block without inner sphere must be
reduce the strength of the oscillations such that the the- considered separately. Only in that case does the equi-
oretical value for R is reached faster. By contrast, the librium distribution of the energy share depend on the
rectangularblocksshowalmostnooscillationsatall,and initial state (for homogeneous initial conditions we get
quickly go over to a smooth approach of different limit R=9.5, otherwise R=3.2).
values for the energy ratio. In doing so, the rectangle
with inner sphere comes much closer to the theoretical
value for R than the rectangle without inner sphere. B. P-mode occupation time statistics
Fig.6showssimulationsforthetworectangularblocks
with initial conditions of type (ii) (on the surface) and In this section, we restrict ourselves to homogeneous
type (iii) (homogeneous). At small times, the energy initialconditions. Inoursimulations,everytrajectoryat
ratio shown starts at zero for initial conditions on the each reflection makes a random choice whether to do a
surface, and at the value of the equipartition ratio R, mode conversion or not. This results in different trajec-
for homogeneous initial conditions. As t increases, the torypathsanddifferentamountsoftimes,therayspends
curves first show some minor fluctuations and then start ineachmode. Inwhatfollows,westudy the distribution
to converge to different equilibrium values – in all cases of the P-mode occupation time t ; obviously, the corre-
P
clearlybelowtheequipartitionratio. Fortherectangular spondingS-modeoccupationtimewouldprovideexactly
block with inner sphere, we observe that independent of the same information. We have chosen this particular
the type of initial conditions, the energy ratio converges quantity because it may be linked to experimental re-
tothesameequilibriumvalue. Fortherectanglewithout sultsin[8,13]withthetheoreticalmodelsproposedin[1]
7
0.1 0.3
0.095 0.27
0.24
0.09
0.21
0.085
/t 0.18
i/t 0.08 )βP 0.15
htP 0.075 var(t 0.12
0.07 0.09
0.065 0.06
0.06 0.03
0.055 0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
t/[ms]
t/[ms]
FIG. 8. For homogeneous initial conditions, the variance of
FIG.7. Forhomogeneousinitialconditions,themeanP-mode
the P-mode occupation times for the different geometries,
occupation time htPi for the different geometries, divided by
rescaledbyβ/t. Thecolorcodingforthedifferentgeometries
thetotaltraveltimet. Wecomparewiththetheoreticalesti-
isthesame asin Fig. 7. Wecompare with thetheoretical es-
mate, Eq. (10), where we have inserted the numerical values
timate, Eq. (10), where β has dropped out (black horizontal
giveninEq.(11). Thecolorcodingisasbefore: Tetrahedron
line). Correcting the theoretical expectation for the β-values
with inner sphere (red points), without inner sphere (green
observed in Fig. 4, leads to the blue horizontal line for the
points), Rectangular block with inner sphere (blue points),
rectangular block with innersphere.
without inner sphere (pink points). The corresponding es-
timates for the asymptotic value are depicted as horizontal
lines of thesame color.
(Coda wave interferometry) and [8] (random mode con-
version model).
We start with the average P-mode occupation time block with inner sphere. In the case of the rectangular
ht i,forwhichtherandommodeconversionmodelmakes block without sphere, we find a steep linear increase of
P
the prediction (10). In Fig. 2(a) we have seen that this var(t )β/t, such that a comparison with a theoretical
P
predictionisindeedveryaccurate,providedtheensemble limit value makes no sense.
ofraysis in the standardequilibriumstate. Fig. 7 shows
ht i/t for homogeneous initial conditions which accord-
P
ingtotheoryshouldconvergeto(1+R)−1 atsufficiently In this last figure we find the largest deviations from
large times. Indeed for the tetrahedron as well as the the random mode conversion model. The result for the
rectangle with inner sphere, the prediction is fulfilled. tetrahedronwithsphereisabout25%belowthetheoreti-
The ratio ht i/t converges to the predicted value, if we calvalue;withoutsphere,itisabout5%below,whilethe
P
replace the theoretical equipartition ratio with the one result for the rectangle with sphere is more than 100%
obtained numerically from our simulations (these values above the theoretical value. The result for the rectangle
are given in Eq. (11)). In the case of the rectangle with- without sphere is totally off the scale.
outinnersphere(outofrange),theconvergenceofht i/t
P
ismuchslowerbutthelimitvaluestillconsistentwiththe
numerical S/P energy partitioning. Finally, we show in Fig. 9 the distribution of P-mode
Fig. 8 shows the variance of the P-mode occupation occupationtimest fortrajectoriesofdurationt=5ms.
P
times, divided by t/β with β calculated from Eq. (9). Inthecaseofthetetrahedrons,wefindthatthedistribu-
The black solid horizontal line shows the theoretically tionsareclosetoGaussians,withthedistributionforthe
expectedvalueaccordingtotherandommodeconversion tetrahedron with inner sphere being a bit narrowerthan
model (10), where β drops out (black horizontal line): for the tetrahedron without. In the case of the rectan-
gular blocks, the distributions are clearly non-Gaussian.
β 2R2 For the rectangle with inner sphere the distribution is
var(t ) → ≈0.1022. (12)
P t (1+R)3 strongly asymmetric, while for the rectangle without in-
ner sphere, the distributions is surprisingly close to a
Correcting the theoretical expectation for the numeri- Lorentzian. This latter observation can explain the fact
callydeterminedβ-values showninFig. 4doesnotmake that the variance of the P-mode occupation times scales
2
a noticeable difference in the case of the Tetrahedrons with t rather than t, since the variance of a Lorentzian
and leads to the blue horizontal line for the rectangular distribution is infinite.
8
by β−1, which should then also give a unique value; see
Eq. (12). For the Aluminum blocks of different shapes,
analysed in Ref. [8], the result was a wide spreading of
experimental values where, the smallest values roughly
20 agreedwiththetheoreticalexpectation,whilethelargest
values where about four times larger [23]. On a qual-
15 itative basis, one could observe the tendency that less
30
1 chaotic geometries lead to larger values for the slope.
×
10
nts
u
o
c 5
B. Our results
0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 1.5 2 The overall result of the present analysis, shown in
tP/[ms] Fig. 8, is similar in this respect. We also find that the
slopes can be very different for different geometries. It
however shows that the relevant quantity is not really
FIG. 9. For homogeneous initial conditions, the histograms
chaos as measured by Lyapunov exponents but possi-
for the tP statistics for t =5ms. The result for the tetrahe-
bly rather ergodicity. From the four different geometries
dronson theleft, with/without innersphere(red/green solid
considered, the tetrahedron with inner sphere has the
line); for the rectangular blocks on the right with/without
innersphere (blue/pink solid lines). smallest slope. It is notable in that case, that the slope
is about 25% below the theoretical value, which shows
that the theory doesn’t provide a lower bound as one
V. DISCUSSION might have been conjectured from the experimental re-
sults. Next comes the tetrahedronwithout inner sphere.
A. Related experimental results Because all its surfaces are plane, the Lyapunov expo-
nentmustbe zeroinanycase. However,thetetrahedron
has most likely ergodic dynamics. And we find indeed
InRefs.[8,13],theauthorsmeasuredtheacousticlong
that the slope is quite close to the theoretical expecta-
timeresponseofshortinitialultrasoundpulsesappliedto
tion. Onlyataconsiderabledistance,wefindtherectan-
different Aluminum bodies. They studied the cross cor-
gle with inner sphere. The inner sphere clearly leads to
relationsbetweenthesesignalstakenatdifferenttemper-
chaoticdynamicswithpositiveLyapunovexponents,but
atures. Since the temperature change induces a change
there are also regionsof integrable dynamics and bounc-
ofthe propagationspeeds ofP-andS-wavesbydifferent
ing ball orbits. That is apparently sufficient to increase
amounts, a temperature change results in the distortion
the slope to more than twice the theoretical value. The
oftheacousticsignalandareductionofthecrosscorrela-
rectangle without inner sphere is clearly the most reg-
tions. The reduction of the cross correlations,which has
ular body. For that geometry, the var(t ) curve rather
alsobeenidentifiedasascatteringfidelity[9,10,12],can P
shows a quadratic dependence on t. This is in line with
be described quantitatively within the theory of Coda
the probabilitydensity for the P-mode occupationtimes
wave interferometry [1].
shown in Fig. 9. Its shape is almost a Lorentzian which
As this theory shows, scattering fidelity (or “distor-
would imply an infinite variance.
tion”howthatquantityhasbeencalledinRef.[8])ises-
Our simulations have shown that many of the as-
sentiallygivenby the distribution ofP-mode occupation
sumptionswithintherandommodeconversionmodelare
times – the quantity studied in the previous Sec. IVB.
rathernicely met. This is true for the P-to-S conversion
AssumingGaussianstatisticsforthesetimes,thescatter-
rate β (Fig. 4), andalso for the equipartitionratio if the
ing fidelity or distortion may be related to the variance
dynamics is ergodic(Fig. 5). We couldalso confirmthat
of the P-mode occupation times as follows:
the averageP-modeoccupationtime isalwaysrelatedto
D(t)=2π2 (∆T f)2 (δ −δ )2 var(t ), (13) the S/P-energy ratio (Fig. 7). Thus, the weak point of
p S p
the randommode conversionmodelis clearlythe in gen-
where f is the carrier frequency of the elastic wave, ∆T eralinaccurateorwrongpredictionofthevarianceofthe
is the change in temperature, while δ = −1.685 × P-mode occupation times. This shows that the succes-
P
10−4K−1, and δ = −2.9×10−4K−1 are the thermal sion of P-mode and S-mode segments is in general not
S
dilation coefficients for P-waves and S-waves, respec- well described by a Poissonian process. We believe that
tively [8]. therearetwoeffectscomingintoplay. (i)Thedurationof
In Ref. [8], Lobkis and Weaver measured the slope of the individual P-mode segments (these are the shortest
D(t) scaled by ∆T2, f2 and β−1, which according to ones) cannot be really exponentially distributed because
therandommodeconversionmodelshouldalwaysbethe the rays have to travel a certain distance before having
same. For our simulations, this is equivalent to deter- the possibility to undergo a mode conversion. (ii) The
mining the slope of var(t ) as a function of time, scaled durations of subsequent P-mode and S-mode segments
P
9
b′
might be correlated.
The first effect would lead to a smaller variance of
the duration of the P-mode segments and thereby to a
smaller variance of the total P-mode occupation times.
Thus if the dynamics destroys correlations sufficiently
rapidly, var(t ) should be smaller than expected theo-
P
retically. Our findings suggest that this is indeed the ϑ
case for the tetrahedron with inner sphere. The tetra-
d
hedron without inner sphere is ergodic but needs more
Θ
time to destroy correlations. It thus seems that the sec-
ondeffectrelatedtocorrelationstendstoenlargevar(t ).
P
For the rectangle with inner sphere, correlations are not
efficiently destroyed and the variance of the P-mode oc-
cupation times becomes much larger still. At last, we
have the rectangle without inner sphere, where var(t )
P b
scales with t2.
FIG. 10. Transformation of the ray widths under mode con-
VI. CONCLUSIONS version. Iftheexitangleisequaltotheangleofincidence(no
modeconversion),thewidthofthereflectedrayisunchanged.
Otherwise the widths are proportional to the cosines of the
We presentedsimulations of the propagationof elastic
angles ϑ and Θ.
wavesinthree-dimensionalbodiesofdifferentgeometries
inthe limitofclassicalrays. Becauseofmodeconversion
at the reflections on the surface of the bodies, one has ofthedecayfunctionofscatteringfidelitywhichhasbeen
to deal with an exponential proliferation of branches of showntofollowverycloselyuniversalrandommatrixpre-
the elastic ray, which is dealt with using Monte Carlo dictions. However, our results indicate that the pertur-
sampling. bationstrengthisanequallyinterestingquantity. Infact
Our simulations have shown that there is no unique it reveals much more system specific information, which
universalequilibriumdistributionforelasticrays,atleast may even have practicalapplications, for example in the
not for bodies with a sufficiently simple geometry. This analysis and verification of the properties of mechanical
is true even if the dynamics must be considered as com- components. It would be interesting to design new ex-
pletelychaotic. Forthetetrahedronwithinternalsphere, periments, which allowed for a quantitative comparison
the equilibrium limit of the S/P energy ratio was close between experiment and numerical simulations.
but not exactly equal to the theoretical value. Even
more surprisingly, we found that the homogeneous and
isotropic distribution of elastic rays with an S/P energy
ACKNOWLEDGMENTS
ratioequaltothetheoreticalequipartitionrationeednot
beanequilibriumdistributionatall. Whenchosenasthe
We are grateful to discussions with O. I. Lobkis and
initial condition of an ensemble of elastic rays, we found
R. L. Weaver at an early stage of the project. We are
thatinthecaseoftherectangularblocks,theS/P energy
also grateful to the hospitality of the Centro Interna-
ratio changes to a different value.
cional de Ciencias (Cuernavaca, M´exico) where this dis-
The main purpose of the present work was to analyse
cussionstookplace. Finally,weacknowledgethesupport
some of the conjectures made in Ref. [8] and to investi-
by CONACyT under grant number 129309.
gate, whether the observed deviations can be explained
withamodelbasedonclassicalrays. Whilemanyconjec-
turesandapproximationswereindeedjustified, themost
problematic assumption was that of random mode con- Appendix A: Energy conservation during ray
splitting
versions with given conversion rates. There, we found
two counteracting effects: One the one hand the length
distribution of individual P-mode and S-mode segments InSec.IIB,wediscussthesplittingofthetotalenergy
is not exponential due to the fact that mode conversions ofareflectedwavepacketintheraylimit. Inordertode-
cantakeplaceonlyduringreflections. Ontheotherhand, termine the energy share of the mode-converted branch,
there are correlations in time where the time scale de- we assumed energy conservation, which is an important
pends on the dynamics. In comparison to the random issue for the self-consistencyofour raymodel. It implies
mode conversion model, the former tends to reduce the that the reflection coefficients defined in Eqs. (2-6) are
variance of the P-mode occupation time while the latter related. Here,weshowexemplarilyforanincidentrayin
leads to an increase. P-mode that the total energy is indeed conserved.
Previouspublications[11,13]havefocusedontheform Theamplitude ofanelasticrayis physicallyrelatedto
10
the displacement of infinitesimal volume elements in the conserved, the following relation must hold:
solid. Therefore, amplitudes u or u (for P-mode and
P S ̺w2 ̺w2
S-mode, respectively) have units of length. The period c u (~r)2 dσ(~r)=
d P
averaged energy density of the respective wave field is 2 ZS 2
then given by
× c R2 u′ (~r)2dσ(~r)+c R2 u′′(~r)2dσ(~r)
(cid:20) d PP ZS′ P s SP ZS′′ P (cid:21)
̺w2 (A4)
2
ε = u , (A1)
P,S 2 {P,S} whereS isasurfaceperpendiculartotheincidentray,S′
isasurfaceperpendiculartothereflectedP-modebranch
where ̺ is the density of the medium and w the angular andS′′isasurfaceperpendiculartothereflectedS-mode
frequency of the wave (see Sec. 1.7 of [22]). To demon- branch. We have assumed that these surfaces are suffi-
strate the energy conservation, we show that the energy ciently to the reflection point such that the wave ampli-
flux is conserved in a quasi-stationary situation where a tudesmaybeconsideredconstantalongthelinesegments
verylongwavepacketisreflectedata freeplane surface. towards and away from the reflection point. The reflec-
It means that the energy flux along the propagation di- tionisschematicallydepictedinFig.10separatelyforthe
rection of the incident P-mode ray must be equal to the P-mode branch (left hand side) and the S-mode branch
sum of the fluxes along the reflectedP-mode andthe re- (right hand side). Without mode conversion, the geo-
flected and mode-converted S-mode branch: In general, metrical properties of the ray do not change, such that
the amount of energy flowing through a transversal sur- the amplitude u′ relative to S′ is just equal to u rela-
P P
face S is given by tive to S, which means that the corresponding integrals
coincide. With mode conversion, the width of the ray
increases by the factor cosΘ/cosϑ > 1. Thus, the inte-
F = ε(~r)v (~r)dσ(~r), (A2) gral over S′′ must be scaled by that factor. Therefore,
n
ZS dividing by the common factor ̺w2/2,
cosΘ
2 2 2
where v is the projectionof the ray velocity on the sur- c (1−R ) u (~r) dσ(~r)=c R
n d PP Z P s SP cosϑ
face normal. If we choose S to be normal to the propa- S
gation direction: 2
× u (~r) dσ(~r) (A5)
P
Z
S
Now, since the integrals are the same, we arrive at
̺w2
2
FP/S = 2 cd/sZSuP/S(~r) dσ(~r). (A3) 1−R2 = cs R2 cosΘ , (A6)
PP c SP cosϑ
d
Just as the propagation speed, also the energy flux de- which may be easily verified. Thus the energy flux is
pendsonthewavemode. Now,iftheenergyfluxisreally indeed conserved.
[1] R.Snieder, Pureappl. geophys. 163, 455 (2006) [12] T. Gorin, T. Prosen, T. H. Seligman, and M. Zˇnidariˇc,
[2] R.Courtland, Nature453, 146 (2008) Phys. Rep.435, 33 (2006)
[3] E. Larose, J. De Rosny, L. Margerin, D. Anache, [13] O.I.LobkisandR.L.Weaver,Phys. Rev.E78,066212
P. Gouedard, M. Campillo, and B. van Tiggelen, (2008)
Phys.Rev.E 73, 016609 (2006) [14] R. L. Weaver,J. Acoust. Soc. Am.71, 1608 (1982)
[4] A.Kohler and R. Blu¨mel, Ann.Phys.267, 249 (1998) [15] A. D. Pierce, Acoustics: An Introduction (McGraw-Hill,
[5] M. V.Berry, Proc. R. Soc. Lond.A 400, 229 (1985) 1981)
[6] P. Jacquod, P. G. Silvestrov, and C. W. J. Beenakker, [16] K.F.Graff,WaveMotioninElasticSolids (Dover,1991)
Phys.Rev.E 64, 055203(R) (2001) [17] H. Schomerus and P. Jacquod, J. Phys. A: Math. Gen.
[7] N. R. Cerruti and S. Tomsovic, Phys. Rev. Lett. 88, 38, 10663 (2005)
054103 (2002) [18] M. C. Gutzwiller, Chaos in classical and quantum me-
[8] O. I. Lobkis and R. L. Weaver, Phys. Rev. Lett. 90, chanics (Springer, 1990)
254302 (2003) [19] L. Kaplan, New J. Phys.4, 90 (2002)
[9] R. Sch¨afer, T. Gorin, H.-J. St¨ockmann, and T. H. Selig- [20] E. Bogomolny, N. Pavloff and C. Schmit, Phys. Rev. E
man, NewJ. Phys.7, 152 (2005) 61, 3689 (2000)
[10] R. Sch¨afer, H.-J. St¨ockmann, T. Gorin, and T. H. Selig- [21] J.S.Hersch,M.R.HaggertyandE.J.Heller,Phys.Rev.
man, Phys.Rev.Lett. 95, 184102 (2005) E 62, 4873 (2000)
[11] T.Gorin,T. H.Seligman,andR.L.Weaver,Phys.Rev. [22] J. D. Achenbach, Wave propagation in elastic solids,
E 73, 015202(R) (2006) Applied mathematics and mechanics, Vol. 16 (North-