text:
# Ultrafast spectroscopy of propagating coherent acoustic phonons in GaN/InGaN heterostructures

## I Introduction

Heterostructures of GaN and InGaN are important materials owing to their applications to blue laser diodes and high power electronics.Nakamura et al. (2000) Strong coherent acoustic phonon oscillations have recently been detected in InGaN/GaN multiple quantum wells.Sun et al. (2000); Özgür et al. (2001a) These phonon oscillations were much stronger than folded acoustic phonon oscillations observed in other semiconductor superlattices.Yamamoto et al. (1994); Bartels et al. (1999); Mizoguchi et al. (1999)

InGaN/GaN heterostructures are highly strained at high In concentrations giving rise to large built-in piezoelectric fields,Lefebvre et al. (2001); Takeuchi et al. (1998); Im et al. (1998); Chichiba et al. (1998) and the large strength of the coherent acoustic phonon oscillations was attributed to the large strain and piezoelectric fields.Sun et al. (2000)

In this paper, we report the generation of strong localized coherent phonon wavepackets in strained layer InxGa1-xN/GaN epilayers and heterostructures grown on GaN and Sapphire substrates.Yahng et al. (2002) By focusing high repetition rate, frequency-doubled femtosecond Ti:Sapphire laser pulses onto strongly strained InGaN/GaN heterostructures, we can, through the carrier-phonon interaction, generate coherent phonon wavepackets which are initially localized near the epilayer/surface but then propagate away from the surface/epilayer and through a GaN layer. As the wavepackets propagate, they modulate the local index of refraction and can be observed in the time-dependent differential reflectivity of the probe pulse. There is a sudden drop in the amplitude of the reflectivity oscillation of the probe pulse when the phonon wavepacket reaches a GaN-sapphire interface below the surface. Theoretical calculations as well as experimental evidences support this picture; the sudden drop of amplitude when the wave encounters the GaN-sapphire interface cannot be explained if the wavepacket had large spatial extent. When the wavepacket encounters the GaN-sapphire interface, part of the wave gets reflected while most of it gets transmitted into the sapphire substrate, depending on the interface properties and the excess energy of the exciting photons. This experiment illustrates a non-destructive way of generating high pressure tensile waves in strained heterostructures and using them to probe semiconductor structure below the surface of the sample. Since the strength of this non-destructive wave is determined by the strain between GaN and InGaN, it is likely that even stronger coherent phonons can be generated in InGaN/GaN digital alloys grown on a GaN substrate.

## II Experimental Results

In the experiments, frequency-doubled pulses of mode-locked Ti:sapphire lasers are used to perform reflective pump-probe measurements on four different sample types which are shown in Fig. 1): Type (I) InGaN epilayers; Type (II) InGaN/GaN double quantum wells (DQWs); Type (III) InGaN/GaN single quantum well (SQW); and Type (IV) InGaN/GaN light-emitting diode (LED) structures. The peak pump power is estimated to be 400 MW/cm2, corresponding to a carrier density of 1019 cm-3 and the doubled pulse width is 250 fs. All samples were grown on a $`c`$-plane sapphire substrate by metal organic chemical vapor deposition. The InGaN epilayers shown in Fig. 1(a) consist of 1 $`\mu `$m GaN grown on a sapphire substrate and capped with 30 nm of InxGa1-xN with In composition, $`x`$, varying from 0.04 to 0.12. The DQW sample shown in Fig. 1(b) consists of GaN (1 $`\mu `$m), double quantum wells of In0.12Ga0.88N (1–16 nm), a barrier of GaN (1–16 nm), and a GaN cap layer (0.1 $`\mu `$m). The SQW sample in Fig. 1(c) consists of GaN (2 $`\mu `$m), an In0.12Ga0.88N well (24 nm), and a GaN cap layer (0.1 $`\mu `$m). The blue LED structure shown in Fig. 1(d) consists of $`n`$-GaN (4 $`\mu `$m), 5 quantum wells of In0.15Ga0.85N (2 nm) and 4 barriers of GaN (10 nm), $`p`$-Al0.1Ga0.9N (20 nm), and $`p`$-GaN (0.2 $`\mu `$m).

Differential reflection pump-probe measurements for the InxGa1-xN epilayers are shown in Fig. 1(a). Figure 2 shows the oscillatory component of the measured probe differential reflectivity for the InGaN epilayers (type I) with various In concentrations, x. For comparison purposes, we performed differential reflectivity measurements on a pure GaN HVPE grown sample in order to show that no differential reflectivity oscillations are present in the absence of strain and an epilayer. The energy of the pump laser was varied between 3.22 and 3.35 eV, to keep the excess carrier energy above the InGaN band gap but below the GaN band gap. We note that if the laser energy was below the InGaN band gap, no signal was detected. Therefore, carrier generation is essential to observing the oscillations, unlike a recent coherent optical phonon experiment in GaNYee et al. (2002). The inset shows the pump-probe signal prior to the background subtraction for $`x=0.12`$. The background results from the relaxation of the photoexcited electrons and holes. The oscillations are quite large, on the order of $`10^2`$$`10^3`$ and the period is 8–9 ps, independent of the In composition but dependent on the probe photon energy. The amplitude of the oscillation is approximately proportional to the In concentration indicating that the strain at the InGaN/GaN interface is important. The observed period is approximately $`\tau =\lambda /2C_\text{s}n`$, where $`\lambda `$ is the probe beam wavelength, $`C_\text{s}`$ the longitudinal acoustic sound velocity, and $`n`$ the refractive index of GaN Thompsen et al. (1986).

We performed two-color pump-probe experiments for a type III InGaN SQW sample as shown in Fig. 1(b). Fig. 3 shows the differential reflectivity oscillations for different probe energies. Note that the period of the oscillation changes and is proportional to the probe wavelength. In addition, the amplitude of the differential reflectivity oscillation decreases as the detuning (with respect to the pump) becomes larger. The inset shows the oscillation amplitude as a function of the probe energy in a logarithmic scale and there is an $``$2-order-of magnitude decrease in differential reflectivity when the probe energy changes from 3.26 eV to 1.63 eV.

Interestingly, Fig. 4 shows that the amplitude of the oscillatory component of the differential reflectivity is independent of the bias voltage, even though the carrier lifetime changes dramatically with voltage bias. Fig. 4 shows the bias dependent acoustic phonon differential reflectivity oscillations in a type IV blue LED structure (see Fig. 1(d)) at a pump energy of 3.17 eV. The lifetime of the background signal drastically decreases as the bias increases as shown in Fig. 4(a). This is due to the carrier recombination time and the decrease in the tunneling escape time in the strong external bias regimeJho et al. (2001). On the other hand, the amplitude and frequency of the oscillatory component of the differential reflectivity doesn’t change much with bias voltage \[Fig. 4(b)\]. Since the observed reflectivity oscillation is independent of the carrier lifetime for lifetimes as short as 1 ps due to ultrafast tunneling \[bottom curve of Fig. 4(a)\], it implies that once the source that modulates the experimentally observed reflectivity is launched by the sub-picosecond generation of carriers, the remaining carriers do little to affect the source. This suggests that the reflectivity oscillation is due to the strain pulse which is generated at short times once the pump excites the carriers and modulates the lattice constant.

Figure 5 shows the well- and barrier-dependent acoustic phonon differential reflectivity oscillations in the type II DQW samples (see Fig. 1(b)) at a pump energy of 3.22 eV. The amplitude increases as the well width increases. However, the oscillation amplitude of the differential reflectivity doesn’t change much with the barrier width. This means that the generation of the acoustic phonons is due to the InGaN well and not the GaN barrier. This also verifies that the oscillation is due to the strain in InGaN layer.

Interesting results are seen in the long time behavior of the reflectivity oscillations shown in Fig. 6(a). The long time scale reflectivity oscillation is plotted for the epilayer (I), DQW (II), SQW (III), and the blue LED structure (IV) at 3.29 eV (below the GaN band gap). Astonishingly, the oscillation amplitude abruptly decreases within one cycle of an oscillation at a critical time which appears to scale with the thickness of the GaN layer in each sample. In addition, the slope of the GaN thickness vs. the critical time is very close to the known value of the sound velocity in GaN \[inset of Fig. 6(a)\] Sun et al. (2000). Fig. 6(b) shows results when the probe laser energy is changed to 3.44 eV which is above the GaN band gap. Then the laser probe is sensitive to coherent phonon oscillations only within an absorption depth of the surface. We see that the amplitude of the reflectivity oscillation exponentially decays with a decay time of 24.2 ps corresponding to a penetration depth in GaN of about 0.17 micron (=24.2 ps$`\times `$7000 m/s). The oscillation reappears at 260 ps for epilayer (I) and 340 ps for DQW (II). This is twice the critical time for the oscillations to disappear when the photon energy is 3.29 eV. This further shows that the probe pulse is sensitive to the coherent acoustic wave. The ”echo” in the probe signal results from the partial reflection of the coherent phonon off the GaN/sapphire interface.

## III Theory

To explain the experimental results discussed in the last section, we have developed a theoretical model of the generation, propagation and detection of coherent acoustic phonons in strained GaN/InGaN heterostructures. The pump laser pulse generates a strain field that propagates through the sample which, in turn, causes a spatio-temporal change in the index of refraction. This change is responsible for the oscillatory behavior seen in the probe-field reflectivity in various semiconductor heterostructures. An approximate method of solving Maxwell’s equations in the presence of spatio-temporal disturbances in the optical properties and obtaining the reflectivity of the probe field in thin films excited by picosecond pump pulses can be found in Thompsen et. al.Thompsen et al. (1986)

The spatio-temporal disturbance of the refractive index is caused by the propagating coherent phonon wavepackets. Thus, an essential ingredient in understanding the probe reflectivity is a model for the generation and propagation of the very short strain pulse in the sample. Recently, a microscopic theory explaining the generation and propagation of such strain pulses was reported by Sanders et. al.Sanders et al. (2001, 2002) where it was shown that propagating coherent acoustic phonon wavepackets are created by the nonequilibrium carriers excited by the ultrafast pump pulse. The acoustic phonon oscillations arise through the electron-phonon interaction with the photoexcited carriers. Both acoustic deformation potential and piezoelectric scattering were considered in the microscopic model. It was found that under typical experimental conditions, the microscopic theory could be simplified and mapped onto a loaded string model. Here, we use the string model of coherent phonon pulse generation to obtain the strain field seen by the probe pulse.

First, we solve Maxwell’s equations to obtain the probe reflectivity in the presence of a generalized spatio-temporal disturbance of the index of refraction. Let $`n_\text{b}`$ be the index of refraction without the strain which is real because initially the absorption can be neglected and let $`\delta \stackrel{~}{n}=\delta n+i\kappa `$ be the propagating change in the index of refraction due to the strain. When the effect of the change of the index of refraction is taken into account, the probe field with energy $`\omega `$ can be described by the following generalized wave equation

$$\frac{^2E(z,t)}{z^2}+\frac{\omega ^2}{c^2}[n_\text{b}+\delta \stackrel{~}{n}(z,t)]^2E(z,t)=0$$
(1)
where $`E(z,t)`$ is the probe field in the slowly varying envelope function approximation and $`\omega `$ is the central frequency of the probe pulse. Eqn. (1) is obtained from Maxwell’s equations assuming that the polarization response is instantaneous and that the probe pulse obeys the slowly varying envelope function approximation. Since $`|\delta \stackrel{~}{n}|n_\text{b}`$ under typical conditions, Eqn. (1) can be cast into

$$\frac{^2E(z,t)}{z^2}+(n_\text{b}k)^2E(z,t)=2n_\text{b}k^2\delta \stackrel{~}{n}(z,t)E(z,t)$$
(2)
where $`k=\omega /c`$ is the probe wavevector. To relate the change of the index of refraction to the strain field, $`\eta (z,t)U(z,t)/z`$, we assume $`|\delta \stackrel{~}{n}|n_\text{b}`$ and adopt the linear approximation Thompsen et al. (1986)

