text:
# 1 Introduction

## 1 Introduction

In the seminal paper  F. Calogero described the cosmic origin of quantization. In paper  the tremor of the cosmic particles is the origin of the quantization and the characteristic acceleration of these particles $`a10^{10}`$ m/s2 was calculated. In our earlier paper  the same value of the acceleration was obtained and compared to the experimental value of the measured spacetime acceleration. In this paper we define the cosmic force –Planck force, $`F_{Planck}=M_Pa_{Planck}(a_{Planck}a)`$ and study the history of Planck force as the function of the age of the Universe.

Masses introduce a curvature in spacetime, light and matter are forced to move according to spacetime metric. Since all the matter is in motion, the geometry of space is constantly changing. A Einstein relates the curvature of space to the mass/energy density:

$$𝐆=k𝐓.$$
(1)
$`𝐆`$ is the Einstein curvature tensor and $`𝐓`$ the stress-energy tensor. The proportionality factor $`k`$ follows by comparison with Newton’s theory of gravity: $`k=G/c^4`$ where $`G`$ is the Newton’s gravity constant and $`c`$ is the vacuum velocity of light; it amounts to about $`2.10^{43}\mathrm{N}^1`$, expressing the rigidity of spacetime.

In paper  the model for the acceleration of spacetime was developed. Prescribing the $`G`$ for spacetime and $`+G`$ for matter the acceleration of spacetime was obtained:

$$a_{Planck}=\frac{1}{2}\left(\frac{\pi }{4}\right)^{1/2}\frac{(N+\frac{3}{4})^{1/2}}{M^{3/2}}A_P,$$
(2)
where $`A_P`$, Planck acceleration equal, viz.:

$$A_P=\left(\frac{c^7}{\mathrm{}G}\right)^{1/2}=\frac{c}{\tau _P}10^{51}ms^2.$$
(3)
As was shown in paper  the $`a_{Planck}`$ for $`N=M=10^{60}`$ is of the order of the acceleration detected by Pioneer spacecrafts .

Considering $`A_P`$ it is quite natural to define the Planck force $`F_{Planck}`$,

$$F_{Planck}=M_PA_P=\frac{c^4}{G}=k^1,$$
(4)
where

$$M_P=\left(\frac{\mathrm{}c}{G}\right)^{1/2}.$$

From formula (4) we conclude that $`F_{Planck}^1=`$ rigidity of the spacetime. The Planck force, $`F_{Planck}=c^4/G=\mathrm{1.2\hspace{0.33em}10}^{44}`$ N can be written in units which characterize the microspacetime, i.e. GeV and fm.

In that units

$$k^1=F_{Planck}=\mathrm{7.6\hspace{0.33em}10}^{38}\mathrm{GeV}/\mathrm{fm}.$$

## 2 The Planck, Yukawa and Bohr forces

As was shown in paper  the present value of Planck force equal

$$F_{Planck}^{Now}(N=M=10^{60})\frac{1}{2}\left(\frac{\pi }{4}\right)^{1/2}10^{60}\frac{c^4}{G}=10^{22}\frac{\mathrm{GeV}}{\mathrm{fm}}.$$
(5)
In paper  the Planck time $`\tau _P`$ was defined as the relaxation time for spacetime

$$\tau _P=\frac{\mathrm{}}{M_Pc^2}.$$
(6)
Considering formulae (4) and (6) $`F_{Planck}`$ can be written as

$$F_{Planck}=\frac{M_Pc}{\tau _P},$$
(7)
where $`c`$ is the velocity for gravitation propagation. In paper  the velocities and relaxation times for thermal energy propagation in atomic and nuclear matter were calculated:

$`v_{atomic}`$ $`=`$ $`\alpha _{em}c,`$ (8)
$`v_{nuclear}`$ $`=`$ $`\alpha _sc,`$
where $`\alpha _{em}=e^2/(\mathrm{}c)=1/137,\alpha _s=0.15`$. In the subsequent we define atomic and nuclear accelerations:

$`a_{atomic}`$ $`=`$ $`{\displaystyle \frac{\alpha _{em}c}{\tau _{atomic}}},`$ (9)
$`a_{nuclear}`$ $`=`$ $`{\displaystyle \frac{\alpha _sc}{\tau _{nuclear}}}.`$
Considering that $`\tau _{atomic}=\mathrm{}/(m_e\alpha _{em}^2c^2)`$, $`\tau _{nuclear}=\mathrm{}/(m_N\alpha _s^2c^2)`$ one obtains from formula (9)

$`a_{atomic}`$ $`=`$ $`{\displaystyle \frac{m_ec^3\alpha _{em}^3}{\mathrm{}}},`$ (10)
$`a_{nuclear}`$ $`=`$ $`{\displaystyle \frac{m_Nc^3\alpha _s^3}{\mathrm{}}}.`$
We define, analogously to Planck force the new forces: $`F_{Bohr}`$, $`F_{Yukawa}`$

$`F_{Bohr}`$ $`=`$ $`m_ea_{atomic}={\displaystyle \frac{(m_ec^2)^2}{\mathrm{}c}}\alpha _{em}^3=\mathrm{5\hspace{0.33em}10}^{13}{\displaystyle \frac{\mathrm{GeV}}{\mathrm{fm}}},`$ (11)
$`F_{Yukawa}`$ $`=`$ $`m_Na_{nuclear}={\displaystyle \frac{(m_Nc^2)^2}{\mathrm{}c}}\alpha _s^3=\mathrm{1.6\hspace{0.33em}10}^2{\displaystyle \frac{\mathrm{GeV}}{\mathrm{fm}}}.`$
Comparing formulae (5) and (11) we conclude that gradients of Bohr and Yukawa forces are much large than $`F_{Planck}^{Now}`$, i.e.:

$`{\displaystyle \frac{F_{Bohr}}{F_{Planck}^{Now}}}`$ $`=`$ $`{\displaystyle \frac{5.10^{13}}{10^{22}}}10^9,`$ (12)
$`{\displaystyle \frac{F_{Yukawa}}{F_{Planck}^{Now}}}`$ $`=`$ $`{\displaystyle \frac{10^2}{10^{22}}}10^{20}.`$
The formulae (12) guarantee present day stability of matter on the nuclear and atomic levels.

As the time dependence of $`F_{Bohr}`$ and $`F_{Yukawa}`$ are not well established, in the subsequent we will assumed that $`\alpha _s`$ and $`\alpha _{em}`$ do not dependent on time. Considering formulae (8) and (11) we obtain

$`{\displaystyle \frac{F_{Yukawa}}{F_{Planck}}}`$ $`=`$ $`{\displaystyle \frac{1}{(\frac{\pi }{4})^{1/2}}}{\displaystyle \frac{(m_Nc^2)^2}{M_Pc^2}}{\displaystyle \frac{\alpha _s^3}{\mathrm{}}}T,`$ (13)
$`{\displaystyle \frac{F_{Bohr}}{F_{Planck}}}`$ $`=`$ $`{\displaystyle \frac{1}{(\frac{\pi }{4})^{1/2}}}{\displaystyle \frac{(m_ec^2)^2}{M_Pc^2}}{\displaystyle \frac{\alpha _{em}}{\mathrm{}}}T.`$ (14)
As can be realized from formulae (13), (14) in the past $`F_{Planck}F_{Yukawa}`$ (for $`T=0.002`$ s) and $`F_{Planck}F_{Bohr}`$ (for $`T10^8`$ s), $`T=`$ age of universe. The calculated ages define the limits for instability of the nuclei and atoms.

## 3 The Planck, Yukawa and Bohr particles

In 1900 M. Planck  introduced the notion of the universal mass, later on called the Planck mass

$$M_P=\left(\frac{\mathrm{}c}{G}\right)^{1/2},F_{Planck}=\frac{M_Pc}{\tau _P}.$$
(15)
Considering the definition of the Yukawa force (11)

$$F_{Yukawa}=\frac{m_Nv_N}{\tau _N}=\frac{m_N\alpha _{strong}c}{\tau _N},$$
(16)
the formula (16) can be written as:

$$F_{Yukawa}=\frac{m_{Yukawa}c}{\tau _N},$$
(17)
where

$$m_{Yukawa}=m_N\alpha _{strong}147\frac{\mathrm{MeV}}{\mathrm{c}^2}m_\pi .$$
(18)
From the definition of the Yukawa force we deduced the mass of the particle which mediates the strong interaction – pion mass postulated by Yukawa in .

Accordingly for Bohr force:

$`F_{Bohr}`$ $`=`$ $`{\displaystyle \frac{m_ev}{\tau _{Bohr}}}={\displaystyle \frac{m_e\alpha _{em}c}{\tau _{Bohr}}}={\displaystyle \frac{m_{Bohr}c}{\tau _{Bohr}}},`$ (19)
$`m_{Bohr}`$ $`=`$ $`m_e\alpha _{em}=3.7{\displaystyle \frac{\mathrm{keV}}{\mathrm{c}^2}}.`$ (20)
For the Bohr particle the range of interaction is

$$r_{Bohr}=\frac{\mathrm{}}{m_{Bohr}c}0.1\mathrm{nm},$$
(21)
which is of the order of atomic radius.

Considering the electromagnetic origin of the mass of the Bohr particle, the planned sources of hard electromagnetic field, e.g. free electron laser (FEL)at TESLA accelerator (DESY)  are best suited to the investigation of the properties of the Bohr particles.

## 4 Possible interpretation of $`F_{Planck}`$, $`F_{Yukawa}`$ and $`F_{Bohr}`$.

In an important work, published already in 1951 J. Schwinger  demonstrated that in the background of a static uniform electric field, the QED vacuum is unstable and decayed with spontaneous emission of $`e^+e^{}`$ pairs. In the paper  Schwinger calculated the critical field strengths $`E_S`$:

$$E_S=\frac{m_e^2c^3}{e\mathrm{}}.$$
(22)
Considering formula (22) we define the Schwinger force:

$$F_{Schwinger}^e=eE_S=\frac{m_e^2c^3}{\mathrm{}}.$$
(23)
Formula (23) can be written as:

$$F_{Schwinger}^e=\frac{m_ec}{\tau _{Sch}},$$
(24)
where

$$\tau _{Sch}=\frac{\mathrm{}}{m_ec^2}$$
(25)
is Schwinger relaxation time for the creation of $`e^+e^{}`$ pair. Considering formulae (11) the relation of $`F_{Yukawa}`$ and $`F_{Bohr}`$ to the Schwinger force can be established

$`F_{Yukawa}`$ $`=`$ $`\alpha _s^3\left({\displaystyle \frac{m_N}{m_e}}\right)^2F_{Schwinger}^e,\alpha _s=0.15,`$ (26)
$`F_{Bohr}`$ $`=`$ $`\alpha _{em}^3F_{Schwinger}^e,\alpha _{em}={\displaystyle \frac{1}{137}},`$
and for Planck force

$$F_{Planck}=\left(\frac{M_P}{m_e}\right)^2F_{Schwinger}^e.$$
(27)
In Table 1 the values of the $`F_{Schwinger}^e`$, $`F_{Planck}`$, $`F_{Yukawa}`$ and $`F_{Bohr}`$ are presented, all in the same units GeV/fm. As in those units the forces span the range $`10^{13}`$ to $`10^{38}`$ it is valuable to recalculate the Yukawa and Bohr forces in the units natural to nuclear and atomic level. In that case one obtains:

$$F_{Yukawa}\frac{16\mathrm{MeV}}{\mathrm{fm}}.$$
(28)
It is quite interesting that $`a_v16`$ MeV is the volume part of the binding energy of the nuclei (droplet model).

For the Bohr force considering formula (11) one obtains:

$$F_{Bohr}\frac{50\mathrm{eV}}{0.1\mathrm{nm}}.$$
(29)
Considering that the Rydberg energy $`27`$ eV and Bohr radius $`0.1`$ nm formula (29) can be written as

$$F_{Bohr}\frac{\mathrm{Rydberg}\mathrm{energy}}{\mathrm{Bohr}\mathrm{radius}}$$
(30)
## 5 Concluding remarks

In this paper the forces: Planck, Yukawa and Bohr were defined. It was shown that the present value of the Planck force (which is the source of the universe acceleration) $`10^{22}`$ GeV/fm is much smaller than the Yukawa ($`10^2`$ GeV/fm) and Bohr ($`10^{13}`$ GeV/fm) forces respectively. This fact guarantees the stability of the matter in the present. However in the past for $`T`$ (age of the universe), $`T<0.002`$ s, $`F_{Yukawa}<F_{Planck}`$ (0.002 s) and $`F_{Bohr}<F_{Planck}`$ ($`10^8`$ s). In this paper the relation of the Schwinger force (for the vacuum creation of the $`e^+e^{}`$ pairs) to the Planck, Yukawa and Bohr force was obtained.


format: text

text:
# I Structural parameters calculated for the orthorhombic phases of ZrO2 and HfO2. Lattice constants 𝑎, 𝑏 and 𝑐, in Å are given in the first row. Triplets of coordinates indicate the Wyckoff positions of the various atoms.

First-principles Study of Electronic and Dielectric Properties of ZrO2 and HfO2

Xinyuan Zhao and David Vanderbilt
Department of Physics and Astronomy,
Rutgers University,
Piscataway, NJ 08854-8019

ABSTRACT

Using density-functional theory with ultrasoft pseudopotentials, we previously investigated the structural and electronic properties of the low-pressure (cubic, tetragonal, and monoclinic) phases of ZrO2 and HfO2, in order to elucidate phonon modes, Born effective charge tensors, and especially the lattice dielectric response in these phases. We now extend this previous work by carrying out similar calculations on the two high-pressure orthorhombic phases, and by providing density-of-states and band-gap information on all polymorphs. Our results show that the electronic structures and dielectric responses are strongly phase-dependent. In particular, the monoclinic phases of ZrO2 and HfO2 are found to have a strongly anisotropic dielectric tensor and a rather small orientational average ($`\overline{ϵ_0}`$) compared to the two other low-pressure phases. Our calculations show that $`\overline{ϵ_0}`$ is even smaller in the orthorhombic phases.

INTRODUCTION

Dielectrics comprise a class of materials of great technical importance. The most widely-used dielectric in electronics is probably SiO2, which has been persistently exploited in the computer industry for decades due to the excellent quality of the Si/SiO2 interface. The dimensions of the gate-dielectric layer (currently SiO2) have been rapidly reduced to the nanometer scale in order to enhance the speed of CMOS-related devices. Such reductions, however, have uncovered severe fundamental problems that imposed a physical limit on the use of SiO2 films as gate dielectrics in the current Si/SiO2 semiconductor technology. A great deal of effort has been directed to finding a way to overcome these limits. One particular direction is to use high-$`K`$ metal oxides with certain desired properties  to replace SiO2 as the gate dielectric, since the higher dielectric constants of the oxides will allow one to make a physically thicker film while at the same time maintaining (or even increasing) the gate capacitance.

ZrO2, HfO2 and their structure-modified derivatives (e.g., Zr and Hf silicates) have emerged as important candidates for this purpose because they have much higher dielectric constants than SiO2 ($``$20 compared with 3.5 for SiO2) and because of their thermodynamic stability in contact with Si. Experiments have demonstrated their promise for this purpose (see Ref.  for a general review). It would be very desirable, therefore, to conduct a systematic theoretical study of the structural, electronic, and dielectric properties of ZrO2 and HfO2, and of the dependence of these properties on the choice of crystalline polymorphs or possible modified structures. In two previous works, we studied the structural and vibrational properties of the three low-pressure phases (cubic, tetragonal, and monoclinic) of ZrO2  and HfO2 , especially the zone-centered phonon modes and the related dielectric properties. These previous works are extended here in two ways. First, we provide density-of-states (DOS) information on all polymorphs of ZrO2 and HfO2, which may be helpful in the interpretations of experiments on these materials. Second, we extend our previous calculations of dielectric properties to include the high-pressure orthorhombic phases using the ABINIT package . We have also been engaged in simulating amorphous ZrO2 ($`a`$-ZrO2) using ab-initio molecular dynamics as implemented in the VASP  package, in a “melt-and-quench” fashion. The resulting amorphous ZrO2 structures were analyzed in order to understand their local bonding geometry and other structural properties, and then used as inputs for a first-principles investigation of the vibrational and dielectric properties of the amorphous ZrO2 systems using ABINIT. The results will be reported elsewhere.

CRYSTAL STRUCTURES

ZrO2 (zirconia) is surprising similar to HfO2 (hafnia) in many physical and chemical properties. It is well established that both oxides exhibit three crystalline phases at ambient pressure, i.e., monoclinic, tetragonal and cubic, as temperature is increased. With increasing pressure, phase transitions involving two orthorhombic phases have been reported, although the exact structures of these high-pressure phases has not been established unambiguously (, and references therein). The first orthorhombic structure used in our simulation (conventionally labeled as Ortho-I) has space group $`Pbca`$; the other orthorhombic phase (Ortho-II) has space group $`Pnma`$. We study these phases at zero pressure, where they are apparently metastable.

The cubic phase has the fluorite (CaF2) structure, with 3 atoms (one chemical formula) per primitive unit cell. The tetragonal structure is a slight distortion of the cubic phase obtained by shifting alternating pairs of oxygen atoms up and down along $`\widehat{z}`$, and the number of atoms per primitive cell is doubled. The monoclinic phase, which has 12 atoms (four formulas units) per primitive cell, has a more complex structure and lower symmetry. In our calculations of the DOS functions, a 12-atom cell has been used for all of these phases, as in Refs.  and , in order to facilitate comparison. The monoclinic phase has two non-equivalent oxygen atoms, which evidently will have different site-projected DOS functions. Among the eight oxygens in the unit cell, four are threefold coordinated (denoted as O1), while the other four (labeled as O2) are fourfold coordinated. Each threefold oxygen is bonded to three nearest-neighbor metal atoms in an almost planar configuration, while each four-fold oxygen forms a distorted tetrahedron with its four nearest metal neighbors. For the cubic and tetragonal phases, all oxygen atoms are equivalent, although we still use the notation O1 and O2 for convenience. Detailed discussion of the structures of these phases can be found in Refs.  and .

A primitive unit cell of the Ortho-I phase (for both ZrO2 and HfO2) consists of 24 atoms, while the Ortho-II phase has 12 atoms per unit cell. We used the lattice parameters given in  (Table III there) and  (Tables II and III there) as the starting point for our ground-state calculations for Ortho-I and Ortho-II, respectively. The Wyckoff coordinates define two non-equivalent sets of oxygen atoms, labeled by O1 and O2.

THEORETICAL METHODOLOGY

The DOS functions are computed using a plane-wave implementation of density functional theory (DFT) in the local-density approximation (LDA) as parametrized by Ceperley and Alder , using ultra-soft pseudopotentials . Although a 4$`\times `$4$`\times `$4 k-point sampling was found to provide sufficient precision in calculating total energies and forces, we used a refined mesh of k-points and higher 36-Ry cut-off to obtain smoother and more accurate DOS functions. The k-point meshes used in the DOS calculations are 20$`\times `$20$`\times `$20 for cubic and tetragonal, 10 $`\times `$10$`\times `$10 for monoclinic, 12$`\times `$12$`\times `$12 for Ortho-I, and 16$`\times `$16$`\times `$16 for Ortho-II phases. In calculating the site-projected partial DOS, spherical radii of 0.9 Å and 1.25 Å are used for the metal ions (Zr or Hf) and oxygen ions respectively. For cubic and tetragonal phases, the larger 12-atom unit cell was used in the calculations, so the DOS values should be reduced accordingly if the DOS per primitive cell is desired. Further details can be found in  and .

The dielectric response of orthorhombic phases of ZrO2 was calculated by specialized linear-response techniques as implemented in ABINIT , taking into account the coupling between phonons and the homogeneous electric field. A 4$`\times `$4$`\times `$4 k-point mesh, Troullier-Matins type pseudopotentials  (4 and 6 valence electrons for Zr and O respectively), and an energy cutoff of 35.0 Hartree provide satisfactory convergence.

RESULTS

In Table I we present our calculated lattice parameters for the two orthorhombic phases, which agree well with the previous theoretical calculations . Because the computed DOS functions for HfO2 and ZrO2 are so similar, we only present those for HfO2 here, shown in Figs. (1)-(3). Some comparisons are given later in Table II for both oxides. Figure 1 shows the total DOS (lowest panel) and the site-projected partial DOS (upper three panels) for monoclinic HfO2. We find this material to be an insulator with a band gap of 3.5 eV, in good agreement with previous theoretical results (3.48 eV in  and 3.6 eV in ). The bands in the energy range from $``$20 eV to $``$15 eV are mostly due to O 2$`s`$ orbitals. The next bands, spanning between $``$6 eV and 0, are mainly O 2$`p`$ in character, and Hf 5$`d`$ orbitals dominate the unoccupied conduction bands ranging from 3 eV up to 11 eV. The two non-equivalent oxygen sites have noticeably different DOS features, as expected from their different surroundings.