$$\delta \stackrel{~}{n}(z,t)=\frac{\stackrel{~}{n}}{\eta }\eta (z,t).$$
(3)
We view Eqn. (2) as an inhomogeneous Helmholtz equation and obtain the solution using the Green’s function technique. The desired Green’s function is determined by solving

$$\frac{^2G(z,z^{})}{z^2}+(n_\text{b}k)^2G(z,z^{})=\delta (zz^{})$$
(4)
and the result is

$$G(z,z^{})=\frac{i}{2n_\text{b}k}\mathrm{exp}\left(in_\text{b}k|zz^{}|\right).$$
(5)
Then, the solution to Eqn. (2) can be written as

$`E(z,t)`$ $`=`$ $`E_\text{h}(z,t)+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z^{}G(z,z^{})\left\{2n_\text{b}k^2\delta \stackrel{~}{n}(z^{},t)E(z^{},t)\right\}`$ (6)
$``$ $`E_\text{h}(z,t)+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z^{}G(z,z^{})\left\{2n_\text{b}k^2\delta \stackrel{~}{n}(z^{},t)E_\text{h}(z^{},t)\right\}`$ (7)
where we have chosen the lowest Born series in the last line. In Eqn. (7), $`E_\text{h}`$ is the homogeneous solution which takes the form $`E_\text{h}(z,t)=\stackrel{~}{E}_\text{h}(z,t)\mathrm{exp}\left\{i(n_\text{b}kz\omega t)\right\}`$ where $`\stackrel{~}{E}_\text{h}(z,t)`$ is a slowly varying envelope function. This describes the probe pulse moving to the right in the sample without optical distortion.

We now apply this approximate solution to our structure where the interface between air and the sample is chosen at $`z=0`$. In the air, where $`z0`$, there is an incident probe pulse traveling toward the sample as well as a reflected pulse. The electric field in the air, $`E_<(z,t)`$, can thus be written as the sum

$$E_<(z,t)=\stackrel{~}{E}_\text{i}(z,t)e^{i(kz\omega t)}+\stackrel{~}{E}_\text{r}(z,t)e^{i(kz+\omega t)}$$
(8)
where $`\stackrel{~}{E}_\text{i}(z,t)`$ and $`\stackrel{~}{E}_\text{r}(z,t)`$ are the slowly varying envelope functions of the incident and reflected probe fields, respectively. Inside the sample, $`z0`$, the solution is given as

$`E_>(z,t)`$ $`=`$ $`E_\text{h}(z,t)+{\displaystyle _0^{\mathrm{}}}𝑑z^{}\mathrm{exp}\left(in_\text{b}k(z^{}z)\right)ik\delta \stackrel{~}{n}(z^{},t)E_\text{h}(z^{},t)`$ (9)
$`=`$ $`\stackrel{~}{E}_\text{t}(z,t)e^{i(n_\text{b}kz\omega t)}+\left\{{\displaystyle _0^{\mathrm{}}}𝑑z^{}\mathrm{exp}\left(2in_\text{b}kz^{}\right)ik\delta \stackrel{~}{n}(z^{},t)\right\}\stackrel{~}{E}_\text{t}(z,t)e^{i(n_\text{b}kz+\omega t)}`$
$`=`$ $`\stackrel{~}{E}_\text{t}(z,t)e^{i(n_\text{b}kz\omega t)}+𝒜(n_\text{b}k,t)\stackrel{~}{E}_\text{t}(z,t)e^{i(n_\text{b}kz+\omega t)}`$ (10)
where we used $`E_\text{h}(z,t)=\stackrel{~}{E}_\text{t}(z,t)\mathrm{exp}(i(n_\text{b}kz\omega t))`$ in the second step, assuming that $`\stackrel{~}{E}_\text{t}(z,t)`$ is nearly constant within the slowly varying envelope function approximation, and we define a reflected amplitude function

$$𝒜(n_\text{b}k,t)_0^{\mathrm{}}𝑑z^{}\mathrm{exp}\left(2in_\text{b}kz^{}\right)ik\delta \stackrel{~}{n}(z^{},t).$$
(11)
The expression for the reflected amplitude function, $`𝒜(n_\text{b}k,t)`$, in Eqn. (9) says there is a frequency-dependent modulation of the amplitude in the reflected wave in the sample due to the propagating strain.

Having determined the waves on both sides of the interface, we can now calculate the reflectivity. We apply the usual boundary conditions to the slowly varying envelope functions and the results are written compactly as

$$\left(\begin{array}{cc}1& 1+𝒜\\ 1& n_\text{b}(1𝒜)\end{array}\right)\left(\begin{array}{c}\stackrel{~}{E}_\text{r}\\ \stackrel{~}{E}_\text{t}\end{array}\right)=\left(\begin{array}{c}\stackrel{~}{E}_\text{i}\\ \stackrel{~}{E}_\text{i}\end{array}\right).$$
(12)
We solve this equation to obtain

$$\frac{\stackrel{~}{E}_\text{r}}{\stackrel{~}{E}_\text{i}}=\frac{r_0+𝒜}{1+r_0𝒜}r_0+𝒜$$
(13)
where $`r_0=(1n_\text{b})/(1+n_\text{b})`$. To the same order, we find that $`\stackrel{~}{E}_\text{t}/\stackrel{~}{E}_it_0(1r_0𝒜)`$ where $`t_0=2/(1+n_\text{b})`$. It is now straightforward to calculate the differential reflectivity as

$$\frac{\mathrm{\Delta }R}{R}=\frac{|r_0+\mathrm{\Delta }r|^2|r_0|^2}{|r_0|^2}\frac{2}{r_0}\mathrm{Re}𝒜.$$
(14)
Finally, by substituting the linear law Eqn. (3) into Eqn. (11) and using Eqn. (14) we get

$$\frac{\mathrm{\Delta }R}{R}=_0^{\mathrm{}}𝑑z(z,\omega )\frac{U(z,t)}{z}$$
(15)
where the sensitivity function, $`(z,\omega )`$, is defined as

$$(z,\omega )=\frac{2k}{r_0}\left[\frac{n}{\eta }\mathrm{sin}(2n_\text{b}kz)+\frac{\kappa }{\eta }\mathrm{cos}(2n_\text{b}kz)\right]$$
(16)
In Eqn. (15), the differential reflectivity is expressed in terms of the lattice displacement, $`U(z,t)`$, due to propagating coherent phonons. Sanders et. al.Sanders et al. (2001, 2002) developed a microscopic theory showing that the coherent phonon lattice displacement satisfies a driven string equation,

$$\frac{^2U(z,t)}{t^2}C_\text{s}^2\frac{^2U(z,t)}{z^2}=S(z,t),$$
(17)
where $`C_\text{s}`$ is the LA sound speed in the medium and $`S(z,t)`$ is a driving term which depends on the photogenerated carrier density. The LA sound speed is related to the mass density, $`\rho `$, and the elastic stiffness constant, $`C_{33}`$, by $`C_\text{s}=\sqrt{C_{33}/\rho }`$. The LA sound speed is taken to be different in the GaN/InGaN heterostructure and Sapphire substrate. For simplicity, we neglect the sound speed mismatch between the GaN and InxGa1-xN layers.

The driving function, $`S(z,t)`$, is nonuniform and is given by

$$S(z,t)=\underset{\nu }{}S_\nu (z,t),$$
(18)
where the summation index, $`\nu `$, runs over carrier species, i.e., conduction electrons, heavy holes, light holes, and crystal field split holes, that are created by the pump pulse. Each carrier species makes a separate contribution to the driving function. The partial driving functions, $`S_\nu (z,t)`$, in piezoelectric wurtzite crystals depend on the density of the photoexcited carriers. Thus,

$$S_\nu (z,t)=\pm \frac{1}{\rho }\left\{a_\nu \frac{}{z}+\frac{4\pi |e|e_{33}}{ϵ_{\mathrm{}}}\right\}\rho _\nu (z,t),$$
(19)
where the plus sign is used for conduction electrons and the minus sign is used for holes. Here $`\rho _\nu (z,t)`$ is the photogenerated electron or hole number density, which is real and positive, $`\rho `$ is the mass density, $`a_\nu `$ are the deformation potentials for each carrier species, $`e_{33}`$ is the piezoelectric constant, and $`ϵ_{\mathrm{}}`$ is the high frequency dielectric constant.

To be more specific, we will consider a SQW sample of the type shown in Fig. 1(c). We adopt a simple model for photogeneration of electrons and holes in which the photogenerated electron and hole number densities are proportional to the squared ground state particle in a box wavefunctions. The exact shape of the electron and hole number density profile is not critical in the present calculation since all that really matters is that the electrons and holes be strongly localized. The carriers are assumed to be instantaneously generated by the pump pulse and are localized in the InxGa1-xN quantum well. Having obtained a model expression for $`\rho _\nu (z,t)`$, it is straightforward to obtain $`S(z,t)`$ using Eqns. (18) and (19).

To obtain $`U(z,t)`$, we solve driven string equations in the GaN epilayer and the Sapphire substrate, namely

$$\frac{^2U(z,t)}{t^2}C_0^2\frac{^2U(z,t)}{z^2}=S(z,t)(0zL)$$
(20a)
and
$$\frac{^2U(z,t)}{t^2}C_1^2\frac{^2U(z,t)}{z^2}=0(LzZ_\text{s})$$
(20b)
where $`C_0`$ and $`C_1`$ are LA sound speeds in the GaN and Sapphire substrate, respectively. In Eqn. (20b), the Sapphire substrate has finite thickness. To simulate coherent phonon propagation in an infinite Sapphire substrate, $`Z_\text{s}`$ in Eqn. (20b) is chosen large enough so that the propagating sound pulse generated in the GaN epilayer does not have sufficient time to reach $`z=Z_\text{s}`$ during the simulation. If $`T_{sim}`$ is the duration of the simulation, this implies $`Z_\text{s}L+C_1T_{sim}`$.

Equations (20a) and (20b) are solved subject to initial and boundary conditions. The initial conditions are

$$U(z,0)=\frac{U(z,0)}{t}=0.$$
(21)
At the GaN-air interface at $`z=0`$, we assume the free surface boundary condition

$$\frac{U(0,t)}{z}=0$$
(22a)
since the air exerts no force on the GaN epilayer. The phonon displacement and the force per unit area are continuous at the GaN-Sapphire interface so that
$$U(Lϵ,t)=U(L+ϵ,t)$$
(22b)
and
$$\rho _0C_0^2\frac{U(Lϵ,t)}{z}=\rho _1C_1^2\frac{U(L+ϵ,t)}{z}.$$
(22c)
The boundary condition at $`z=Z_\text{s}`$ can be chosen arbitrarily since the propagating sound pulse never reaches this interface. For mathematical convenience, we choose the rigid boundary condition
$$U(Z_\text{s},t)=0.$$
(22d)
To obtain $`U(z,t)`$ for general $`S(z,t)`$, we first need to find the harmonic solutions in the absence of strain, i.e. $`S(z,t)=0`$. The harmonic solutions are taken to be

$$U_n(z,t)=W_n(z)e^{i\omega _nt}(\omega _n0)$$
(23)
and it is easy to show that the mode functions, $`W_n(z)`$, satisfy

$$\frac{d^2W_n(z)}{dz^2}+\frac{\omega _n^2}{C_0^2}W_n(z)=0(0zL)$$
(24a)
and
$$\frac{d^2W_n(z)}{dz^2}+\frac{\omega _n^2}{C_1^2}W_n(z)=0(LzZ_\text{s})$$
(24b)
Applying the boundary conditions (22) we obtain the mode functions