In Fig. 2, we provide the DOS functions calculated for cubic and tetragonal HfO2. It can be readily seen that these two phases have rather similar DOS features. Cubic HfO2 has a band gap of 3.2 eV, while for tetragonal HfO2 it is 3.8 eV. Ref.  gives 3.4 eV and 3.77 eV for the cubic and tetragonal phases, respectively. Again, both theoretical works agree quite well. The most significant difference between these two phases and the monoclinic phase is that the Hf $`d`$-like states separate into two bands in the cubic and tetragonal phases, but form a single band in the monoclinic phase. The DOS functions for the two high-pressure orthorhombic phases are presented in Fig. 3, from which one readily sees that the oxygen atoms exhibit two quite different types of DOS, as in the monoclinic phase. In fact, there is a surprising degree of resemblance of the monoclinic DOS with that of the Ortho-I and Ortho-II phases, especially the former.

The DOS functions of the ZrO2 polymorphs (not shown) are qualitatively very similar to the corresponding ones for HfO2. The most important differences are in the band gaps. Because of LDA error, the absolute values of the band gaps are not reliable, but trends should be meaningful. As shown in Table II, we find that the band gaps are systematically larger (by $``$0.5 eV) for HfO2 than for ZrO2, and that variations of crystal structure can lead to band gap changes of order 1 eV. Also tabulated in Table II are our calculated dielectric constants (orientational average $`\overline{ϵ_0}`$) for ZrO2 and HfO2 polymorphs. For the low-pressure cubic, tetragonal, and monoclinic phases, the dielectric constants are quoted from  and , where the full dielectric tensors are given. For the Ortho-I and II phases of ZrO2, the linear-response scheme  is utilized in the calculation of the static dielectric tensors, which includes the lattice contributions as well as the electronic screening. In orthorhombic ZrO2, $`\overline{ϵ_0}`$ becomes even smaller, approximately 20 and 18 for the Ortho-I and II phases, respectively. The dielectric tensors for the orthogrombic phases are diagonal, with elements (22.6, 18.1, 19.6) and (18.8, 18.9, and 17.8) for Ortho-I and Ortho-II respectively. The average principal values of the Born effective charge tensors are about 5.0 and 2.5 for Zr and O atoms respectively for both orthorhombic phases.

CONCLUSIONS

We studied the structural, electronic, vibrational and dielectric properties for nearly all the currently well-recognized phases of ZrO2 and HfO2, using advanced first-principles techniques, with a special focus on the dielectric properties. It is found that the crystalline structure can have a substantial effect on the band gaps and dielectric constants. Our results show that the dielectric response depends strongly on the crystal phase, and can span a wide range of values ($`\overline{ϵ_0}`$ = 18 – 40 for ZrO2). Monoclinic ZrO2 (or HfO2) has a rather anisotropic dielectric tensor and a smaller averaged dielectric constant than the tetragonal and cubic phase. The $`\overline{ϵ_0}`$ of the Ortho-I and II phases appear even smaller. We hope that these theoretical results will shed light on the potential for ZrO2, HfO2 or their derivatives to replace SiO2 as the gate dielectric in CMOS technology, and will provide hints for manipulating these oxides to serve this purpose.

ACKNOWLEDGMENTS

This work has been supported by NSF Grant DMR-9981193. We wish to thank E. Garfunkel and S. Sayan for useful discussions. One of authors (X.Z.) thanks G.-M. Rignanese for assistance with ABINIT.

REFERENCES

1. G. D. Wilk, R. M. Wallace, and J. M. Anthony, J. Appl. Phys. 89, 5243 (2001).
2. X. Zhao and D. Vanderbilt, Phys. Rev. B 65, 75105 (2002).
3. X. Zhao and D. Vanderbilt, Phys. Rev. B 65 233106 (2002).
4. X. Gonze et al., Mater. Sci. 25, 478 (2002).
5. G. Kresse and J. Hafner, Phys. Rev. B 47, R558 (1993); 54, 11169 (1996).
6. J. K. Dewhurst and J. E. Lowther, Phys. Rev. B 57, 741 (1998).
7. G. Jomard, T. Petit, A. Pasturel, L. Magaud, G. Kresse, and J. Hafner, Phys. Rev. B 59, 4044 (1999).
8. D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
9. D. Vanderbilt, Phys. Rev B 41, 7892 (1990).
10. N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
11. A. A. Demkov, Phys. Stat. Sol. B 226, 57 (2001).
12. P. K. Boer and R. A. de Groot, J. Phys.: Condens. Matter 10, 10241 (1998).


text:
# 1 Introduction

## 1  Introduction

The time series of heartbeat intervals (RR intervals), used in various analyses, are usually Holter type data or data from steady state ambulatory measurements. The great wealth of data about the dynamics of the heart that is contained in such ECG records is usually reduced to characterize only the mean heart rate and the presence and frequency of some abnormal electrocardiographic complexes. The analysis of short-range and long-range fluctuations and its universality are largely ignored.1) The normal human heartbeats display considerable variability and nonstationarity over time. The result is that fluctuations of heartbeat intervals around the mean value are irregular and complex, even under resting conditions, Fig. 1. A number of new methods have been proposed to quantify these complex fluctuations, mainly based on the nonlinear and fractal analysis. 2),3),4),5)

In this paper we examine the possibility that the distribution of heartbeat intervals exhibits universality of the type predicted by the RMT. 6)

## 2 Statistical analysis of heartbeat intervals

The heartbeat time series of an human represents an ordered sequence of R-waves (QRS-complexes) occuring at times $`\{t_1,t_2,\mathrm{},t_i,\mathrm{}\}`$. In this form the heartbeat dynamics can be regarded as a fractal stochastic point process on a time-line. The point process is completely described by the set $`\{t_i\}`$, or equivalently by the set of RR-intervals defined as $`\{RR(i)=t_{i+1}t_i,i=1,2,\mathrm{}\}`$.

To see the universality in beat-to-beat RR interval fluctuations we need to magnify the heartbeat spectrum $`\{t_i\}`$ so that mean spacing is unity. Let us first determine the number of heartbeats up to time $`t`$. It is given by the staircase function $`N(t)`$:

$$N(t)=\underset{i}{}\theta (tt_i)$$
(1)
where $`\theta (t)`$ is the step function.

In order to be able to compare statistical properties of heartbeat spectra, in particular the fluctuations around the mean value, with the predictions of the RMT it is necesary to normalize (unfolde) the time spectrum $`\{t_i\}`$ to $`\{\tau _i\}`$ in such a way that the mean spacing between neighbouring $`\tau _i`$ is unity. This is achieved by the mapping

$$\tau _i=\overline{N}(t_i),S_i=\tau _{i+1}\tau _i$$
(2)
where $`\overline{N}(t)`$ denotes the smoothed staircase which account for the mean trend in $`N(t)`$ and $`S_i`$ is the normalized heartbeat spacing. The ”local” mean spacing is defined as $`\mathrm{\Delta }(t)=[d\overline{N}(t)/dt]^1`$ so that $`S_i`$ becomes

$$S_i=\frac{t_{i+1}t_i}{\mathrm{\Delta }[(t_{i+1}+t_i)/2]}.$$
(3)
During the real time measurements, we always deal with the finite number of heartbeats, say $`n`$. In that case, we define the ”global” mean spacing as

$$\mathrm{\Delta }=\frac{1}{n1}\underset{i=1}{\overset{n1}{}}RR(i)=\frac{t_nt_1}{n1}$$
(4)
so that $`S_i=RR(i)/\mathrm{\Delta }`$.

There are two important statistical measures that we shall use to characterize the heartbeat spectra $`\{\tau _i\}`$. One is the nearest-neighbor heartbeat spacing distribution $`P(S)dS`$ which gives the probability to find $`\tau _i`$ in the interval ($`S,S+dS`$). $`P(S)`$ is a short range statistics, it measures the short range correlations between heartbeats. Its behaviour for small spacing, $`S0`$,

$$P(S)const\times S^\beta $$
(5)
determines the type of short range correlations between heartbeat intervals. Two cases can be distinguished: a) when $`\beta =0`$ meaning that the short range fluctuations of $`\{S_i\}`$ are just those of uncorrelated random numbers, i.e. Poissonian with a well-known phenomenom of attraction (clustering), and b) when $`\beta 0`$ meaning that the long range fractal fluctuations are present in the sequence of $`S_i`$ that is characterized by the phenomenom of repulsion. The universal statistics (no free parameters) predicted by the RMT 6) are given by

$`P_{Poisson}(S)`$ $`=`$ $`exp(S),\beta =0,`$
$`P_{GOE}(S)`$ $`=`$ $`{\displaystyle \frac{\pi }{2}}Sexp({\displaystyle \frac{\pi }{4}}S^2),\beta =1,`$ (6)
$`P_{GUE}(S)`$ $`=`$ $`{\displaystyle \frac{32}{\pi ^2}}S^2exp({\displaystyle \frac{4}{\pi }}S^2),\beta =2,`$
The GOE (Gaussian orthogonal ensemble) distribution, known as a Wigner distribution, should be valid for irregular (chaotic) systems with an antiunitary symmetry, whereas GUE (Gaussian Unitarian Ensemble) should apply if the system has no such symmetry.

To avoid the smoothing procedure on a measured data, it is easier and less ambiguous to study the inegrated or the cumulative RR intervals

$$I(S)=_0^SP(S)𝑑S$$
(7)
which yields a useful statistics even for small sample of RR intervals. This gives

$$I_{Poisson}(S)=1exp(S),I_{Wigner}(S)=1exp(\frac{\pi }{4}S^2),$$
(8)
and a nonanalytical expression for $`I_{GUE}(S)`$.

Another important statistics, not considered here, is the spectral rigidity $`\mathrm{\Delta }_3`$ of Dyson and Metha 7) which characterizes the long-range correlations of heartbeats. It can be related to the Hurst exponent 8),9) associated with the heartbeat time series.

## 3 Results and discussion

Time series of RR intervals were measured at the Institute for Cardiovascular Disease and Rehabilitation in Zagreb. Each masurement of ECG was over a time duration of about 15 min. ($``$ 2000 heartbeats), in a controlled ergometric regime. These type of measurements are used as routine in everyday clinical diagnostic practice, because some heart deseases, such as stable angine pectoris (SAP), usually become transparent under physical activities.

The ECG data were digitized by the WaveBook 512 (Iotech. Cal. USA), and transferred to a computer. The RR interval series was passed through a filter that eliminates noise and artefacts. All R-wave peaks were first edited automatically, after which a careful manual editing was performed by visual inspection of the each RR interval. After this, all questionable portions were excluded manually, and only segments with $`>90\%`$ sinus beats were included in the final analysis.

Each measurement consisted of stationary state part (pretrigger Pt), a few stages of running (P1-P4) and a period of relaxation (Re), Fig. 1.

Here we report on the universality analysis of RR interval fluctuations in the pretrigger and P1 period only. 10) Our patients were divided in two groups: one with the evidence of ishemic ST-segment depression (SAP subjects), and the control group of healthy subjects. Selection of subjects was performed by a cardiologist according to the generally accepted medical knowledge.

The integrated RR interval distributions $`I(S)`$ are shown in Fig.2 together with $`I_W(S)`$.

Clearly the universal, parameter free, distributions $`P_{Poisson},P_W`$ and $`P_{GUE}`$ of RMT are excluded. A possible reason for this is that heart is a very complex system of nonlinearly coupled biological oscillators whose motion is partly regular and partly irregular (chaotic) or mixed . Another possibility is that universality is only local.10)

The first theoretical step towards the understanding of the regular-irregular coexistence in a dynamical system was the work by Berry and Robnik 11) in which they offered a semiclassical model for $`P(S)`$ known as the Berry-Robnik formulae. It is to be applied at large intervals $`S>1`$. At small and intermediate intervals $`S<1`$ the Brody distribution 12)

$$P_B(S)=aS^\beta exp(bS^{\beta +1}),a=(\beta +1)b,b=\{\mathrm{\Gamma }(\frac{\beta +2}{\beta +1})\}^{\beta +1}$$
(9)
gives better fitting results. The Brody distribution becomes Poisson at $`\beta =0`$ and Wigner for $`\beta =1`$.

Figure 3 shows that the Brody distribution fits very well the ECG measured distribution of RR intervals. Each subject has its own value of $`\beta `$ which might be an indicator of a dynamical state of the heart. Indeed, Fig. 4 clearly shows that $`\beta `$ could be a useful parameter for separating SAP subjects from healthy ones.

## 4 Conclusion

In this paper we have investigated the possibility that fluctuations of RR intervals exhibit universal statistical laws. It is found, Fig. 2, that this is not so. Heart is a complex nonlinear system with many coexisting regular and irregular motions which manifest themselves through the nonuniversal fluctuations of RR intervals. It is also found, Fig. 3 and 4, that the one-parameter Brody distribution could be used successfuly to describe the fluctuation pattern of RR intervals. The parameter $`\beta `$ of the Brody distribution is seen, Fig. 4,to be connected with the dynamical state of the heart.

Further studies in larger populations are needed to confirm these interesting results.

## Acknowledgements

This work was supported by the Ministry of Science and Technology of the Republic of Croatia.


text:
# Cherenkov radiation emitted by ultrafast laser pulses and the generation of coherent polaritons

## I Introduction

Polaritons are hybrid modes of solids and, more generally, dense systems associated with elementary excitations that carry electromagnetic polarization.Born and Huang (1954); Mills and Burstein (1974) Such excitations couple to light, resulting in a new mixed particle called the polariton. The polariton description is, in some sense, an alternative to that involving the index of refraction; polaritons are dressed photons. Phonon polaritons, lattice vibrations coupled to infrared light, were first studied experimentally by spontaneous Raman scattering.Henry and Hopfield (1965) Soon after, mixing of coherent polaritons with visible laser light was reported for GaP.Faust and Henry (1966) The generation, dispersion and decay of phonon polaritons have been extensively investigated using various coherent optical methods, particularly, coherent anti-Stokes Raman scattering (CARS)Shen and Bloembergen (1965) followed by time-domain CARSGale et al. (1986) and, more recently, the closely related technique of impulsive stimulated Raman scattering (ISRS).

This work centers on ISRS, i.e., coherent scattering using pulses whose bandwidth exceeds the polariton frequency.Yan and Nelson (1987) In ISRS experiments, a pump pulse imparts an impulsive force on the phonon component of the polariton starting a coherent oscillation which in turn, through the inverse process, perturbs the index of refraction.Dhar et al. (1994); Merlin (1997) This perturbation is measured by a probe pulse that follows behind the pump pulse at a controllable time delay. The most common ISRS approach to generate coherent polariton fields is to use an intensity grating obtained from crossing two pump beams. This geometry leads to polaritons of nearly well-defined wavevector even in the pulsed case where the (transient) grating travels with the pulses.Wiederrecht et al. (1995); Bakker et al. (1998a) An alternative method, associated with the nonlinear susceptibility $`\chi ^{(2)}`$, is to excite the electromagnetic component of the polariton through difference frequency generation (DFG), or optical rectification.Bass et al. (1962) Here, the pump pulse induces a nonlinear polarization proportional to its intensity envelope. If the pulse is tightly focused and its group velocity is greater than the phase velocity in the infrared, its polarization leads to emission of infrared light in much the same way that a relativistic dipole emits (coherent) Cherenkov radiation (CR). This interpretation, originally proposed for second-harmonic generation,Tien et al. (1970) was put forth for DFG by Auston and coworkersAuston et al. (1984); Kleinman and Auston (1984); Auston and Nuss (1988) whose technique became later the standard for generating short THz pulses.Hu et al. (1990). Recently, the CR interpretation was revisited in the analysis of pump-probe experiments on ZnSeStevens et al. (2001) and ZnTe.Wahlstrand et al. (2002) In these and other polar materials, strong dispersion near the frequency of the transverse-optical (TO) phonon results in two qualitatively distinct regimes for CR depending on whether the velocity of the source is larger or smaller than the phase velocity of light at zero frequency. In the former, superluminal regime, the CR field is qualitatively the same as in the experiments of Auston and co-workers Auston and Nuss (1988) whereas the latter subluminal regime displays new features. Afanasiev et al. (1999)

In the ultrafast literature, the relationship between the two, transient grating and tightly-focused single pump methods is rarely discussed. While the generation mechanisms are closely related, ISRS is often described as four-wave mixing and, as such, as a technique that distinguishes itself from DFG. Here, we present a unified approach encompassing the two methods. We show that the generation of THz radiation by single or multiple pulses of arbitrary shape can always be described as CR. We also report experimental results on three materials, ZnTe, GaP, and LiTaO3, in three geometries (point, planar, and transient grating) using an electro-optic field sampling technique based on that of Auston and co-workers.Auston and Nuss (1988) This method allows us to map out the polariton field as a function of time delay and lateral position for polarization sources of arbitrary shape and, simultaneously, to image the shape of the source pulse itself.

## II Phonon polaritons and optical susceptibilities

As mentioned earlier, polaritons can be described either as excitations that carry polarization through the coupling to the electromagnetic field, or as light perturbed by the refractive index of the medium, which exhibits strong dispersion near the resonant TO frequency. In the Cherenkov picture of phonon polariton generation we take the latter view, incorporating the effect of the infrared-active phonon into the dielectric function and the second-order nonlinear susceptibility. For a diatomic cubic lattice, and neglecting damping, the dielectric function is of the Lorentz form

$$ϵ(\mathrm{\Omega })=ϵ_{\mathrm{}}+\frac{ϵ_0ϵ_{\mathrm{}}}{1(\mathrm{\Omega }/\mathrm{\Omega }_{\text{TO}})^2},$$
(1)
where $`\mathrm{\Omega }_{\text{TO}}`$ is the TO frequency of the three-fold degenerate infrared-active mode, $`ϵ_{\mathrm{}}`$ is the dielectric constant including the effects of higher-lying resonances, and $`ϵ_0`$ is the static dielectric constant. The Lyddane-Sachs-Teller (LST) formula $`\mathrm{\Omega }_{\text{LO}}^2=\frac{ϵ_0}{ϵ_{\mathrm{}}}\mathrm{\Omega }_{\text{TO}}^2`$ gives the relationship between $`\mathrm{\Omega }_{\text{LO}}`$, the longitudinal optical (LO) phonon frequency, and $`\mathrm{\Omega }_{\text{TO}}`$. The refractive index $`n(\mathrm{\Omega })=\sqrt{ϵ(\mathrm{\Omega })}`$, phase velocity $`c/n(\mathrm{\Omega })`$, and polariton dispersion relation $`\mathrm{\Omega }=cq/n(\mathrm{\Omega })`$ are plotted in Fig. 1 for a material with $`ϵ(\mathrm{\Omega })`$ given by Eq. (1). In this work we deal mainly with frequencies below $`\mathrm{\Omega }_{\text{TO}}`$. Near zero frequency, where the dispersion curve is linear, the polariton is mostly light-like whereas, near $`\mathrm{\Omega }_{\text{TO}}`$, it is phonon-like. When there is more than one infrared-active phonon, as in LiTaO3, each mode contributes to the polarization and the dielectric function reflects a sum of these contributions (poles).

Nonlinear interaction of visible light with phonon polaritons involves both the lattice and the THz electromagnetic field. The nonlinear tensor that couples polaritons and visible light is thus composed of two parts: the Raman tensor $`\frac{d\chi }{dQ}=\chi ^\mathrm{R}`$, which couples visible light to the phonon mode with amplitude $`Q`$, and the second-order nonlinear susceptibility $`\frac{d\chi }{dE}=\chi ^{(2)}`$, which couples visible to infrared electric fields. Both of these susceptibilities are almost entirely electronic in nature. Counting both contributions, the total nonlinear polarization isFaust and Henry (1966)

$$P_i^{\text{NL}}=\chi _{ijk}^{(2)}E_j(\mathrm{\Omega })E_k(\omega _L)+\chi _{ik}^\text{R}Q_jE_k(\omega _\text{L})$$
(2)
where $`Q_j`$ is the amplitude of one of the three-fold degenerate phonon modes, $`\omega _L`$ is the frequency of the visible laser light and $`\mathrm{\Omega }`$ is the frequency of the phonon polariton. As is well-known, the Raman contribution can be included in an effective nonlinear susceptibility $`d_{ijk}`$, each tensor element of which is of the form

$$d(\mathrm{\Omega })=\chi ^{(2)}[1+C(1\frac{\mathrm{\Omega }^2}{\mathrm{\Omega }_{\mathrm{TO}}^2})^1],$$
(3)
where the Faust-Henry coefficient $`C=\chi ^\mathrm{R}(\chi ^{(2)})^1e^{}\mu ^1(\mathrm{\Omega }_{\mathrm{TO}})^2`$ gives the ratio of the Raman (often termed ionic) to the electronic contribution.Faust and Henry (1966) Here, $`e^{}`$ is the transverse effective charge and $`\mu `$ is the reduced mass of the lattice mode. Because $`\chi ^\mathrm{R}`$ and $`\chi ^{(2)}`$ do not change appreciably with frequency in the THz range, $`C`$ can be approximated as a constant for a given laser wavelength and material. If there is more than one phonon mode, there is a Raman contribution to $`d`$ and to the Faust-Henry coefficient from each one. In our experiments we do not need to consider the dispersion of $`d`$, for we mainly deal with low frequency polaritons. Thus, in the following we assume that the nonlinear susceptibility is constant, and approximately given by $`\chi ^{(2)}(1+C)`$.