$$W_n(z)=\{\begin{array}{cc}\mathrm{cos}\left(\omega _nz/C_0\right)\hfill & \text{if }0zL\hfill \\ B_n\mathrm{sin}\left(\omega _n(Z_\text{s}z)/C_1\right)\hfill & \text{if }LzZ_\text{s}\hfill \end{array}$$
(25a)
with
$$B_n=\frac{\mathrm{cos}\left(\omega _nL/C_0\right)}{\mathrm{sin}\left(\omega _n(Z_\text{s}L)/C_1\right)}.$$
(25b)
The mode frequencies, $`\omega _n`$, are solutions of the transcendental equation

$$\frac{1}{\rho _0C_0}\mathrm{cot}\left(\frac{\omega _nL}{C_0}\right)=\frac{1}{\rho _1C_1}\mathrm{tan}\left(\frac{\omega _n(Z_\text{s}L)}{C_1}\right)$$
(26)
which we solve numerically to obtain the mode frequencies, $`\omega _n`$ ($`n=0,1,2,\mathrm{}`$). The index, $`n`$, is equal to the number of nodes in the mode functions, $`W_n(z)`$.

A general displacement can be expanded in terms of the harmonic modes as

$$U(z,t)=\underset{n=0}{\overset{\mathrm{}}{}}q_n(t)W_n(z).$$
(27)
Substituting equation (27) for $`U(z,t)`$ into (20) and taking the initial conditions (21) into account, we find that the expansion coefficients, $`q_n(t)`$, satisfy a driven harmonic oscillator equation

$$\ddot{q}_n(t)+\omega _n^2\dot{q}_n(t)=Q_n(t),$$
(28)
subject to the initial conditions $`q_n(0)=\dot{q}_n(0)=0`$. The harmonic oscillator driving term $`Q_n(t)`$ is given by

$$Q_n(t)=\frac{_0^{Z_\text{s}}𝑑zW_n(z)S(z,t)}{_0^{Z_\text{s}}𝑑zW_n(z)^2}.$$
(29)
In our simple displacive model for photogeneration of carriers, $`S(z,t)=S(z)\mathrm{\Theta }(t)`$ where $`\mathrm{\Theta }(t)`$ is the Heaviside step function. In this case, the lattice displacement is explicitly given by

$$U(z,t)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{S_n}{\omega _n^2}(1\mathrm{cos}(\omega _nt))W_n(z)$$
(30)
with $`S_n`$ defined as

$$S_n=\frac{_0^{Z_\text{s}}𝑑zW_n(z)S(z)}{_0^{Z_\text{s}}𝑑zW_n(z)^2}.$$
(31)
Using the lattice displacement (30), we obtain the time-dependent differential reflectivity at the probe frequency, $`\omega `$, from equation (15). The result is

$$\frac{\mathrm{\Delta }R}{R}(\omega ,t)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{S_n}{\omega _n^2}(1\mathrm{cos}(\omega _nt))R_n(\omega )$$
(32)
where

$$R_n(\omega )=_0^{Z_\text{s}}𝑑z(z,\omega )\frac{dW_n(z)}{dz}.$$
(33)
can be evaluated analytically.

With the above formalism, we solve for the lattice displacement, $`U(z,t)`$, for a coherent LA phonon pulse propagating in a multilayer structure consisting of a 1.124 $`\mu `$m thick GaN epilayer grown on top of an infinitely thick Sapphire substrate with the growth direction along $`z`$. We take the origin to be at the GaN-air interface and the GaN-Sapphire interface is taken to be at $`z=L=1.124`$ $`\mu `$m. We assume that carriers are photogenerated in a single 240 Å thick InxGa1-xN quantum well embedded in the GaN layer 0.1 $`\mu `$m below the GaN-air interface and 1 $`\mu `$m above the Sapphire substrate. Our structure thus resembles the SQW sample shown in Fig. 1(c). In the GaN epilayer, we take $`C_{33}=379`$ GPa and $`\rho _0=6.139`$ gm/cm3 Ambacher et al. (1999) from which we obtain $`C_0=7857`$ m/s. For the Sapphire substrate, we take $`C_{33}=500`$ GPa and $`\rho _1=3.986`$ gm/cm3 Every et al. (1984) from which we find $`C_1=11200`$ m/s.

The results of our simulation are shown in Fig. 7. A contour map of the strain, $`U(z,t)/z`$, is shown in Fig. 7(a). We plot the strain as a function of the depth below the surface and the probe delay time. Photoexcitation of electrons and holes in the InGaN quantum well generates two coherent LA sound pulses traveling in opposite directions. The pulses are totally reflected off the GaN/air interface at $`z=0`$ and are partially reflected at the Sapphire substrate at $`z=1.124\mu `$m. Approximately 95% of the pulse energy is transmitted and only 5% is reflected at the substrate. The speed of the LA phonon pulses is just the slope of the propagating wave trains seen in Fig. 7(a) and one can clearly see that the LA sound speed is greater in the Sapphire substrate.

From the strain the differential reflectivity can be obtained from equation (15). From Fig. 7(a), the strains, $`U(z,t)/z`$, associated with the propagating pulses are highly localized and travel at the LA sound speed. Each pulse contributes a term to the differential reflectivity that goes like

$$\frac{\mathrm{\Delta }R}{R}(\omega ,t)(C_0t,\omega )\frac{\omega }{c}\mathrm{sin}\left(\frac{2n_\text{b}\omega }{c}C_0t+\varphi \right)$$
(34)
The period of the oscillations of $``$ depends on the probe wavelength, $`\lambda =2\pi c/\omega `$, with the result that the observed differential reflectivity oscillates in time with period, $`T=\pi c/(n_\text{b}C_0\omega )=\lambda /(2n_\text{b}C_0)`$, where $`n_\text{b}=2.4`$ is the index of refraction, and $`C_0=7857`$ m/s is the LA sound speed in GaN. For $`\lambda =377`$ nm ($`\mathrm{}\omega =3.29`$ eV), this gives us $`T=10`$ ps. The sensitivity function, $`(z,\omega )`$, defined in equation (16) is an oscillating function in the GaN/InGaN epilayer and is assumed to vanish in the Sapphire substrate. Our computed differential reflectivity is shown in Fig. 7(b) for a probe wavelength of $`\lambda =377\text{nm}`$. We find that the reflectivity abruptly attenuates when the strain pulse enters the Sapphire substrate at $`t=170`$ ps. The reflected strain pulses give rise to the weaker oscillations seen for $`t>170`$ ps. These oscillations are predicted to continue until the reflected pulses are again partially reflected off the Sapphire substrate at $`t=430`$ ps.

## IV Simple model

Since the coherent oscillation observed in the differential reflectivity stems essentially from the strain pulse propagating into the layers, most phenomena can be understood by a simple macroscopic model that is presented in this section.

Instead of solving the loaded string equations for the strain to obtain a propagating disturbance in the refractive index, the propagating strain pulse at a given moment can be viewed as a thin strained layer in the sample, where the index of refraction is assumed to be slightly different from the rest of the sample. This situation is schematically depicted in Fig. 8 where a fictitious, thin GaN strained layer is located at $`z`$ in the thick host GaN layer. The thickness of the strained layer, $`d`$, is approximately the width of the traveling coherent phonon strain field, $`U(z,t)/z`$, and is to be determined from the microscopic theory. From the last section, it was seen that the propagating strain field is strongly localized so that $`d`$ is small. In the example of the last section, $`d`$ is approximately equal to the quantum well width. Here we assume the strain pulse has been already created near the air/GaN interface and do not consider its generation procedure. We treat $`d`$ as a phenomenological constant and also assume that the change in the index of refraction is constant. This strained GaN layer travels into the structure with the speed of the acoustic phonon wavepacket $`C_0=7\times 10^3m/s`$, so the location of the stained layer is given as $`z=C_0\tau `$ where $`\tau `$ is the pump-probe delay time.

Within the slowly varying envelope function approximation, the solutions to the Maxwell equation can be written as plane waves in each region,

$$E_i(z,t)=a_ie^{ik_iz_ii\omega t}+b_ie^{ik_iz_ii\omega t}$$
(35)
where $`E_i`$ is the electric field in the left layer of $`i`$th interface, and $`a_i`$ and $`b_i`$ are the slowly varying amplitudes. The magnetic field is given by $`B_iE_i/z`$. For a normally incident probe field, we apply the usual matching conditions on $`E`$ and $`B`$ to obtain

$$\left(\begin{array}{c}a_i\\ b_i\end{array}\right)=_i\left(\begin{array}{c}a_{i+1}\\ b_{i+1}\end{array}\right),$$
(36)
where the transfer matrix, $`_i`$, is given explicitly as

$`_i=`$ (37)
$`{\displaystyle \frac{1}{2}}\left(\begin{array}{cc}(1+\frac{k_{i+1}}{k_i})e^{i(k_{i+1}k_i)z_i}& (1\frac{k_{i+1}}{k_i})e^{i(k_{i+1}+k_i)z_i}\\ (1\frac{k_{i+1}}{k_i})e^{i(k_{i+1}+k_i)z_i}& (1+\frac{k_{i+1}}{k_i})e^{i(k_{i+1}k_i)z_i}\end{array}\right).`$
To apply this formula to our configuration, we normalize the incident amplitude to $`1`$, let $`r`$ be the reflected amplitude in the air, and impose the boundary condition that there is only a transmitted wave into the Sapphire substrate with amplitude $`t`$ and no reflected wave from the GaN-Sapphire interface back into the GaN epilayer. The latter assumption is reasonable since the microscopic theory of the previous section shows that only 5% of the pulse energy is reflected from the interface between the GaN and the Sapphire substrate. The total reflection and transmission amplitudes $`r`$ and $`t`$ for the GaN epilayer structure are determined by

$$\left(\begin{array}{c}1\\ r\end{array}\right)=\left(\begin{array}{cc}_{11}& _{12}\\ _{21}& _{22}\end{array}\right)\left(\begin{array}{c}t\\ 0\end{array}\right)$$
(40)
where $`=_1_2\mathrm{}_n`$. The reflection amplitude is readily found to be

$$r=_{21}/_{11}=r_0+\mathrm{\Delta }r,$$
(41)
where $`r_0`$ is the background contribution without the strained layer. The total transmission amplitude is given as $`t=1/_{11}`$. We can now numerically calculate the differential reflectivity by substituting $`r`$ into Eqn. (14).

In Fig. 9 we compare our theoretical model with experimental results. In Fig. 9(a), the oscillatory part of the differential reflectivity is plotted as a function of probe delay for three values of the probe photon energy. This is the same data that was shown in Fig. 3 where it was seen that the oscillation period of the coherent phonon reflectivity oscillations are proportional to the probe wavelength. In Fig. 9(b), we plot the corresponding theoretical differential reflectivities obtained from Eqn. (14). Comparing Figs. 9(a) and 9(b), we see that our geometrical optics model successfully explains the observed relation between the coherent phonon oscillation period and the probe wavelength. However, the calculated amplitudes of oscillation is inconsistent with the experimental data because in this case we have not taken into account of the change of index of refraction with respect to the probe wavelength which we will do later in the paper.

In Fig. 10 we plot the differential reflectivity as a function of the probe delay for three different values of the change in the index of refraction in the strained layer, $`\delta \stackrel{~}{n}`$. As expected, the greater the change in index of refraction, the greater the amplitude of the differential reflectivity oscillations. A larger change in the index of refraction implies that more electron-hole pairs are excited near the air-GaN interface, which acts as a stronger source for the coherent acoustic phonon reflectivity oscillations. These results are qualitatively consistent with the experimental results shown in Fig. 2. In Fig. 11 we fix the change in the index of refraction, $`\delta n=0.01`$, and varied the thickness of the strained layer, $`d`$. The result shows a larger amplitude for the differential reflection in wider strained layers, which is consistent with what was experimentally observed in Fig. 5.

Both Fig. 10 and Fig. 11 suggest that the amplitude of the differential reflectivity seems to increase monotonically with $`\delta n`$ and $`d`$. All these features can be understood more easily using the single-reflection approximation to the full formula for the reflectivity in Eqn. (41). Assuming the change in the index of refraction is small, $`|\delta \stackrel{~}{n}|1`$, we may select contributions from only terms proportional to $`\delta \stackrel{~}{n}`$ in the infinite Fabry-Perot series for the total reflection amplitude. The relevant reflection processes selected are schematically depicted in Fig. 12. In this case the total reflection amplitude is given by the leading terms in the Fabry-Perot series

$$r=r_0+r_1+r_2+𝒪\left(\delta \stackrel{~}{n}^2\right)$$
(42)
where $`r_0`$ is the background reflection amplitude, and the first order terms in $`\delta \stackrel{~}{n}`$ are

$$r_1e^{2ikz}\delta \stackrel{~}{n},r_2e^{i2k(z+d)}\delta \stackrel{~}{n}.$$
(43)
To linear order in $`\delta \stackrel{~}{n}`$, the differential reflectivity Eqn. (14) becomes

$`{\displaystyle \frac{\mathrm{\Delta }R}{R}}`$ $`=`$ $`{\displaystyle \frac{8\mathrm{sin}(kd)}{n^21}}[\delta n\mathrm{sin}(2kz+kd)+\kappa \mathrm{cos}(2kz+kd)]`$ (44)
$`=`$ $`{\displaystyle \frac{8\mathrm{sin}(kd)}{n^21}}|\delta \stackrel{~}{n}|\mathrm{sin}(kd)\mathrm{sin}(2kz+\varphi ),`$
where the phase $`\varphi `$ is given by

$$\varphi =kd+\mathrm{arctan}\left(\frac{\kappa }{\delta n}\right).$$
(45)
Note that for small values of $`kd`$, we can expand Eq. (44) in a Taylor series to obtain

$$\frac{\mathrm{\Delta }R}{R}=\frac{8kd}{n^21}\left[\delta n\mathrm{sin}(2kz)+\kappa \mathrm{cos}(2kz)\right]$$
(46)
to first order in $`kd`$.

From Eqn. (44) one can see that the amplitude of the oscillation scales linearly with the change of the real part and imaginary part of the index of refraction respectively.

Since the strain-pulse-front moves with the speed of sound, $`z=C_0t`$, one can rewrite Eqn. (44) as

$$\frac{\mathrm{\Delta }R}{R}(\omega ,t)|\delta \stackrel{~}{n}|\mathrm{sin}(kd)\mathrm{sin}\left(2\pi \frac{t}{T}+\varphi \right)$$
(47)
where $`T=\pi c/(n_\text{b}C_0\omega )=\lambda /(2C_0n_\text{b})`$. Note that in the limit $`kd0`$, Eqn. (47) is equivalent to Eqn. (34) of the microscopic theory described in the last section.

Equation (47) shows that the period of oscillations in the differential reflectivity is given by $`T`$ and hence is determined by the wavelength of the probe pulse. This is consistent with what was seen in these experiments, but differs from previous experiments where the frequency of coherent phonon oscillation was determined by the sample and not the external probe. Note also that Eq. (47) explains the experimental dependence of the the amplitude of the oscillation. It is proportional to the change in the index of refraction and varies with $`\mathrm{sin}(kd)`$, and hence will increase linearly with the thickness of the strained layer as well as the probe frequency $`\omega `$ for $`kd=n_\text{b}(\omega /c)d1`$.

The change of index of refraction in eqn. (47) comes from two effects. The first one, which is small, is the band gap change (renormalization) due to the deformational potential phonon interaction. The other effect results from the change in the built-in piezoelectric field which leads to non-zero absorption below the band gap because of the Franz-Keldysh effect as shown in figure 13.

In general, it is difficult to determine the Franz-Keldysh field experimentally. Here we provide a way to estimate the order-of-magnitude of the built-in piezoelectric field from the strength of the measured differential reflectivity. In Fig. 13 we display the bulk absorption coefficient as a function of photon energy, $`\mathrm{}\omega `$, in bulk GaN (dotted line) at a fixed piezoelectric field, $`F=0.93MV/cm`$ (solid line). The inset shows the absorption coefficient at a fixed photon energy $`\mathrm{}\omega =3.29eV`$ as a function of the piezoelectric field, $`F`$. The change in the absorption at $`\mathrm{}\omega =3.29eV`$ (below the band-gap $`E_g=3.43eV`$), can be used to estimate the piezoelectric field. For a probe energy below band-gap, we read off the amplitude of the differential reflectivity from the experimental data. From this value, and an estimate for $`d`$ based on information about the sample geometry, we obtain the change in the index of refraction which, in turn, gives an estimate for the absorption coefficient $`\alpha (\omega )`$ at the probe photon energy. In the Franz-Keldysh effect, the electro-absorption $`\alpha (\omega )`$ is related to the piezoelectric field, $`F`$, via Haug and Koch (1990)

$$\alpha (\omega )\frac{1}{n_\text{b}}\frac{\omega }{c}\frac{f}{E_g\mathrm{}\omega }\mathrm{exp}\left\{\frac{4}{3f}\left(\frac{(E_g\mathrm{}\omega )}{E_0a_0^2}\right)^{3/2}\right\}.$$
(48)
Here $`f=eF/(E_0a_0^2)`$ where $`E_0=3.435\text{eV}`$ and $`a_0=39.68`$ Å are the excitonic energy and length scales in GaN. From this formula, for instance, we estimate that a piezoelectric field on the order of $`F0.93`$ MeV/cm is responsible for a 1$`\%`$ change in the index of refraction near the band-edge.

The corresponding changes in the real and imaginary parts of the dielectric function are shown in Fig. 14.

Finally, we would like to mention the possibility of using multiple pump pulses to coherently control the generated coherent phonons in the epilayers. There have already been quite a few coherent control experiments with multiple quantum wells from several different groups such as Sun Sun et al. (2001), Özgür Özgür et al. (2001b), and Nelson Feurer et al. (2003). In the InGaN/GaN epilayer systems, the effect of multiple pump pulses should be more prominent then in AlGaAs/GaAs systems because of the stronger piezoelectric field. This might lead to novel applications in imaging or possible new devices based on coherent acoustic phonons.

## V Conclusion

We have presented a theory for the detection of a new class of large amplitude, coherent acoustic phonon wavepackets in femtosecond pump-probe experiments on InxGa1-xN/GaN epilayers and heterostructures. The InGaN/GaN structures are highly strained and at high In concentrations have large built in piezoelectric fields which account for the large amplitude of the observed reflectivity oscillations. This new class of coherent acoustic phonons is generated near the surface and propagates into the structure. The frequency of the reflectivity oscillations is found to be proportional to the frequency of the probe. These coherent phonon wavepackets can be used as a powerful probe of the structure of the sample. We are able to model the generation and propagation of these acoustic phonon wavepackets using a simple string model which is derived from a microscopic model for the photogeneration and propagation of coherent acoustic phonon wavepackets in InGaN/GaN multiple quantum wells. Our model successfully predicts the observed dependence of the coherent phonon reflectivity oscillations on probe wavelength and epilayer thickness.

###### Acknowledgements.

This work was supported by the National Science Foundation through grant DMR 9817828. CSK acknowledges the support from the Korea Research Foundation. We gratefully acknowledge useful discussions with Chi-Kuang Sun.


format: text

text:
# GRB941017: A Case Study of Neutrino Production in Gamma Ray Bursts

## 1 Introduction

Gamma-Ray Bursts (GRBs) are the most powerful objects in the universe, typically emitting luminosities between $`10^{50}`$ and $`\mathrm{\hspace{0.17em}10}^{54}`$ erg/s over 0.1-100 seconds. Their energetics suggest that they may be the sources of the highest energy cosmic rays (Waxman, 1995; Vietri, 1995). It has been pointed out that, if high-energy cosmic rays are accelerated in GRBs, photomeson interactions between accelerated protons and target photons are inevitable, producing very high energy neutrinos and gamma-rays (see, for instance, Waxman & Bahcall 1997). Independent of any association of cosmic rays with GRBs, observable neutrino rates are predicted in models where similar energy goes into the acceleration of electrons and protons in the expanding fireball.

In October of 1994, the Burst And Transient Source Experiment (BATSE) and the Energetic Gamma-Ray Experiment Telescope (EGRET) each observed an exceptionally powerful example of such an object. GRB941017 has the eleventh highest fluence observed in the nine years of BATSE observations. More interesting, however, is the fact that this burst displays not only the typical synchotron-inverse Compton spectrum at $``$30-1000 keV, but also shows a power-law high-energy component extending at least to 200 MeV in energy (González et al. 2003). Other than this burst, EGRET has seen four GRBs at energies of $`100`$ MeV. These were each consistent with an extension of the synchotron-inverse Compton spectrum (Dingus, 2001). Additionally, Milagrito has also observed evidence for emission near $`100`$ GeV in one burst. This observation lacked the ability to reveal significant spectral information, however (Atkins et al., 2000). In a burst-by-burst analysis of the complete BATSE catalogue, Guetta et al. (2003) identified GRB941017 as the most powerful neutrino emitter with, unfortunately, no neutrino telescope to observe it in 1994. Although other explanations may be possible, GRB941017 provides the best evidence to date for high-energy proton interactions with source photons in GRBs.

If protons are accelerated to energies sufficient to produce the gamma-ray spectrum observed in GRB941017, high-energy neutrinos are a necessary consequence (Waxman & Bahcall, 1997). Although no high-energy neutrino telescope of sufficient volume was operational in October 1994, it is interesting to consider the prospects for detection of a burst with similar characteristics in future experiments.

In this letter, we calculate the neutrino spectrum predicted for such an event, and discuss the prospects for detection in future experiments such as the kilometer scale neutrino telescope IceCube (Ahrens et al., 2003). Our conclusions are: i) An event like this is likely to be observable with $`0.55`$ events from a single burst and, ii) a handful of events of this type can produce the diffuse flux of order 10 events per year predicted by fireball phenomenology (Guetta et al., 2003).