## III Theory of Cherenkov emission in an infinite crystal

In non-centrosymmetric media, an optical pulse induces through $`\chi ^{(2)}`$ a low-frequency polarization proportional to its intensity. This is DFG, also known as optical rectification.Bass et al. (1962) The dipolar charge distribution moves at the optical pulse’s group velocity $`v_\text{g}`$ and, like a charged particle, it can emit CR if $`v_g`$ is greater than the phase velocity of the THz radiation in the medium. The bandwidth of the THz radiation emitted by a pulse is limited to the pulse bandwidth. The spatial and temporal shape of the pulse envelope, together with the optical properties of the medium, determine the THz field.

Kleinman and Auston adapted the theory of CR to an optical pulse traveling through a nonlinear isotropic medium.Kleinman and Auston (1984) Their approach, which we build on here, is essentially to separate the nonlinear mechanism by which the polarization in the medium is generated from the Cherenkov mechanism by which infrared light is emitted. This approach works well in the absence of higher-order effects such as self-focusing and self phase modulation, which affect the spatiotemporal shape of the pulse as it propagates, and depletion, which diminishes the induced polarization. Neglecting these effects, the emission process is entirely analogous to the emission of CR by a continuum of dipoles traveling at the group velocity of the pulse. In the following, we discuss electro-optic CR by pulses of a variety of spatiotemporal shapes.

### III.1 Point source

For a beam focused to a small waist compared to the polariton wavelength, the polarization induced by the optical pulse can be approximated by a point dipole. This geometry, the 3D case discussed by Kleinman and Auston,Kleinman and Auston (1984) is the closest electro-optic analogue to conventional CR (emitted by a relativistic monopole), which is coherent. The theory of CR was largely worked out more than 60 years ago by Tamm and Frank.Frank and Tamm (1937); Tamm (1939) For an infinite medium and pulses propagating along the $`z`$ axis, the fields are a function of $`zv_gt`$ and the cylindrical coordinates $`\rho `$ and $`\varphi `$. This approximation, valid for a thick sample, greatly simplifies the theory.

The radiation field due to a propagating dipole can be derived from the radiation field due to a monopole by taking a derivative in the direction of the dipole orientation.Zrelov (1970) In the experiments discussed here, the induced dipole is oriented perpendicular to the $`z`$ axis. For this orientation, the fields are proportional to $`\mathrm{cos}\varphi `$, where $`\varphi `$ is measured from the direction of the dipole, dropping for simplicity the $`\varphi `$ dependence, we haveFrank and Tamm (1937); Tamm (1939); Zrelov (1970)

$`E_z(\rho ,zv_gt)`$ $``$ $`{\displaystyle e^{i\mathrm{\Omega }(t\frac{z}{v_g})}(\frac{c^2}{v_g^2n^2(\mathrm{\Omega })}1)\frac{a}{\rho }\mathrm{\Omega }𝑑\mathrm{\Omega }}`$ (4)
$`E_\rho (\rho ,zv_gt)`$ $``$ $`{\displaystyle }ie^{i\mathrm{\Omega }(t\frac{z}{v_g})}{\displaystyle \frac{1}{n^2(\mathrm{\Omega })}}{\displaystyle \frac{^2a}{\rho ^2}}d\mathrm{\Omega },`$ (5)
where $`E_z`$ and $`E_\rho `$ are the $`z`$ and $`\rho `$ components of the electric field, and

$$a(\rho ,\mathrm{\Omega })=H_0^{(1)}[s(\mathrm{\Omega })\rho ]$$
(6)
where $`H_0^{(1)}`$ is a Hankel function and $`s(\mathrm{\Omega })=\mathrm{\Omega }(ϵv_g^2/c^21)^{1/2}/v_g=\mathrm{\Omega }\mathrm{tan}\theta _C/v_g`$.

Radiation is emitted at frequencies for which the Cherenkov condition $`v_g>c/n(\mathrm{\Omega })`$ is satisfied. Integrating Eq. (5) over these frequencies, two qualitatively different radiation patterns emerge, depending on $`v_g`$.Stevens et al. (2002) For $`v_g>c/n(0)`$, the laser pulse propagates faster than the entire lower polariton branch $`0<\mathrm{\Omega }<\mathrm{\Omega }_{\mathrm{TO}}`$. We call this the superluminal regime. For $`v_g<c/n(0)`$, the laser pulse only travels faster for frequencies in the band $`\mathrm{\Omega }_C<\mathrm{\Omega }<\mathrm{\Omega }_{\mathrm{TO}}`$, where the velocity matched frequency $`\mathrm{\Omega }_C`$ is defined by $`v_g=c/n(\mathrm{\Omega }_C)`$. We have referred to this as the subluminal regime.gin

Shown in Fig. 2 are plots of $`E_\rho `$, calculated using Eq. (5), in a medium with refractive index given by Eq. (1). For the superluminal regime (Fig. 2a), there is a sharp shock front with small amplitude ripples behind it. The weak ripples are due to dispersion, but the shock wave is defined by $`\mathrm{\Omega }=0`$. On the other hand, the subluminal regime (Fig. 2b) is dominated by dispersion effects. The angle of the phase front depends on $`\rho `$, and there is a large field near the axis of motion $`\rho =0`$. The radiation pattern is still confined within a cone, but the cone angle becomes larger as $`v_g`$ decreases, in contrast to the superluminal regime. The special features of the subluminal regime are apparent in the calculations reported by Afanasiev and co-workers.Afanasiev et al. (1999)

In the CR literature, it is common to use the asymptotic form of the fields for large $`\rho `$ to find an approximate analytical solution. In this limit, $`E_\rho `$ becomes

$$E_\rho _{\mathrm{\Omega }_C}^{\mathrm{\Omega }_{\mathrm{TO}}}\frac{i\sqrt{2s^3(\mathrm{\Omega })}}{n^2(\mathrm{\Omega })\sqrt{\pi \rho }}e^{i[\mathrm{\Omega }(tz/v_g)s(\mathrm{\Omega })\rho +\pi /4]}𝑑\mathrm{\Omega }.$$
(7)
For large $`\rho `$, the integral is small except for values of $`\mathrm{\Omega }`$ for which the integrand does not oscillate. Using the stationary phase method, Eq. 7 can be approximated as

$$E_\rho \underset{\mathrm{\Omega }_i}{}\frac{2\sqrt{s^3(\mathrm{\Omega }_i)}}{n^2(\mathrm{\Omega }_i)\rho }\frac{1}{\sqrt{|s^{\prime \prime }(\mathrm{\Omega }_i)|}}e^{i[\mathrm{\Omega }_i(tz/v_g)s(\mathrm{\Omega }_i)\rho ]},$$
(8)
where $`\mathrm{\Omega }_i`$ are the frequencies within the range of integration for which the phase in the exponential is stationary. These are the solutions to

$$v_g\frac{ds}{d\mathrm{\Omega }}=\frac{v_gtz}{\rho }=\mathrm{cot}\theta ,$$
(9)
which defines a cone of angle $`\theta `$ with respect to the $`z`$ axis. The cone angle depends on frequency because $`s(\mathrm{\Omega })`$ contains $`n(\mathrm{\Omega })`$. Explicitly, Eq. (9) becomes

$$v_g\frac{ds}{d\mathrm{\Omega }}=\mathrm{\Omega }^2\frac{v_g}{sc^2}\frac{dϵ}{d\mathrm{\Omega }}+\frac{v_gs}{\mathrm{\Omega }}=\mathrm{cot}\theta .$$
(10)
The dependence of $`\theta `$ on $`\mathrm{\Omega }`$ is plotted for the two regimes in Fig. 3. In the superluminal case, $`\theta `$ is a monotonic function of frequency, whereas in the subluminal regime two frequencies map onto the same angle, resulting in the complex beating behavior shown in Fig. 2b.

The angle $`\theta _C`$, defined as the smallest angle at which radiation is emitted, is shown as a dashed line in Fig. 3. The frequencies associated with this angle are solutions to $`s^{\prime \prime }(\mathrm{\Omega })=0`$. In the superluminal regime, $`ϵ(\mathrm{\Omega })v_g^2/c^21>0`$ for all frequencies within the range of integration, and $`\mathrm{\Omega }=0`$ is the only solution. Thus, from Eq. (10), $`\mathrm{cos}\theta _C=v_gn(0)/c`$, which is the standard Cherenkov formula. In the subluminal regime, there are two solutions. For a material with dispersion given by Eq. (1), $`\theta _C=\mathrm{tan}^1(4\gamma ^1\sqrt{(\eta \zeta )/\zeta ^3})`$, where $`\gamma =(1ϵ_{\mathrm{}}v_g^2/c^2)^{1/2}`$, $`\zeta =2(43\eta )^{1/2}`$, and $`\eta =1(\mathrm{\Omega }_C/\mathrm{\Omega }_{\mathrm{TO}})^2`$.Afanasiev et al. (1999)

So far, we have considered isotropic materials, and therefore, the considerations above apply to ZnTe and GaP but, in the strictest sense, not to the uniaxial LiTaO3. The theory of CR in an anisotropic material is considerably more complex than that discussed here, because extraordinary and ordinary waves may be emitted. Other complications result for the case of an extremely tightly focused pulse, for which there exists a small longitudinal component of the laser electric field, leading to the breakdown of the point dipole approximation. In addition, the confocal parameter comes into play resulting in the breakdown of the infinite crystal approximation. We note in passing that exact integrals for an electric dipole oriented arbitrarily with respect to the axis of motion, as well as those of a magnetic dipole, are given in the work of Zrelov.Zrelov (1970)

### III.2 Planar source

When the beam waist $`w`$ is much larger than the wavelength of the polariton ($`w>>2\pi cn/\mathrm{\Omega }`$), we have another important limiting case. Here, the polarization induced by the optical pulse is shaped like a pancake. This planar geometry is often used for the generation and detection of THz pulses by nonresonant optical rectification.Wu and Zhang (1996); Nahata et al. (1996) When applied to this problem, the CR picture is useful because it quantifies phase matching issues encountered in thick crystals due to dispersion in the infrared.

For an infinitely thick crystal, and approximating the optical pulse as infinitely short (this is valid as long as the pulse is short compared to the period of the radiation being generated) the polarization induced is $`𝐏\zeta \delta (zv_gt)𝐞_d`$, where $`𝐞_d`$ is a unit vector pointing in the direction of the dipole orientation and $`\zeta `$ is the areal polarization. It is possible to obtain the CR field from the point dipole solution by convolution, but here we use the corresponding Hertz potential

$$𝚷(\mathrm{\Omega })=\frac{4\pi \zeta e^{i\mathrm{\Omega }(tz/v_g)}}{v_gϵ(\mathrm{\Omega })s^2(\mathrm{\Omega })}𝐞_d.$$
(11)
The electric field is thenStevens et al. (2001)

$`𝐄`$ $`=`$ $`{\displaystyle _0^{\mathrm{\Omega }_{\mathrm{TO}}}}{\displaystyle \frac{ϵ(\mathrm{\Omega })\mathrm{\Omega }^2}{c^2}}𝚷(\mathrm{\Omega })𝑑\mathrm{\Omega }`$ (12)
$`=`$ $`{\displaystyle \frac{4\pi \zeta }{v_gc^2}}𝐞_d{\displaystyle _0^{\mathrm{\Omega }_{\mathrm{TO}}}}{\displaystyle \frac{\mathrm{\Omega }^2e^{i\mathrm{\Omega }(tz/v_g)}}{s^2(\mathrm{\Omega })}}𝑑\mathrm{\Omega }.`$ (13)
It can be shown that the field vanishes unless $`s^2(\mathrm{\Omega })=0`$. In the superluminal regime, $`s^2(\mathrm{\Omega })>0`$ for all $`\mathrm{\Omega }`$, and $`𝐄=0`$. In the subluminal regime, $`s^2(\mathrm{\Omega }_C)=0`$. Note that the phase matching condition is *only* met in the presence of dispersion. In the superluminal regime, no radiation is emitted from the bulk because the pulse is traveling faster than the phase velocities of all infrared frequencies.

Practically, a beam cannot be focused to a point, nor is it really planar. If cylindrical lenses are used to focus the beam, it can be extended in one direction and tightly focused in another. The extreme case of the 2D geometry was treated by Kleinman and Auston.Kleinman and Auston (1984) The point geometry solution is the Green’s function of the problem, in that the radiation pattern of an arbitrary source can be obtained by convolution of the point solution with the source shape. Two examples are shown in Fig. 4. This is a powerful machinery for calculating the radiation field generated by ultrafast optical pulses which have been shaped spatially and temporally, for example, in the manner recently demonstrated by Koehl, Adachi, and Nelson.Koehl et al. (1999) To the best of our knowledge, our approach has not been previously applied to this problem.

### III.3 Periodic source (transient grating)

The last source shape we consider is a periodic distribution, or transient grating, implemented by interfering two pump beams inside the sample. This geometry is useful because it excites a unique polariton wavevector determined by the grating period. If we cross two pump beams of equal intensity at an angle $`\alpha `$, Gaussian in space with waists $`w_x`$ in the $`x`$ direction and $`w_y`$ in the $`y`$ direction and Gaussian in time with pulsewidth $`\tau `$, the intensity grating has the form

$$I(x,y,zv_gt)e^{x^2/w_x^2+y^2/w_y^2}\mathrm{cos}(\frac{k_Lx}{\mathrm{sin}\alpha })e^{(z/v_gt)^2/\tau ^2}$$
(14)
where $`k_L`$ is the wavevector of the laser light inside the medium. Here, we have assumed that $`\alpha `$ is small, which makes $`v_g`$ approximately equal to the group velocity of each separate pulse. This is a good approximation for the experiments described here, since $`k_L`$ is much larger than the wavevector of the gratings. The Cherenkov angle converts the spatial periodicity of the transient grating into a temporal periodicity in the polariton field, resulting in a checkerboard pattern and a well-defined polariton frequency $`\mathrm{\Omega }`$, as shown in Fig. 5.

In the impulsive limit, a defocused probe pulse traveling at the same speed as the pump pulse and, thus, at the speed of the polariton field, samples a slice of the checkerboard field in the $`xy`$ plane. The electric field in this slice is periodic with the same wavevector $`k_x=k_L/\mathrm{sin}\alpha `$ as the intensity grating. Through the electro-optic effect, the probe pulse diffracts from the periodically varying index of refraction, resulting in an oscillatory signal at $`2\mathrm{\Omega }`$ as a function of time delay, because the diffracted intensity is proportional to $`|E|^2`$. This technique is known as homodyne detection in the ISRS literature.Wiederrecht et al. (1995) Heterodyne detection,Bakker et al. (1998a); Gautier et al. (1996); Crimmins et al. (2002) in which one mixes the diffracted signal with undiffracted probe light, gives a much larger signal at $`\mathrm{\Omega }`$, because the signal is proportional to $`E`$. Because it is sensitive to the phase of the polariton field, which defines the phase of the diffracted probe pulse, heterodyne detection works just as well outside the checkerboard, in what we might call the Cherenkov wing of the field. As we discuss in detail below, heterodyne detection is often discussed in terms of four-wave mixing.Gautier et al. (1996)

It is important to emphasize that the wavevector of the grating is not the same as the polariton wavevector, since the polariton propagates in the direction normal to the wavefront. The magnitude of the polariton wavevector is given by $`q=k_x\mathrm{sin}\theta `$, where $`\theta `$ is the Cherenkov angle. In many early ISRS papers, it is assumed that the wavevector of the grating is the same as the polariton wavevector, tantamount to assuming that $`\theta `$ is very large. This is generally a poor approximation, and leads to incorrect polariton dispersion measurements, as recently pointed out by Hebling.Hebling (2002)

The simple interpretation above works best for large $`w`$, in which case the polariton wavevector is well-defined. For small $`w_x`$, edge effects make the polariton field more complex but, by convolving the point dipole solution with the grating (as in Fig. 5), we can calculate the polariton field due to a grating of arbitrary shape.

## IV Phase matching

One source of confusion when comparing the transient grating, the tightly focused, and the planar geometries is that, at first glance, the phase matching conditions seem entirely different. Here, we clarify this point. An ultrafast pulse contains a continuum of frequencies and wavevectors, and the generated polariton wavevectors are determined by conservation of energy and momentum. The Cherenkov interpretation is, of course, consistent with phase matching. In all geometries, the nonlinear polarization contains terms of the form $`e^{i(𝐤_2𝐤_1).𝐫i(\omega _2\omega _1)t}`$, where $`𝐤_1`$, $`𝐤_2`$, $`\omega _1`$, and $`\omega _2`$ denote wavevector and frequency components of the optical pulse. For the planar geometry, the wavevector components all point along the $`z`$ axis, and conservation of momentum requires $`q=k_2k_1=[\omega _2n(\omega _2)\omega _1n(\omega _1)]/c`$, where $`q`$ is the polariton wavevector. Since $`\mathrm{\Omega }=\omega _1\omega _2`$,

$`q`$ $`=`$ $`{\displaystyle \frac{\omega _2n(\omega _2)(\omega _2+\mathrm{\Omega })n(\omega _2+\mathrm{\Omega })}{c}}`$ (15)
$``$ $`{\displaystyle \frac{\mathrm{\Omega }}{c}}[n(\omega _2)+\omega _2n^{}(\omega _2)]={\displaystyle \frac{\mathrm{\Omega }}{v_g(\omega _2)}},`$ (16)
making the connection to the group velocity matching argument derived in Section IIIB.

A tightly focused optical pulse has wavevector components pointing in many directions. The dependence of the polarization on $`zv_gt`$ forces $`q_z=\mathrm{\Omega }/v_g`$. The polariton wavevector and frequency must also satisfy $`\mathrm{\Omega }=|𝐪|c/n(\mathrm{\Omega })`$. Let $`q_\rho `$ be the component of the polariton wavevector perpendicular to the $`z`$ axis. Momentum conservation gives

$$\frac{q_\rho }{q_z}=\frac{\sqrt{q^2q_z^2}}{\mathrm{\Omega }/v_g}=\sqrt{\frac{n^2(\mathrm{\Omega })v_g^2}{c^2}1},$$
(17)
and the polariton is emitted at an angle $`\mathrm{tan}^1q_\rho /q_z`$ to the $`z`$ axis, which is identical to the Cherenkov angle.

In the transient grating geometry, the angle between $`𝐤_1`$ and $`𝐤_2`$ is well-defined, and the wavevector $`𝐪=𝐤_1𝐤_2`$ is emitted in the direction given by Eq. (17) because the polarization is a function of $`zv_gt`$. The transient grating configuration allows one to excite wavevectors that are higher than those accessible in the planar geometry.

A recently proposed technique developed by Hebling and co-workers relies on a tilted pulse front.Hebling et al. (2002) This can also be encompassed in a Cherenkov radiation framework as a planar source which propagates at an angle. Tuning the tilt angle changes the wavevector of the generated polariton, and this shows promise as an easily tunable THz source.

As is well-known, it is possible to interpret ISRS as four wave mixing.Gautier et al. (1996) Here one hides both the generation and detection of the polariton in $`\chi _{\mathrm{eff}}^{(3)}`$, much as we included the Raman contribution to the generation of the polariton in $`\chi ^{(2)}`$. The four-wave mixing polarization is

$`P(\omega _s)`$ $`=`$ $`\chi _{\mathrm{eff}}^{(3)}(\omega _p,\mathrm{\Omega }\omega _p,\omega _p)E_1(\omega _1)E_2^{}(\omega _1\mathrm{\Omega })E_p(\omega _p)`$ (18)
$`=`$ $`d(\mathrm{\Omega })[d(\mathrm{\Omega })E_1(\omega _1)E_2^{}(\omega _1\mathrm{\Omega })]E_p(\omega _p)`$
where $`E_1`$ and $`E_2`$ are the pump fields, $`E_p`$ is the probe field, and $`d`$ is the effective nonlinear susceptibility defined in Section II. All propagation and damping effects can be included in $`\chi _{\mathrm{eff}}^{(3)}`$, and this formalism is often used.Gale et al. (1986)

Finally, we discuss briefly the advantages and disadvantages of the various scattering geometries. The experiments with tightly focused beams produce the highest peak polariton amplitude (for a given pulse peak power), but also a wide range of wavevectors and frequencies. ISRS by a transient grating has the advantage of defining the polariton wavevector and producing a large amplitude polariton at that particular wavevector. CARS has the advantage of yielding the best frequency and wavevector resolution.

## V Materials

Generation of electro-optic CR requires a non-zero $`\chi ^{(2)}`$, present only in crystals without inversion symmetry. In our experiments, we used three readily available and well-characterized materials: the zincblende ZnTe and GaP, and LiTaO3. The parameters relevant to our work are listed in Table 1.