## 2 Neutrinos From Photomeson Interactions In GRB Fireballs

Accelerated protons in GRB fireballs produce parent pions via the processes

$`p\gamma \mathrm{\Delta }n\pi ^+`$ (1)
and,

$`p\gamma \mathrm{\Delta }p\pi ^0`$ (2)
which have very large cross sections of $`\sigma _\mathrm{\Delta }5\times 10^{28}\mathrm{cm}^2`$. The charged $`\pi `$’s subsequently decay producing charged leptons and neutrinos, while the neutral $`\pi `$’s decay into high-energy photons. To have sufficient center-of-mass energy for these processes to take place, protons must meet the threshold condition:

$`\epsilon _p^{}{\displaystyle \frac{m_\mathrm{\Delta }^2m_p^2}{4\epsilon _\gamma ^{}}}.`$ (3)
Primed and unprimed quantities refer to values measured in the comoving and observer’s frames, respectively. In the observer’s frame,

$$\epsilon _p1.4\times 10^{16}\frac{\mathrm{\Gamma }_{2.5}^2}{\epsilon _{\gamma ,\mathrm{MeV}}}\mathrm{eV},$$
(4)
resulting in a neutrino energy

$$\epsilon _\nu =\frac{1}{4}x_{p\pi }\epsilon _p7\times 10^{14}\frac{\mathrm{\Gamma }_{2.5}^2}{\epsilon _{\gamma ,\mathrm{MeV}}}\mathrm{eV},$$
(5)
where $`\mathrm{\Gamma }_{2.5}=\mathrm{\Gamma }/10^{2.5}`$ is the bulk Lorentz factor and $`\epsilon _{\gamma ,\mathrm{MeV}}=\epsilon _\gamma /1\mathrm{M}\mathrm{e}\mathrm{V}`$ is the typical target photon energy. $`x_{p\pi }0.2`$ is the average fraction of energy transferred from the initial proton to the produced pion. The factor of 1/4 comes from the assumption that the 4 final state leptons in the decay chain $`\pi ^+\nu _\mu \mu ^+\nu _\mu e^+\nu _e\overline{\nu _\mu }`$ equally share the pion energy.

Typical GRBs display a broken power-law spectrum consistent with synchotron and inverse Compton emission. This spectrum can be parameterized by:

$$F_\gamma =\epsilon _\gamma dn_\gamma /d\epsilon _\gamma \{\begin{array}{cc}\epsilon _\gamma ^\alpha ;\hfill & \epsilon _\gamma <\epsilon _\gamma ^b\hfill \\ \epsilon _\gamma ^\beta ;\hfill & \epsilon _\gamma >\epsilon _\gamma ^b\hfill \end{array}.$$
(6)
In the case of GRB941017, $`\alpha `$ and $`\beta `$ appear to change with time varying between 1.2 and 2.1 for $`\beta `$ and between $`0.2`$ and $`0.5`$ for $`\alpha `$. In the first tens of seconds of the burst, the break energy is about 500 keV. It steadily drops, falling well below 100 keV after about 100 seconds.

For each proton energy, the resulting neutrino spectrum traces the GRB photon spectrum, but with a much higher energy break:

$$\epsilon _\nu ^b=7\times 10^{14}\frac{1}{(1+z)^2}\frac{\mathrm{\Gamma }_{2.5}^2}{\epsilon _{\gamma ,\mathrm{MeV}}^b}\mathrm{eV}.$$
(7)
$`z`$ is the redshift of the GRB. For this particular burst, the redshift has been estimated to be $`z0.6`$, although its precise value is uncertain (Guetta et al., 2003).

At very high energies, pions lose energy by synchrotron emission before decaying. This affects the spectrum of neutrinos from pion decay above the energy:

$$\epsilon _{\nu _\mu }^s=\frac{10^{17}}{1+z}L_{\gamma ,52}^{1/2}\mathrm{\Gamma }_{2.5}^4t_{v,2}\mathrm{eV},$$
(8)
where $`L_{\gamma ,52}`$ is the gamma-ray luminosity of the burst in units of $`10^{52}`$ erg/s. For GRB941017, $`L_{\gamma ,52}1`$. $`t_{v,2}`$ is the time scale of fluctuations in the GRB lightcurve in units of $`10^2\mathrm{s}`$ Neutrinos from muon decay have an energy cutoff that is 10 times smaller due to the longer lifetime of the muon compared to that of the pion. Above this energy, the slope of the neutrino spectrum steepens by two.

The fraction of energy in accelerated protons converted to pions in the fireball is estimated from the ratio of the size of the shock, $`\mathrm{\Delta }R^{}`$, and the proton mean free path:

$`f_\pi {\displaystyle \frac{\mathrm{\Delta }R^{}}{\lambda _{p\gamma }}}x_{p\pi }.`$ (9)
Here, the proton mean free path is given by $`\lambda _{p\gamma }=1/n_\gamma \sigma _\mathrm{\Delta },`$ where $`n_\gamma `$ is the number density of photons. The photon number density is given by the ratio of the photon energy density and the photon energy in the comoving frame:

$$n_\gamma =\frac{U_\gamma ^{}}{\epsilon _\gamma ^{}}\left(\frac{L_\gamma t_v/\mathrm{\Gamma }}{4\pi R^2\mathrm{\Delta }R^{}}\right)/\left(\frac{\epsilon _\gamma }{\mathrm{\Gamma }}\right).$$
(10)
Using these equations, and the relationship $`R2\mathrm{\Gamma }^2ct_v`$, it is found

$`n_\gamma \left({\displaystyle \frac{L_\gamma }{16\pi c^2t_v\mathrm{\Gamma }^5\mathrm{\Delta }R^{}}}\right)/\left({\displaystyle \frac{\epsilon _\gamma }{\mathrm{\Gamma }}}\right)={\displaystyle \frac{L_\gamma }{16\pi c^2t_v\mathrm{\Gamma }^4\mathrm{\Delta }R^{}\epsilon _\gamma }},`$ (11)
which leads to

$$f_\pi \frac{L_\gamma }{\epsilon _\gamma }\frac{1}{\mathrm{\Gamma }^4t_v}\frac{\sigma _\mathrm{\Delta }x_{p\pi }}{16\pi c^2}0.2\times \frac{L_{\gamma ,52}}{\mathrm{\Gamma }_{2.5}^4t_{v,2}\epsilon _{\gamma ,\mathrm{MeV}}^b}.$$
(12)
Thus far, this calculation describes protons at the break energy. In general,

$$f_\pi (\epsilon _p)0.2\frac{L_{\gamma ,52}}{\mathrm{\Gamma }_{2.5}^4t_{v,2}\epsilon _{\gamma ,\mathrm{MeV}}^b}\times \{\begin{array}{cc}(\epsilon _p/\epsilon _p^b)^\alpha \hfill & \epsilon _p>\epsilon _p^b\hfill \\ (\epsilon _p/\epsilon _p^b)^\beta \hfill & \epsilon _p<\epsilon _p^b\hfill \end{array},$$
(13)
where $`\epsilon _p^b`$ is given by Eq.(4).

This equation indicates that $`f_\pi `$ and, therefore, the spectrum’s normalization and the final event rates, vary significantly from burst-to-burst (Halzen & Hooper, 1999; Alvarez-Muñiz et al., 2000). Such fluctuations are constrained, however, and similar conclusions are reached when fixing $`f_\pi =0.2`$, its typical value (Guetta et al., 2001).

To estimate the bulk Lorentz factor, $`\mathrm{\Gamma }`$, we relate it to the peak energy of the gamma-ray spectrum:

$$\epsilon _\gamma ^b\frac{L_{\gamma ,52}^{1/2}}{\mathrm{\Gamma }_{2.5}^2t_{v,2}}\mathrm{MeV}.$$
(14)
For the case of GRB941017, we estimate $`\mathrm{\Gamma }`$ 250 and this is again somewhat uncertain (Guetta et al., 2003).

Finally, we obtain the neutrino spectrum:

$`\epsilon _\nu ^2{\displaystyle \frac{dN_\nu }{d\epsilon _\nu }}{\displaystyle \frac{1}{8}}{\displaystyle \frac{F_\gamma }{\mathrm{ln}(10)}}f_\pi `$ (15)
where, $`F_\gamma `$ and $`f_\pi `$ are determined by equations 6 and 13, respectively. The neutrino spectrum found for GRB941017 in the time bins in (González et al. 2003), is shown in figure 1.

## 3 Event Rates in Neutrino Telescopes

Currently, the AMANDA-II detector (Andres et al., 2001), located at the South Pole, provides the strongest limits for high-energy neutrinos from GRBs (Barouch & Hartdke, 2001; Stamatikos et al., 2003). Considerably larger experiments such as the kilometer-scale neutrino observatory IceCube (Ahrens et al., 2003), are presently under construction. For a review of high-energy neutrino astronomy, see Halzen & Hooper 2002 or Learned & Mannheim 2000.

High-energy neutrinos are observed as muons or showers created in interactions inside or near a Cherenkov detector embedded in an optically transparent medium such as ice or water. Muons, created in charged current interactions travel several kilometers before losing the majority of their energy, thus enhancing their prospects for observation. The probability for observing an energetic muon is given by

$$P_{\nu _\mu \mu }=N_\mathrm{A}R_\mu \sigma _{\mathrm{CC}},$$
(16)
where $`N_\mathrm{A}`$ is Avogadro’s number, $`R_\mu `$ is the muon range (Dutta et al., 2001) and $`\sigma _{\mathrm{CC}}`$ is the neutrino-nucleon charged current cross section. For GRB941017, with a zenith angle of 94.5 (for experiments at the South Pole), there is sufficient ice to allow for very long muon ranges. In our calculations, we take into account the absorption of neutrinos in the Earth as well as the effect of oscillations.

All three flavors of neutrinos may interact producing showers within the detector volume. The probability of making such an observation is similar to the relation for muons, although shower events do not benefit from long muon ranges, and typically they have to be produced inside the detector to be observed. More details on how the rate is calculated can be found, for instance, in Appendix C of Guetta et al. 2003.

The rates for a GRB similar to GRB941017 are shown in table 1. We have made two estimates of the rate. The first estimate uses Eq. 15 for the neutrino flux, obtaining $`F_\gamma `$ from a broken power-law fit to the observed gamma-ray spectrum. It is important to remark that this fit cannot accommodate the late-time high-energy feature of the gamma-ray spectrum. We predict 0.3 events per square kilometer, or roughly 1 event in IceCube with an effective area for muons and showers exceeding 1 km 2. Bear in mind that for an observation over a narrow temporal and angular window in coincidence with the optical display, the background is entirely negligible (Guetta et al., 2003). A single high-energy neutrino event is a conclusive observation. This rate is, of course, subject to the ambiguities associated with $`z`$ and $`\mathrm{\Gamma }`$.

To accommodate the late-time, high-energy feature seen in the gamma-ray spectrum, instead of using the method described above, we can, alternatively, estimate the neutrino flux and event rate by relating the observed high-energy gamma-ray component of GRB941017 to a corresponding high-energy neutrino component. This method assumes that the observed high-energy gamma-ray component is the result of photomeson interactions. The observed high-energy gamma-ray component of GRB941017 can be parameterized (González et al. 2003) by

$$\frac{dN_\gamma }{dE_\gamma }A_\gamma E_\gamma ^\delta $$
(17)
where $`\delta `$ is measured to be 0.96 to 1.10 (González et al. 2003) and $`A_\gamma `$ is a normalization constant fixed to the observed spectrum. The corresponding neutrino flux can be calculated by energy conservation given that the physics of photomeson interactions fixes the ratio of energy going to photons and neutrinos (Alvarez-Muñiz & Halzen, 2003). The main uncertainties are the slope of the neutrino spectrum, and the energy to which the high-energy gamma-ray component extends. Without cascading of the gamma rays in the source both would be uniquely determined by the photon spectrum. This is not to be expected and we therefore calculate the neutrino flux varying the neutrino spectral slope and the high-energy cutoff of the gamma ray spectrum. Note that the latter is unobservable because high-energy gamma rays are absorbed on diffuse infrared background radiation.

The results are shown in figures 2 and 3. For a $`E^1`$ neutrino spectrum and a gamma-ray cutoff near 100 TeV, we predict $`10`$ events per square kilometer. Experiments capable of observing TeV gamma-rays from GRBs, such as MILAGRO (Atkins et al., 2001) and IceCube (Halzen & Hooper, 2003) will be useful in determining the maximum energy to which the photomeson component of the gamma-ray spectrum extends.

Another way of summarizing our results is shown in figure 3 where the rate is plotted as a function of the neutrino spectral slope. The observed event rate increases with the spectral index because more neutrinos of lower energy are produced and, although their detection probability is smaller, their larger number compensates. This trend continues until a neutrino spectral slope of $`2.1`$. Above this slope, a significant amount of the energy goes into neutrinos below the experimental energy threshold of the detector and the rate decreases. The same argument explains why for a flatter neutrino spectrum, similar to the observed photon spectrum, the same event rate is reached for a higher gamma-ray cutoff as seen in figure 2.

Although $`10^3`$ GRBs are observed per year, the total rate from all of these events is typically estimated to be $`10`$ events per square kilometer per year (Guetta et al., 2003). A small number of exceptional bursts, such as the one we study here, can contribute substantially to the total neutrino flux.

## 4 Conclusions

We have pointed out that the high-energy feature observed in GRB941017 provides further support for the expectation of detectable fluxes of high-energy neutrinos in coincidence with gamma ray bursts. We have estimated the neutrino spectrum that would accompany such a burst and discussed the prospects for such an event’s detection in future neutrino observatories. We find that such an event is expected to produce on the order of 1 event in a kilometer scale neutrino telescope and that this would be a conclusive observation since there is no competing background during the time and in the direction of the burst.

J. Alvarez-Muñiz is supported in part by Xunta de Galicia (PGIDT02 PXIC 20611PN) and by MCyT (FPA 2001-3837 and FPA 2002-01161). F. Halzen work is supported in part by DOE grant No. DE-FG02-95ER40896, NSF grant No. OPP-0236449 and in part by the Wisconsin Alumni Research Foundation. D. Hooper is supported by the Leverhulme Trust.


text:
# Two Combinatorial Models with identical Statics yet different Dynamics

## 1 Introduction

Viewing optimization problems that arise in computer science from the perspective of statistical mechanics has led to successful insights . From this point of view, where optimization algorithms are intimately related to the energy or cost function, the features of the energy landscape are crucial in determining the success or otherwise of optimization. This connection is less evident in computer science since algorithms do not necessarily have such a direct connection with the cost function. One of the open questions of the field is to understand the extent to which the performance of algorithms can be determined by the character of the optimization problem and not the details of the algorithm itself .

In this paper we hope to illuminate the issue by introducing a pair of models with the unusual feature that although they have the same static properties determined by the list of energy levels, the energy landscape is different since the notion of which states are close to each other is not the same. We perform dynamical simulations of relaxation after a quench to explore the energy landscape for each of the models following a similar strategy to that used to investigate other computer science problems such as satisfiability (K-SAT) and graph colouring .

One of the best known optimization models is the Traveling Salesman Problem (TSP) . States can be labeled by permutations since the paths in the TSP correspond to various orders of the cities visited. The combinatorial models we consider in this paper have states which are permutations of individual numbers rather than permutations of more general objects such as the set of city coordinates used in the TSP. In this case the natural optimization problem is sorting the numbers and we shall take that problem as sufficient motivation. The problem of sorting has been extensively treated in the computer science literature  and many efficient algorithms are known. Most of this work is concerned with analyzing the time taken to perform the sort, though issues such a memory requirement and ability to use cache are also important.

For some optimization problems, such as TSP, the cost function is clearly the length of the path and this can be taken as the Hamiltonian of the statistical mechanical model. The problem of sorting does not have an uniquely obvious cost function. We require that the lowest energy corresponds to the sorted state but have freedom to decide how to measure the degree other permutations are sorted. In this paper we consider two different energy functions that compute the degree of sortedness in different ways: one is similar to the TSP and measures the cost in terms of the length of a path, the other assumes knowledge of the sorted state and computes the distance from that state in a direct manner.

By investigating these two models we discover that the spectrum of energies is identical. This implies that the partition functions are also identical and all static properties will be the same. Yet the Hamiltonians are different, and the energy each Hamiltonian assigns to a given state is not the same. We make some preliminary observations on the mapping between states with a given energy according to one Hamiltonian and the same energy according to the other Hamiltonian, but do not delve into mathematical details. Of more concern to this presentation is the fact that although no physical distinction between the two models is visible at the static level, it is manifest in the dynamics. We investigate the local dynamics after a deep quench and show how it displays very different behavior for each model: in one case the (sorted) ground state is eventually found, in the other case it is not. The difference in behavior is identified as being due to rather different energy landscapes, with stable local minima in one Hamiltonian but not the other.

From a physical point of view, study of the dynamical relaxation after a quench falls under the topic of phase ordering kinetics . From the perspective of computer science, it is the natural way to study the efficiency of a local search algorithm in finding the optimum solution to a problem.

## 2 Models

### 2.1 State Space

For statistical mechanical approaches to sorting, the states or configurations are the permutations of a set of numbers. In the following, we shall consider the set of $`N`$ integers $`1,2\mathrm{}N`$. A more general model, in which the numbers are taken to have randomly chosen real values, is briefly described in appendix A. Most of the properties investigated in this paper hold for both models, but for the sake of clarity we shall work exclusively with the simpler integer model.

The states correspond to all possible permutations $`P`$ of the set of $`N`$ integers. The permutation label can be regarded either passively as an ordered $`N`$-tuple or actively as the mapping needed to get to that $`N`$-tuple from the identity . Generally we shall employ the first interpretation, so $`P_i`$ signifies the $`i`$’th element of the permutation, ($`P_1,P_2\mathrm{}P_N`$), but the second interpretation will be convenient later when we use it to write the permutation in terms of cycles. We imagine the $`P_i`$ as analogs of spins on a line, and sometimes we will refer to this line as a one dimensional spatial direction.

The size of the state space grows as $`N!`$ in comparison with typical (Ising) spin models where the space grows as $`2^N`$. It is well-known that models such as this, where the state space grows faster than exponentially, have difficulties of interpretation since scaling of the temperature or energy with $`N`$ is necessary to ensure that certain quantities are extensive. In this work, we do not investigate the phase structure, and never need to perform this scaling since we only consider zero (or infinite) temperature. We refer the reader to Mézard and Parisi  and to Anderson and Fu  for further discussion of this matter.

### 2.2 Two Energy Functions

We consider two energy or cost functions that introduce a measure of the distance of a permutation away from the identity permutation. The energy is lowest and vanishes when evaluated for the ordered or identity permutation. In writing these expressions it is convenient to introduce the analog of the Kronecker delta overlap between individual spins as the one dimensional distance metric between the numbers $`P_i`$:

$$d(P_i,P_j)=|P_iP_j|$$
(1)
The first energy is familiar as the cost function for the TSP , here evaluated for the one-dimensional case. The boundary conditions are slightly different from the usual TSP since we insist they be fixed rather than periodic. The connection with sorting is clearer with fixed boundary conditions and we need not consider the degeneracy of all states under cyclic shifts or inversions. However, it should be noted that the bulk part of the energy minimizes to either ascending or descending order, and it is only the boundary conditions that select ascending order. This energy, which henceforth we shall call the TSP energy ($`E_{TSP}`$), is:

$`E_{TSP}(P)`$ $`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=0}{\overset{N}{}}}d(P_{i+1},P_i){\displaystyle \frac{(N+1)}{2N}}`$ (2)
$`=`$ $`{\displaystyle \frac{1}{2N}}\left(P_1P_N\right)+{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}|P_{i+1}P_i|`$
The first form is written in the standard form for TSP and we have additionally assumed: $`P_0=0`$, $`P_{N+1}=N+1`$ for any $`P`$. In the second form the constant term is removed, by writing the end terms of the sum explicitly. Instead of using the function (1), powers (notably quadratic) of the node position differences could be considered, but these energy functions do not appear to be natural in this context and do not obey the properties we will demonstrate below.

The energy function for our second model, which we shall term the displacement energy due to its connection with the “total displacement” of R.W. Floyd, is given as a sum of site distances from the sorted configuration $`P_i^{(0)}=i`$.

$$E_D(P)=\frac{1}{2N}\underset{i=1}{\overset{N}{}}d(P_i^{(0)},P_i)=\frac{1}{2N}\underset{i=1}{\overset{N}{}}|P_ii|$$
(3)
This form of cost function is familiar as an “assignment problem” .

The choice of distance measure relying on site differences is by no means unique. Many other ways of defining the distance between configurations are possible. For example, in the physics literature, an overlap based on counting similar links is common . Another overlap (that happens to yield a tractable model ) is based on the matching problem and is determined by counting the number of sites that are in their correct relative positions.

The different cost functions can be interpreted geometrically as illustrated in figure 1. For the permutation, 1364725 with $`N=7`$, used as the example in the figure, the energies may be computed as:

$`E_{TSP}`$ $`=`$ $`{\displaystyle \frac{1}{2N}}\left(6+2+3+5+6\right){\displaystyle \frac{1}{2N}}\left(8\right)=1`$
$`E_D`$ $`=`$ $`{\displaystyle \frac{1}{2N}}\left(1+3+4+2+2\right)=6/7`$ (4)
These expressions and the diagrams in the figure make it clear that the energy may be computed as a sum over runs (either ascending or descending) in the case of the TSP energy, and as a sum over cycles in the case of the displacement energy. A run is the term used in the combinatorial literature for a subsequence of adjacent elements that are in sorted (or antisorted) order. For the TSP energy each contribution is the difference between the maximum and minimum value contained in the run. However, the contribution of each cycle to the displacement energy is not usually a single term (for more complex permutations than that shown in the figure) and it is necessary to look at runs within a cycle.

## 3 Relation between Models

The energy spectra of $`E_{TSP}`$ and $`E_D`$ are identical. That is, there is a one-one map between the full list of energies for all states computed with $`E_D`$ and the list computed with $`E_{TSP}`$. Of course there is a shuffling in the way the states are associated with energies.

For example, figure 2 shows the spectrum of the $`N=4`$ model, with states labeled for each energy function. For many, but not all of these $`N=4`$ states, the energy is the same whether evaluated with the TSP or the displacement energy. The proportion of states with invariant energy decreases at larger $`N`$. There are many degeneracies between energy levels: in the figure there are a 9-plet, a 7-plet, a 4-plet and a triplet besides the singlet ground state (which is always unique) that together make the $`4!=24`$ states. In appendix B we list the multiplicity structure for small values of $`N`$. This information might be expected to lead to an expression for the partition function. However, we have been unable to obtain a general formula for these multiplicities and according to Knuth, the generating function does not appear to have a simple form.

### 3.1 Mapping between Models

To demonstrate the relation between the two models formally, we present a mapping between the states $`PP^{}`$, (in fact a permutation in the state space) that has the property that $`E_D(P)=E_{TSP}(P^{})`$.

The mapping relies on the representation of the state $`P`$ in terms of a permutation of the ordered state $`(1,2,3\mathrm{}N)`$ using the cyclic representation. This mapping is well-known and is for example described in Knuth .

We write any permutation in terms of $`M`$ distinct cycles (singleton cycles included explicitly). Each cycle is labeled by a superscript $`j`$ and contains $`m_j`$ elements in the cycle as:

$$(i_1^1i_2^1i_3^1\mathrm{}i_{m_1}^1)(i_1^2i_2^2i_3^2\mathrm{}i_{m_2}^2)\mathrm{}(i_1^ji_2^ji_3^j\mathrm{}i_{m_j}^j)\mathrm{}(i_1^Mi_2^Mi_3^M\mathrm{}i_{m_M}^M)$$
(5)
We fix the freedom available in the way the cycle form is written by requiring that:

* The largest element appears at the start of each cycle: $`i_1^j>i_k^j`$ for all other $`k`$ in the cycle $`k=2,3\mathrm{}m_j`$.
* Cycles appear in the order determined by their first elements: $`i_1^{j+1}>i_1^j`$ for $`j=1,2\mathrm{}(M1)`$.

Note that this is slightly different from the convention defined by Knuth.

The action of this permutation on the ordered state $`(1,2,3\mathrm{}N)`$ is to take the value $`i_k^j`$ to $`i_{k+1}^j`$ (or cyclically). The resulting permutation is the original state $`P`$ of the map. The state $`P^{}`$ obtained from the map is given simply by removing brackets from the definition of the element (5) above and regarding the list of numbers as a permutation.

For example with $`N=4`$, the cyclic form $`(1)(423)`$ takes $`1234`$ to $`1342`$. Thus $`P`$ = $`1342`$ and $`P^{}`$ = $`1423`$. In this case, both $`P`$ and $`P^{}`$ lie in the same multiplet according to either the TSP or displacement energy. This is not the case for $`P`$ = $`3412`$ and $`P^{}`$ = $`3142`$ that appear at the top of figure 2 and correspond to the cyclic form $`(31)(42)`$. Nonetheless, $`E_D(3412)`$ = $`E_{TSP}(3142)`$.

The displacement energy in the state $`P`$ is given by a sum over cycles as:

$$E_D=\frac{1}{2N}\underset{j=1}{\overset{M}{}}\left(d(i_1^j,i_2^j)+d(i_2^j,i_3^j)+\mathrm{}+d(i_{m_{j1}}^j,i_{m_j}^j)+d(i_{m_j}^j,i_1^j)\right)$$
(6)
The TSP energy can be evaluated in the state $`P^{}`$ to obtain three terms, one from the elements at the end of the original sum in (2), one from the contribution of each cycle, and a term corresponding to the intercycle contributions.

$`E_{TSP}`$ $`=`$ $`{\displaystyle \frac{1}{2N}}(i_1^1i_{m_M}^M)+{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{M1}{}}}d(i_1^{j+1},i_{m_j}^j)`$ (7)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{M}{}}}\left(d(i_1^j,i_2^j)+d(i_2^j,i_3^j)+\mathrm{}+d(i_{m_{j1}}^j,i_{m_j}^j)\right)`$
The cycle sums appearing in (6) and the last term in (7) are identical except for one additional term in the first equation. The total difference between the energies is given by:

$$E_{TSP}E_D=\frac{1}{2N}(i_1^1i_{m_M}^M)+\frac{1}{2}\underset{j=1}{\overset{M1}{}}|i_1^{j+1}i_{m_j}^j|\frac{1}{2}\underset{j=1}{\overset{M}{}}|i_1^ji_{m_j}^j|$$
(8)
Now, by using the requirements on ordering of the $`i_k^j`$’s stated above, we find that the modulus signs may be removed in each sum and the total vanishes. This does not occur for other choices of distance function, but does continue to hold for the model based on real numbers rather than integers.

The map is invertible due to the requirements listed above. Note that the map must leave the ground state unchanged. However, the multiplet structure is not preserved, states that appear in a certain multiplet according to one energy function may appear in a different multiplet according to the other energy function. This can be seen in the case of $`N=4`$ as the states marked with a “$``$” in figure 2.

Several distinct versions of the map exist, based on different conventions for the way the cycle form is written. An approach to understanding the symmetries of the system would be to combine a map and the inverse of a different version. We do not consider this approach here since it takes us too far afield from the aim of the present paper.

### 3.2 Average Energy

Here we consider statistics of the energy (either TSP or displacement) spectrum: the average energy and its fluctuation. These are averages from a combinatorial point of view in which all states contribute equally; thermodynamically they are effectively at infinite temperature.

An estimate of the large $`N`$ behavior is most easily obtained by writing the TSP energy in terms of $`\mathrm{\Delta }`$, the average distance between nodes:

$$E_{TSP}(N)=\frac{1}{2N}\underset{i=0}{\overset{N}{}}|P_{i+1}P_i|\frac{1}{2}\mathrm{\Delta }=\frac{N}{6}$$
(9)
Where the approximation consists in ignoring correlations between the distances between different pairs of nodes. In this case, the pairs can be imagined as two independent points thrown at random in the interval $`[0,N]`$, so the average distance between them is $`\mathrm{\Delta }=N/3`$.

A more formal computation based on the displacement energy proceeds as follows, where $`_P`$ indicates a sum over all permutations.

$`E(N)`$ $`=`$ $`{\displaystyle \frac{1}{N!}}{\displaystyle \underset{P}{}}{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{N}{}}}|P_ii|`$ (10)
$`=`$ $`{\displaystyle \frac{1}{2NN!}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{P}{}}|P_ii|`$
$`=`$ $`{\displaystyle \frac{(N1)!}{2NN!}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{k=1}{\overset{N}{}}}|ki|`$
The order of summation is exchanged in the second line after which the sum over permutations becomes simple since there are $`(N1)!`$ permutations in which $`P_i`$ takes a given value $`k`$.

The resulting sum can be performed using standard techniques and a similar argument holds for the second moment.

$`E(N)`$ $`=`$ $`{\displaystyle \frac{N^21}{6N}}`$
$`(E(N)E)^2`$ $`=`$ $`{\displaystyle \frac{(N+1)(2N^2+7)}{180N^2}}`$ (11)
The width of the energy distribution, $`\sqrt{(EE)^2}`$, scales as $`N^{1/2}`$ and therefore becomes relatively more peaked at large $`N`$.

## 4 Dynamics

Although the energy spectra of the two models is identical and the partition functions are the same, the models still have distinct properties. This would be apparent by studying the response to some external field that couples to states in the same way in each model. The displacement energy itself could be regarded as an example of this kind of additional term in the cost function. However, in the context of this paper and the problem of sorting, the distinction is best studied by looking at the dynamical properties of the models.

We only consider local dynamics: that is the basic moves are adjacent transpositions. Of course there is no physical basis to these models requiring locality, and in real sorting algorithms non-locality of the elementary moves is essential to obtain efficient sorting. Furthermore, we only consider a dynamics associated with the energy function - namely the Monte Carlo Metropolis algorithm. This is natural from a theoretical physics perspective (even though it does not correspond to any physical dynamics ), but algorithms that have much less clear relationship with the cost function are common in computer science. Indeed we could imagine a reasonable algorithm that selects random sites and transposes with the neighbor if they are out of order.

The usual approach physicists have used to investigate computer science problems is to study the dynamical relaxation under local search algorithms . From the physical point of view, the approach consists in studying the phase ordering kinetics after a deep quench .

### 4.1 Zero Temperature Metropolis Algorithm

The Monte Carlo moves are local adjacent transpositions and at zero temperature this effectively constitutes a randomized steepest descent algorithm. The Metropolis algorithm we shall use selects a trial site at random and transposes with its (right hand) neighbor provided this move reduces the energy (or does not change it). We perform numerical simulations starting from a random configuration which typically has energy very close to the average energy computed in section 3.2.

Figure 3 shows the results of numerical studies of dynamics according to this algorithm. The two curves in the figure correspond to dynamics based on each of the energy functions we have defined. Starting from a randomly chosen initial state, the plot follows the evolution of the energy averaged over many choices of this initial state. This situation corresponds to a quench from very high temperature to zero temperature.

The model based on the displacement energy evolves in the expected manner: the energy slowly reduces and eventually reaches the sorted ground state. The other model, based on the TSP energy starts with a similar initial energy then rapidly reduces to a plateau value at which level it continues indefinitely. No further decrease in energy is evident, even for much longer runs, and the energy never arrives at the sorted ground state. This behavior can be improved somewhat by using simulated annealing but it remains extremely slow and still tends to get stuck.

### 4.2 Timescales for Energy Decay

According to figure 3, the displacement energy appears to reduce at an exponential rate. The log-log plot (not shown) is only able to substantiate this for the first part of the decay. The quality of this initial exponential decay, and the form of the subsequent deviations from this form are shown in figure 4.

Within the initial exponential decay region it is possible to measure the $`N`$ dependence of the exponential timescale. At short times this characteristic time scales as $`N^2`$ (fits using times up to order $`N^2`$ give the exponent with accuracy $`0.1\%`$, and value tending to 2.0 as fewer points are taken), so the time axis of the figure has been scaled as $`t/N^2`$ in order to collapse the data in the initial region.

The deviation from exponential form and absence of data collapse in the later part of the data is a finite size effect. With the help of small quantities of data for very large sizes up to $`N=10^5`$, careful measurements of the time required to reach fixed values of $`E/N`$ can be made. The way these times scale with $`N`$ shows a consistent trend towards an exponent of 2.0 at larger $`N`$.

Certainly for any size that can be simulated, the total time to arrive at zero energy is affected by the finite size effects and has a scaling exponent larger than 2. Not surprisingly, the resulting sort is rather slower than achieved with standard sorting algorithms that have best case behavior increasing as $`N\mathrm{log}N`$.

For the TSP energy, figure 5 shows a scaled plot for different size systems. In this case the time is only scaled by a factor of $`1/N`$ and the energy axis is not logarithmic. The stability of the plateau energy is very clear in this figure. A detailed investigation finds that after the effect of finite size effects are removed, the per site energy on the plateau is $`E_{TSP}/N=0.1092\pm 0.0001`$.

The fact that the evolution timescales of the TSP and displacement models are respectively $`N`$ and $`N^2`$ indicates another surprising distinction between the dynamics of the two models.

### 4.3 Spatial Correlations

Throughout this work we have implicitly considered the indices to form a line. In this section we consider correlations along this line and refer to such correlations as spatial. Studying the evolution of these correlators may illuminate the dynamical properties of the model since dynamical lengthscales may appear. By analogy with the definition of the distance metric (1), we define the spatial correlator at distance $`r`$ as:

$$C(r)=\frac{1}{N(Nr)}\underset{i=1}{\overset{Nr}{}}|P_{i+r}P_i|$$
(12)
Up to terms relating to boundary conditions, $`C(1)`$ is nothing other than the TSP energy. A small value of $`C`$ signifies a strong correlation, and in the fully sorted state it grows linearly $`C(r)=r/N`$. From arguments similar to those used to derive the average energy, it can be shown that the correlator takes a constant value of $`1/3`$ when averaged over randomly chosen configurations (more precisely, $`C(r)=(N+1)/3N;r>0`$). In the figures to follow, we only show $`C(r)`$ up to values of $`r<N/2`$. This is because for larger $`r`$, only a small number of pairs appear in the sum (12) and the indices of these pairs are often near the boundaries, thus making this region excessively dependent on boundary effects.

As the configuration evolves according to the TSP model, the spatial correlator remains similar to its initial constant value $`1/3`$ corresponding to the initial random configuration, but some structure develops at small $`r`$. Once the plateau is reached, there is no further change in the spatial correlator and figure 6 shows its final form. The non-trivial structure has a clear size, of less than 10 units, and neither the size scale nor the shape of the structure depend on $`N`$. This scale indicates the distance over which ordering can take place according to the TSP dynamical process. The mechanism that blocks further growth of the order beyond this scale is discussed below.

In the case of the displacement energy, the correlator evolves as shown in figure 7. Here, no detailed structure ever appears, and the correlator is always a straight line with a gradient that smoothly evolves towards the steepest slope corresponding to the final sorted ground state. There is no characteristic lengthscale over which ordering takes place and then grows. It would rather appear that the system becomes organized on all scales simultaneously.

### 4.4 Time Correlations

In order to investigate correlations between configurations at different times, we use the same overlap that was employed for the displacement energy.

$$C(t,t^{})=\frac{1}{2N}\underset{i=1}{\overset{N}{}}|P(t)_iP(t^{})_i|$$
(13)
Other possibilities are certainly possible. For example, an overlap based on matching counts the number of positions where the number has not changed. The graphs based on this choice do not convey any significantly different information from those given below.

The simplest correlation to measure is against the completely sorted ground state. This is none other than the displacement energy itself, and was shown (for dynamics based on this energy) in figure 3.

Another familiar comparison is against the initial state. For the dynamics according to the displacement energy, the correlation with the initial state disappears rapidly (within about $`20N`$ steps), and we do not show a figure in this case. For the dynamics according to the TSP energy, we show in figure 8 the two-time correlators for different waiting time plotted against $`\mathrm{log}(tt_w)`$ in the conventional manner. For $`t_w=0`$ the correlation against the initial state is included in this figure. We have drawn two-time correlators not because aging occurs in this model - it does not, but as a convenient way of demonstrating two features: that the configurations on the plateau retain some correlation with the original random state, and to show that evolution is not frozen on the plateau. With this aim, the waiting times shown in figure 8 are quite short.

The final value of the $`t_w=0`$ curve (about $`0.092`$) is much less than that associated with correlation between random configurations ($`1/6`$ in this case due to the factor of 2 in the definition (13)), indicating that not all information in the original configuration has been lost by the time the plateau is reached. This agrees with the result of the spatial correlation that indicates that only local modification takes place.

The lowest two-time curve holding for all $`t_w>1000`$ (for size $`N=100`$) corresponds to waiting times that have reached the plateau. This curve is not constantly zero as would be the case if there was no dynamics happening on the plateau. The configuration continues to evolve, though the energy does not change. However, the correlation is bounded: irrespective of how long after the waiting time, the correlator never rises above a certain value (about $`0.021`$). This limiting value is independent of $`N`$. The reason for this behaviour is the limited size of the flat directions that are identified below.

### 4.5 Local Energy Minima

The fact that a simple algorithm such as zero temperature Monte Carlo is unable to find the ground state is hardly a surprise in optimization models. There are many statistical mechanical examples of this state of affairs, and the effect is usually ascribed to the features of the free energy landscape. In the K-Satisfiability problem the difficulty of finding the ground state via a local search procedure is due to the proliferation of states which trap the search into metastable phase. Eventually, a large fluctuation provides a means of reaching the ground state . However, in our TSP model the reason is more prosaic and the origin of the effect is due to the presence of stable local energy minima. In this respect, it is similar to the XOR-SAT model, that also suffers from such minima .

Though the two energy functions have matching energy levels they have very different energy landscape characteristics. The displacement version has no stable local minima, but the TSP version does. Moreover, these minima appear to proliferate as $`N`$ increases and are never avoided.

For small values of $`N`$ the local minima of the TSP model may be found explicitly. No such minima exist for $`N=3`$. For $`N=4`$ there are four minima connected by Monte Carlo moves as shown in figure 9. Note that these minima are precisely the states that appear with a “$``$” in figure 2, and indeed the energy assignments of figure 9 are inverted for the displacement energy, so in that case there are no local minima.

For larger $`N`$, the trapping configurations appear as shown in figure 10. They consist of alternating ascending and descending runs separated by turning pairs. For an upper turning pair, each element of the pair is greater than either of the elements neighbouring the pair, with a similar definition for a lower turning pair. These trapping configurations found on the plateau, are slightly more ordered than random configurations that consist of ascending and descending runs separated by ordinary turning points. A turning point is a requirement on a subsequence of length 3, whereas a tuning pair is a requirement on a longer subsequence of length 4. The dynamical Monte Carlo process performs this small scale ordering (as observed in the spatial correlators) to arrive at the plateau configurations.

The trapping configurations are stable local minima: flat directions correspond to the moves that interchange the two elements of the turning pairs, and all other transpositions raise the energy. The number of states in the trap is therefore $`2^{numberofturningpairs}`$. Since the number of turning pairs grows linearly with the size of the system, the size of the trap grows exponentially with $`N`$. In appendix C, an estimate for the number of turning points (same as the number of turning pairs) for plateau configurations is derived, so the exponential growth is approximately $`2^{N/3}`$. Of course, many different traps exist, each with this typical size.

The picture of the trapping configurations in figure 10 makes the turning pairs appear like domain walls. This is a reasonable interpretation since the bulk part of the TSP energy has two different minima with ascending and descending order and it is these phases that are separated by the turning pairs. The boundary conditions raise the degeneracy of the sorted and antisorted phases, but this effect never comes into play here since the domain walls are frozen and do not move after the initial relaxation.

The local minima provide a basis for attempting to understand the value of the plateau energy observed under the dynamics above. A rather naive argument given in appendix C, based on characterizing the typical length of ascending or descending runs in the permutation defining the local minimum state gives the per site value of $`1/10`$. This should be compared with the numerical value of $`E_{TSP}/N=0.1092\pm 0.0001`$ found in section 4.2.

### 4.6 Interpolating Model

Given the rather different dynamic properties of the two models, it is natural to consider the behavior of models interpolating between them. The interpolating Hamiltonian is:

$$E(p)=(1p)E_{TSP}+pE_D$$
(14)
where $`p`$ is a parameter in (0,1).

The average energy (in the combinatoric sense at infinite temperature) is independent of $`p`$. All other quantities, such as the energy of individual states and the degeneracy pattern are different from either of the models previously discussed.

Depending on the value of $`p`$, the dynamics is expected to be more like one or other of the two lines shown in figure 3. One might hope for a transition to take place for some intermediate value of $`p`$, but in fact the dynamics reaches a plateau for any $`p<1`$. Figure 11 shows how the energy value of the plateau varies with $`p`$. For small $`N`$ this curve is composed of a series of steps, ending in a horizontal ledge near $`p=1`$. For larger $`N`$ the curve becomes smooth and the width of the ledge shrinks, eventually the curve becomes a sawtooth.

The question arises as to why stable local minima appear even for an infinitesimal perturbation away from the displacement energy. Consider the state corresponding to the completely reversed permutation (with the central two elements transposed in the case of $`N`$ even). The displacement energy of this state is flat with respect to any of the local moves, though a series of moves will eventually lead to a lower energy state. On the other hand the same state is part of a local minimum quartet with respect to the TSP energy. An infinitesimal addition of the TSP component is therefore sufficient to make the state a local minimum, albeit with infinitesimally small barriers. Although this argument correctly describes the reason that trapping can take place for $`p`$ so close to 1, the identification of completely reversed states as being responsible is incorrect since their energy is considerably higher than the plateau.

## 5 Conclusion

Considering statistical mechanical models based on states that are permutations of numbers, we have demonstrated a relationship between two models with distinct energy functions. Statically their partition functions are identical since the energy spectrum is the same. Dynamically there are substantial differences since one model has local energy minima and the other does not. This strange situation, that models with identical static properties have distinct local dynamics, is not a paradox. In this case it is clearly due to the shuffling of states that modifies the energy landscape by rearranging the energies of states that are close to each other.

The static analysis showed that the energy spectrum was the same by demonstrating a one-to-one map relating states with the same energy in each model. This map continues to hold for the more general model based on real numbers rather than integers.

Dynamically, we simulated a deep quench and showed that the energy decay proceeds in a completely different way in each of the two models. For the displacement model, no dynamical lengthscale appears as the system is reorganised. Characteristic timescales do appear and vary as $`N^2`$, though with strong finite size effects. For the TSP model, there was rapid (timescale varying as $`N`$) decay to a plateau as a result of reorganisation over small spatial scales of size less than 10 units, that did not completely destroy the corelation with the initial state. Evolution continued to occur on the plateau, but was limited in range. This was interpreted in terms of trapping configurations with turning pairs. Flat movement within the trap was still possible between about $`2^{N/3}`$ local minima states.

Optimization problems of most interest, such as K-Satisfiability, are much richer than the models presented here. Nevertheless, in the context of the one parameter family of interpolating models, we have shown that our simple local search procedure exhibits a transition (at the very edge of the interpolation region) from being able to sort the numbers to becoming trapped by local minima. We hope that these models provide a simple arena for studying the general question of the role of the energy landscape in the performance of local search algorithms.

## Appendix A Appendix: Model based on Real Numbers

Instead of the integer model discussed in the text, this more general model is based on the set of $`N`$ real numbers $`x_i`$, $`i=1,2\mathrm{}N`$, with each $`x_i`$ in the range $`[0,N+1]`$. Without any loss of generality, we take the $`x_i`$’s to be ordered ($`x_i<x_j`$ for all $`i<j`$) as this makes the identity permutation of the indices correspond to the lowest energy, or sorted state. The $`x_i`$’s should be regarded as quenched random variables, and for this reason the model might be regarded as a disordered model in contrast to the integer ordered model. The disorder however, is not of the independent variety that has been considered for both assignment and TSP problems using a replica approach in , but rather corresponds to Euclidean distances in one dimension. Here, the replica approach becomes intractable since the average over disorder couples sites.

For the disordered model, the definitions of the two energy functions remain exactly as given on the first lines of equations (2) and (3), however the distance function is replaced by:

$$d(P_i,P_j)=|x(P_i)x(P_j)|$$
(15)
Irrespective of the choice of $`x`$’s, there are degeneracies in the spectrum, but not to the extent found in the ordered model. For example, for $`N=4`$, the disordered model has energy levels with multiplicities (9,4,3,3,1,1,1,1,1). In the ordered model these combine to give (9,7,4,3,1).

The mapping between the energy levels according to each model continues to hold, however there is a small difference in the moments. Here the moments are defined after additional averaging over the $`x`$’s. The mean energy is the same as in the ordered model by design - indeed it was this criterion that selected the range of the $`x_i`$’s to be $`[0,N+1]`$. The second moment is however slightly larger than quoted in (11).

$$\overline{(E(N)E)^2}=\frac{(N+1)^2(3N^2+10N+2)}{(N+2)180N^2}$$
(16)
Even for individual instances of the disordered model, similar dynamical behaviour to that described in the text for the integer model is observed. After averaging over realizations of the $`x_i`$’s, the similarity becomes even closer.

## Appendix B Appendix: Multiplicities of Degeneracies

The multiplicities are most easily obtained using the displacement energy $`E_D`$ This may be evaluated for every permutation and $`D_N(E)`$ denotes the number of times that displacement $`E`$ appears. For small $`N`$ (up to $`N16`$), it is straightforward to numerically enumerate these coefficients, and the first few are given in the table below.

| $`N`$ | $`D(0)`$ | $`D(1)`$ | $`D(2)`$ | $`D(3)`$ | $`D(4)`$ | $`D(5)`$ | $`D(6)`$ | $`D(7)`$ | $`D(8)`$ | $`D(9)`$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 3 | 1 | 2 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 4 | 1 | 3 | 7 | 9 | 4 | 0 | 0 | 0 | 0 | 0 |
| 5 | 1 | 4 | 12 | 24 | 35 | 24 | 20 | 0 | 0 | 0 |
| 6 | 1 | 5 | 18 | 46 | 93 | 137 | 148 | 136 | 100 | 36 |

The usual arguments intended to lead to recurrence relations that would relate $`D_{N+1}(E)`$ to sums of $`D_N(E^{})`$ are not helpful in this case. Although there are some relations based on a sum over partitions of $`E_S`$, these are only valid in the region below the diagonal of the table. The problem appears to be non-trivial and indeed, according to Knuth, the generating function does not appear to have a simple form.

## Appendix C Appendix: TSP Plateau Energy

The average TSP energy is estimated after the effect of the Monte Carlo moves. It is convenient to evaluate the TSP energy by summing contributions from ascending and descending runs, rather than from adding each individual term which leads to many cancellations. The theory of runs of this kind is presented in , but very little of the general development is necessary here. The main property we use is that their endpoints are characterized by turning points in the permutation sequence. We assume that $`N`$ is large and only consider leading terms.

$$E_{TSP}\frac{1}{2N}\underset{i=0}{\overset{N}{}}|P_{i+1}P_i|\frac{1}{2N}s\mathrm{\Delta }$$
(17)
Where $`s`$ is the average number of runs, and $`\mathrm{\Delta }`$ is the average (absolute) change between the start and end of the run ($`\mathrm{\Delta }=d(P_{start},P_{end})`$). The first approximation in this approach is to ignore the correlations between the number of runs and the value of $`\mathrm{\Delta }`$ for that run in the formula above.

To illustrate this form, let us use it to reproduce the average energy of a random configuration. In that case, $`s=2N/3`$, since by considering sets of three adjacent points, the central one has a maximum or minimum value in 4 out of the 6 equally likely orderings. To deduce $`\mathrm{\Delta }`$ we consider the average value taken at upper turning points separating runs. Since $`max(P_1,P_2,P_3)=3N/4`$ and a symmetric result for the minimum, $`\mathrm{\Delta }=N/2`$. Combining these results in equation (17) reproduces $`E_{TSP}(N)=N/6`$ as was derived in the text by considering contributions from all neighboring pairs.

An estimate of the value of the average in the plateau is obtained by taking into account the effect of the Monte Carlo moves. We consider a subsequence of four points and look at the effect of a transposition on the central pair of the four, in all 24 possible orderings. The four points are labeled 1234, but this only signifies their relative magnitude. Table 1 indicates whether the transposition causes a positive, negative or zero change in energy and lists the change in number of (internal) maxima or minima.

Using this table, we find that moves that do not change the energy do not alter the number of runs. Moreover, there are six configuration that both reduce the energy and number of runs. These are the pair 1324, 4231, and the quartet 2314, 1423, 3241, 4132. If we take their naive weights from the initial random configuration, then such transpositions for every pair of adjacent points leads to an average reduction in $`s`$ of $`2\times N\times 2/24+1\times N\times 4/24`$. So the final value after one Monte Carlo sweep is $`s=N/3`$.

The argument for the average value at a maximum is extended from 3 (given in the example above) to 4 points, where it reproduces the 3 point case when all contributions are included. However, some of these contributions are removed by the Monte Carlo move (we only include those configurations with $`\mathrm{\Delta }E0`$ in table 1) resulting in $`\mathrm{\Delta }=3N/5`$.

Overall we obtain the estimate $`E_{TSP}(N)=1/2N\times N/3\times 3N/5=N/10`$, to be compared with the numerical value of $`E/N=0.1092\pm 0.0001`$. The approximations in this approach are quite drastic. First we ignored correlations between the number of runs and the value of $`\mathrm{\Delta }`$ for that run in formula (17). Then we just looked at the effect of one Monte Carlo sweep, used naive weights and ignored any influence of one Monte Carlo move on another. It is therefore surprising that we obtain such a reasonable estimate. Indeed, numerical studies show that the estimates of both $`s`$ and $`\mathrm{\Delta }`$ are incorrect by about 10%.