### V.1 GaP and ZnTe

Materials with the zincblende structure have a triply degenerate infrared-active optical mode. For light normally incident on a (110) face, the selection rules forbid excitation of the nondispersive LO mode. For the TO mode,

$$R_1=\left(\begin{array}{ccc}0& 0& b\\ 0& 0& b\\ b& b& 0\end{array}\right)\text{ and }R_2=\left(\begin{array}{ccc}0& b& 0\\ b& 0& 0\\ 0& 0& 0\end{array}\right),$$
(19)
where $`R_1`$ and $`R_2`$ are the Raman tensors for the modes polarized along the \[$`1\overline{1}0`$\] and  directions, respectively. The single independent tensor component $`b`$ corresponds to the nonlinear coefficient $`d_{41}`$. The dipolar distribution induced by the pump pulse is always oriented perpendicular to the  direction. We measure the polariton by the inverse process, and we are only sensitive to the electric field component $`E_\rho `$.

In their pioneering work using nonlinear mixing of visible and infrared cw laser light, Faust and Henry measured $`C=0.47`$ in GaP.Faust and Henry (1966) This leads to total destructive interference between the ionic and electronic contributions to $`d`$ at 6.5 THz. Because this frequency is close to $`\mathrm{\Omega }_C`$ for Ti:sapphire pulses, this results in a small signal. Recently Leitenstorfer and co-workers found $`C=0.07`$ in ZnTe by analyzing data from time-domain spectroscopy of THz pulses by electro-optic sampling.Leitenstorfer et al. (1999) This small value leads to $`d=0`$ at 5.2 THz, very close to $`\mathrm{\Omega }_{\mathrm{TO}}`$. For our experiment in ZnTe, we are only able to measure up to roughly 4 THz, not high enough to probe effects of the dispersion of $`d`$.

### V.2 LiTaO3

This material crystallizes in the perovskite structure. LiTaO3 is ferroelectric at room temperature and exhibits 4 $`A_1`$ symmetry and 9 $`E`$ symmetry modes, all of which are Raman and infrared active.Raptis (1988) There has been considerable disagreement over the assignment of several Raman peaks, but the lowest-frequency $`A_1`$ mode considered here is well-characterized. The Raman tensor for $`A_1`$ phonons isRaptis (1988)

$$A_1=\left(\begin{array}{ccc}a& 0& 0\\ 0& a& 0\\ 0& 0& b\end{array}\right).$$
(20)
The independent tensor elements $`a`$ and $`b`$ are associated with nonlinear coefficients $`d_{31}`$ and $`d_{33}`$, respectively. Here, we use $`\chi ^{(2)}=d_{33}`$ and $`\chi ^\mathrm{R}=b`$, which are the largest coefficients, to both generate and detect the polariton. The dipolar distribution points along the optic axis, and the probe pulse is sensitive to the component of the Cherenkov electric field along that axis. Barker, Ballman, and Ditzenberger measured the infrared reflectivity and found that, for the extraordinary index, most of the oscillator strength in the infrared is in the lowest-lying $`A_1`$ TO mode.Barker et al. (1970) This mode is also the strongest Raman scatterer. We performed our experiment on the lowest polariton branch.

The refractive index in the visible is weakly anisotropic. Since the laser pulses in our experiment are polarized along the optic axis, there is no effect of this anisotropy, in contrast with experiments on $`E`$ modes.Albert et al. (1996) The static dielectric constant is also weakly anisotropic. However, since the lowest-lying $`A_1`$ and $`E`$ TO frequencies are different (6 and 4.2 THz, respectively), the dielectric function becomes more anisotropic at large frequencies, where the cubic approximation is expected to break down. For the low frequencies discussed here, the approximation is well obeyed. For the lowest-lying $`A_1`$ mode in LiTaO3, $`C`$ is positive.Boyd and Pollack (1973)

## VI Experiments

### VI.1 Techniques

Our experimental technique is similar to that of Auston and Nuss.Auston and Nuss (1988) A diagram of the experimental setup is shown in Fig. 6. Pump and probe beams are focused and crossed inside the sample. The pump pulse is the Cherenkov source and the tightly focused probe serves as the detector of the Cherenkov radiation field. The probe pulse “surfs” on top of the Cherenkov wakefield generated by the pump pulse and samples the polariton field through both the linear electro-optic and the Raman effects, which both cause changes in the refractive index. We change the time delay ($`t`$) between pump and probe with a motorized delay stage in the probe path and the relative focal position ($`\rho `$) between pump and probe by moving a motorized translation stage upon which the pump’s focusing lens is mounted. By sampling the diffracted probe beam as a function of $`t`$ and $`\rho `$, we map the polariton field.

We used two different detection schemes. The first is a frequency shift measurement, which takes advantage of the fact that the spectrum of the probe is shifted proportional to $`dn/dt`$.Dhar et al. (1994) The polariton induces a change in the index of refraction due to its electromagnetic and vibrational components, so that $`dn/dt`$ contains a term proportional to $`dE/dt`$ and another one to $`dQ/dt`$. We measure the frequency shift by cutting half of the probe spectrum with a bandpass filter before detection with a photodiode. This technique is especially useful when the pump and probe polarizations are chosen perpendicular to each other to suppress the scattered pump light. The second technique, widely used for detecting THz radiation, is polarization detection.Valdmanis et al. (1982); Wu and Zhang (1996) When properly oriented, the polariton’s electric field induces birefringence, which can be probed by measuring the polarization state of the probe beam after the sample. This yields a signal directly proportional to the THz Cherenkov electric field.

The spot size inside the crystal was determined experimentally using the same technique. In ZnTe and GaP, two photon absorption depletes the probe when pump and probe overlap temporally and spatially, creating a large Gaussian peak in the signal. In LiTaO3, the intensity-dependent refractive index due to $`\chi ^{(3)}`$ modulates the probe pulse’s phase when it overlaps with the pump pulse, also producing a large signal (known as the coherent artifact).

Auston and Nuss used thin samples and a collinear geometry to achieve spatial resolution on the order of 1.5 $`\mu `$m.Auston and Nuss (1988) A problem with this setup is that it is difficult to measure the field behind the source, because of strong pump scattering due to the fact that the two beams are collinear. In order to measure the source itself plus the field behind it, we used the nearly-collinear setup described above, which enables the spatial separation of the pump and probe beams even when they overlap spatially inside the sample. This introduces a source of distortion in the data due to spatial walk-off between the pump and probe pulses as they travel through the sample. To minimize this walk-off, we focused the pump and probe beams with two small lenses placed as close together as possible. In order to ensure that the walk-off did not affect the spatial resolution, we made the plane shared by the pump and probe beams perpendicular to the axis along which the pump lens was moved. We accomplished this by directing the probe beam below the pump lens, and moving the relative focal position horizontally. The resolution in $`zvt`$ is determined by the pulsewidth of the probe pulse, and in $`\rho `$ by the focal waist of the probe beam.

Ultrafast pulses were generated by a Ti:sapphire oscillator producing 800 nm, 60 fs pulses at the repetition rate of 82 MHz. Acousto-optic amplitude modulators were used to modulate the pump and probe beams at 3 MHz and 2.95 MHz, respectively, and the probe was detected with an unbiased photodiode at the difference frequency, 50 kHz, with a lock-in amplifier. The maximum pump and probe average power at the sample were $``$140 and 30 mW, respectively, limited mostly by the efficiency of the modulators. We emphasize the importance of the double-modulation technique because it helps reject the scattered pump light. Recently, experiments have been performed using a Zernike method with a CCD camera, to image the phase change of a weakly-focused probe beam caused by the polariton field as a function of position and time delay.Koehl et al. (1999); Koehl and Nelson (2001) This is another powerful technique for imaging the Cherenkov field, although it suffers from unwanted large scattering near the pump pulse.

### VI.2 Point source: Experiments

If the pump pulse is focused to a waist smaller than the wavelength of the polariton, we can approximate it as a point dipole. Because of the large thickness of our samples, we are limited to a focal waist of about 20 $`\mu `$m which corresponds to a maximum polariton frequency of 2-4 THz, depending on the material. For lower frequencies, the pump pulse acts like a point dipole.

First we discuss the data on LiTaO3, a superluminal material; see Fig. 7. Our sample was an $`x`$-cut wafer 2 inches in diameter and 1 mm thick, and we used frequency shift detection. For pulses of wavelength 800 nm, $`\theta _C=69^{}`$. The experiment matches the calculations quite well after taking into account the 20 $`\mu `$m pump waist by convoluting the point dipole solution with a Gaussian of that width. However, the absence of ripples behind the sharp Cherenkov cone indicates that damping, not included in the simulation, is strong for higher frequencies. In addition to polaritons of the lowest branch, we were also able to generate and detect polaritons and coherent phonons from higher branches; see later.

Our measurements on ZnTe, a subluminal material, and the comparison with calculations are shown in Fig. 8. The (110) oriented ZnTe sample was 10$`\times `$10$`\times `$1 mm3, and polarization-sensitive detection was used. The data shows many of the features predicted for subluminal CR. The artifact that appears about 1.6 ps after the pump pulse is due to a second pulse about 1000 times smaller than the pump pulse, due to a reflection from an optical element upstream. We have subracted the contribution from an exponentially decaying signal caused by free carriers generated by two-photon absorption. This signal is partially balanced out by using polarization-sensitive detection, and the remnant has been fitted and subtracted out in Fig. 8a to show the absence of CR near $`\rho =0`$.

As discussed earlier, subluminal CR is expected to exhibit a beating pattern. In Fig. 8a, there is a striking node that separates the shock wave into two and possibly three distinct parts. This can be explained in the calculations by both subluminal beats, or by a convolution artifact due to the finite spatial size of the pump and probe pulses. However, note that the node is not as prominent in the calculation. This discrepancy, and the lack of experimental signal near $`\rho =0`$, is probably due to the finite size of the crystal, as discussed in detail later.

### VI.3 Planar source: Experiments

For the planar geometry, scans were performed only along the $`\rho `$ axis in the $`zv_gt`$ direction. Results for GaP and ZnTe with a pump waist of roughly 100 $`\mu `$m are shown in Fig. 9. The GaP sample, (110) oriented, was a 4$`\times `$1$`\times `$1 mm3 crystal. The time scans show oscillations at the group velocity matched frequencies predicted by theory, with an apparent decay probably due to a combination of the finite length of the sample (see below), and the finite extent of the pump pulse.Stevens et al. (2001) The planar geometry is commonly used for generating THz pulses by nonresonant optical rectification.Nahata et al. (1996) The distortion of the electro-optic sampling signal due to infrared dispersion, especially in a thick crystal, has been the subject of much discussion recently.Bakker et al. (1998b); Gallot and Grischkowsky (1999)

We also performed an experiment on the superluminal material LiTaO3. As expected, there was no signal from the lowest polariton branch. However, our pulses were short enough to see higher branches. The phase matching condition is given by $`\mathrm{\Omega }/v_g`$, and this expression was found to match the higher polariton branches quite well.

### VI.4 Finite crystal effects

Point geometry experiments were also performed on a 300 $`\mu `$m thick ZnTe sample. As is apparent in Fig. 8c, the node we observed in the 1-mm-thick crystal has now disappeared. This is probably a result of the breakdown of the infinite crystal approximation. In the superluminal regime, the infinite crystal calculation works well because the radiation is emitted at large angles at all frequencies down to zero. Therefore, effects of the finiteness of the crystal can be expected only at the very edge of the shock wave, at $`\rho =L\mathrm{tan}\theta _C=Ln_0/c`$, where $`L`$ is the length of the crystal. For example, the edge of the shock wave is at $`\rho 2.6`$ mm in our 1 mm thick LiTaO3 sample, well beyond what we can practically measure due to the fall-off proportional to $`\rho ^{1/2}`$ in the electric field from a point source.

The situation is considerably different in the subluminal regime because polaritons near $`\mathrm{\Omega }_C`$ are emitted in the near forward direction. In the 1 mm thick ZnTe crystal, the Cherenkov angle for 2.5 THz polaritons is about $`9^{}`$, so that this component reaches $`\rho =160`$ $`\mu `$m before the pulse reaches the end of the crystal. In the 300 $`\mu `$m thick crystal, it only reaches $`\rho =50`$ $`\mu `$m. Thus, the subluminal beats, due to interference between frequency components sent in the forward direction, require some propagation time to develop. We believe this is the reason why they do not appear in the thin crystal. This could also explain other discrepancies between theory and experiment in Fig. 8. The large oscillatory signal at $`\mathrm{\Omega }_C`$ near $`\rho =0`$ in the calculation is clearly absent in the experiment, probably because radiation at $`\mathrm{\Omega }_C`$ is emitted in the forward direction The difference in the shape of the cone between theory and experiment is likely due to the same effect. A finite crystal calculation is probably needed to fully understand the details of the experiment. Note that CR in a finite crystal is similar to the so-called “Tamm problem” of CR for a particle traveling along a finite path.Tamm (1939); Afanasiev et al. (2000)

The finiteness of the crystal is expected to play an even more important role in the planar geometry, in which Cherenkov components are emitted exclusively in the forward direction. The CR pattern is the result of interference among these components, so that this pattern should change dramatically with crystal thickness. A very thin crystal would show no propagation effects, simply reproducing the time derivative of the pump pulse envelope as a short THz pulse.

### VI.5 Periodic source: Experiments

We adapted the setup described previously to perform experiments with two pump pulses, creating a periodic dipolar distribution in LiTaO3, and used frequency shift detection to map the polariton field. The pump pulse was split into two equally intense pulses which were focused by a single lens and crossed inside the sample at an angle $`\alpha `$. Using a beam waist such that $`w<<v_g\tau /\mathrm{sin}\alpha `$, we avoided the so-called “pancake effect,”Crimmins et al. (2002) a shrinking of the overlap between weakly-focused pulses due to their short duration. In our experiments, we did not have to measure $`\alpha `$, because the grating period was directly obtained by measuring the coherent artifact as a function of spatial position. As for the point geometry, we mapped the polariton field due to this pump configuration as a function of $`zv_gt`$ and $`\rho `$. The results in Fig. 10 show excellent agreement with the theoretical predictions.

### VI.6 LiTaO3: Phonons and Polaritons

LiTaO3 has a particularly rich history in the ISRS literature.Bakker et al. (1992, 1993); Wiederrecht et al. (1995); Brennan and Nelson (1997); Perry and Dougherty (1997); Bakker et al. (1998a); Bakker (1998); Romero-Rochin et al. (1999); Koehl et al. (1999); Koehl and Nelson (2001); Crimmins et al. (2002) In this section we compare our data to those of other reports. There are two main points of contention in the literature: the origin of frequency-dependent damping effects and avoided crossings at low frequency. Our results for the polariton dispersion are summarized in Fig. 11.

As previously discussed for the transient grating geometry, see below Eq. (17), the Cherenkov angle appears in the expression for the polariton wavevector. Using the correct expression, our data fall almost exactly on the dispersion relation of the lowest branch $`A_1`$ polaritons inferred from the IR data of Barker et al.Barker et al. (1970) Recent ISRS papers have begun to take this into account as well.Hebling (2002); Crimmins et al. (2002) The grating technique is also sensitive to non-dispersive phonons, and the signal at 7.8 THz roughly matches an $`A_1`$ phonon peak seen in Raman spectra.Raptis (1988) Scans in the planar geometry, with wavevector given by the group velocity of the pulse, are also shown in Fig. 11. Here, we are sensitive to even more phonon branches, including a higher-branch polariton above 12 THz. Note that higher branches were also observed in a recent heterodyne experiment.Crimmins et al. (2002)

Our measurements of the polariton damping rate $`\gamma `$ by the grating method,dam  agree very well with those reported by others: $`\gamma =2.4`$ cm-1 at 1 THz, $`\gamma =4.4`$ cm-1 at 2 THz, and $`\gamma =7.5`$ cm-1 at 3 THz. Large damping, peaked at 2.7 THz, was reported by Auston and Nuss in their original experiment.Auston and Nuss (1988) They suggested that this was due to multi-phonon absorption. In their homodyne transient grating experiment, Wiederrecht et al Wiederrecht et al. (1995) measured increased damping peaked at 2.76 THz and 4.38 THz, and attributed this to strain-induced coupling to heavily damped $`E`$-symmetry modes at the same frequencies. Very recently, Crimmins, Stoyanov, and Nelson measured the damping rate using heterodyne detection and suggested that the damping is due to quadratic coupling between the phonon-polariton mode at 2.1 THz and the lowest-lying $`E`$-mode at 4.2 THz.Crimmins et al. (2002)

The inelastic neutron scattering data of Cheng et alCheng et al. (1993) supports the idea that the damping is due to anharmonicity, more precisely, to a process involving the simultaneous creation and destruction of two different phonons. Close to the observed peaks at 2.76 and 4.38 THz, the phonon dispersion shows TO and LA (TA) branches for which the frequency difference at the zone edge, where the density of states is large, is 2.6 THz (4 THz). In this context, we note that decay by two-phonon difference modes has been seen in many materials, including ZnTe.Schall et al. (2001) Furthermore, many spontaneous Raman studiesPenna et al. (1976); Raptis (1988), including one on our particular crystal,Wahlstrand and Merlin  show a broad peak at 2.55 THz, consistent with a second-order feature. An alternative explanation for polariton broadening, supported by the fact that the TA and LA branches at the zone edge are at 3.2 and 4.6 THz, respectively,Cheng et al. (1993) is coupling to high-wavevector phonons due to disorder-induced breakdown of momentum conservation. Given that LiTaO3 is often non-stoichiometric, defect-related processes may contribute strongly to the damping.

Another anomaly reported in LiTaO3 is an extra avoided crossing in the polariton dispersion at 0.9 THz.Bakker et al. (1992, 1993, 1998a) Similar avoided crossings were also reported for LiNbO3,Bakker et al. (1994) and increased absorption was observed at the corresponding frequencies in stimulated Raman gain measurements.Qiu et al. (1995) The avoided crossings were assigned by Bakker and co-workersBakker et al. (1998a) to quantum beats between levels of an intrinsic anharmonic potential. However, our work on LiTaO3, as well as a recent heterodyne ISRS experiment,Crimmins et al. (2002) show no sign of an avoided crossing at 0.9 THz. The observation of several low-lying modes in a CARS study of MgO-doped LiNbO3Schwarz and Maier (1998) suggests that the extra modes are due to defects. If this were the case, the reported discrepancies for LiTaO3 may simply be due to differences in the concentration of defects for different samples.

We take issue with the intrinsic anharmonicity proposed by Bakker et alBakker et al. (1998a) to explain these extra modes. Unlike in molecules, lattice anharmonicity in solids leads to decay, because anharmonic coupling involves a continuum of modes. Hence, unless the anharmonicity is strong enough to support 2-phonon bound states, it cannot introduce new frequencies. To the best of our knowledge, no such bound states have been observed in LiTaO3 or LiNbO3.

An unexplored area is the role of anisotropy in LiTaO3. As briefly discussed earlier, for frequencies approaching the lowest $`E`$-phonon frequency, the index of refraction becomes extremely anisotropic. Anisotropy is known to have an important effect on CRZrelov (1970) but as far as we know, it has been completely ignored in the ISRS literature. For a point charge traveling through a uniaxial material in a direction perpendicular to the optical axis, both ordinary and extraordinary waves are emitted, with the extraordinary waves polarized at a finite angle to the optical axis. Assuming this is also true for a dipole oriented along the optic axis, it should have a large effect on the CR pattern.

Finally, we note the observation of wavevector overtones in transient grating experiments using amplified pulses, which Nelson et alBrennan and Nelson (1997); Romero-Rochin et al. (1999) ascribed to lattice anharmonicity. Since polariton nonlinearity is primarily due to phonon anharmonicity and the polaritons in question are mainly light-like, we consider this interpretation to be questionable. Instead, we propose that the higher-order diffraction peaks simply reflect a non-sinusoidal grating due to higher-order optical nonlinearities. This is an area which needs to be explored further.

## VII Conclusions

In summary, we have used ultrafast optical pulses to generate and image coherent polariton fields in ZnTe, GaP and LiTaO3. For single and multiple pulses of arbitrary shape and, in particular, for the transient-grating geometry, we have shown that the polariton field can be calculated by convolving the slowly-varying pulse envelope with the Cherenkov field due to a dipole. The standard phase matching arguments can also be explained by applying ideas from the theory of Cherenkov radiation. Results for point-like and planar sources reveal important differences between the superluminal and subluminal regimes, especially for the Cherenkov emission by a plane of dipoles which can only occur in the subluminal case.

###### Acknowledgements.

We thank G. Narayanasamy for assistance in the early stages of the experiments and T. E. Stevens for useful discussions. We are grateful to M. DeCamp, C. Herne, and J. Murray for providing the ZnTe samples, and J. Deibel for the LiTaO3 sample. This work was supported by the AFOSR under contract F49620-00-1-0328 through the MURI program, and by the NSF FOCUS Physics Frontier Center.