text:
# Atomic wave packet basis for quantum information

## Abstract

We propose a wave packet basis for storing and processing several qubits of quantum information in a single multilevel atom. Using radially localized wave packet states in the Rydberg atom, we construct an orthogonal basis that is related to the usual energy level basis by a quantum Fourier transform. A transform-limited laser pulse that is short compared with the classical Kepler period of the system interacts mainly with the wave packet state localized near the atomic core, allowing selective control in this basis. We argue that wave packet control in this regime is useful for multilevel quantum information processing.

PACS numbers: 03.67.-a, 32.80.Qk

Quantum information processing is a growing interdisciplinary field that has seen rapid progress in the last few years . However physical implementations of large scale quantum networks are severely limited by decoherence , the loss of coherence in the superposition state of the system. This is especially true of extended systems that involve entanglement of macroscopic degrees of freedom, such as the center-of-mass motion of atoms in a trap .

We consider multiple computational levels in each atom to store information more densely. The radiative lifetime of Rydberg levels scales as $`n^3`$, where $`n`$ is the principal quantum number, and approaches a millisecond for $`n>100`$. In most atomic quantum logic schemes, the bottleneck for the coherence time is the macroscopic entanglement of two or more atoms. The multilevel approach minimizes this decoherence by reducing the number of atoms in the quantum network, with each $`d`$-level atom representing $`\text{log}_2d`$ qubits of information . In a $`d=8`$ level scheme, for example, the number of atoms that need to be entangled is less by a factor of three.

A great deal of research has been carried out in the past few years developing techniques for controlling multilevel Rydberg wave packets using shaped short pulses. It is attractive to exploit this work by developing an alternative wave packet basis for the quantum information stored in the atomic energy levels. Rydberg wave packets are coherent superpositions of atomic energy eigenstates with large principal quantum numbers. We can use the relative phases of the different eigenstates in the superposition to store information. Using the energy-time Fourier kernel in unitary time evolution, we define an orthogonal, discrete wave packet basis that is related to the energy levels by a quantum Fourier transform. This transform plays a key role in quantum computing applications , and we show how it naturally permits complementary bases for representing the information stored in a multilevel quantum system.

In the optical domain, a discrete Fourier transform connects the spectral and temporal modes of the radiation field in a cavity. Iaconis et al.  use this

approach to characterize the Husimi distribution of a multimode quantum field in heterodyne detection. Tittel et al.  exploit the complementarity of energy and time in their scheme for quantum cryptography using entangled photons. In the atomic case, Muller and Noordam  introduce a Fourier relation between the energy quantum number and a continuous phase coordinate to study the photoexcitation of Rydberg states. We show how the quantum Fourier transform of the energy levels corresponds exactly to a discrete temporal basis of wave packets in cyclic motion about the nucleus.

We consider radial wave packets, which are superpositions of Rydberg states with low angular momentum. These states are radially localized in probability and evolve periodically for short times, or in the classical correspondence limit . The wave packet basis we introduce corresponds to discrete times in this classical evolution, when each wave packet state is centered at a different radius from the atomic nucleus. This basis is non-stationary in time, with wave function amplitudes evolving cyclically during a classical period. Due to the unequal spacing of the energy levels, the wave function undergoes dispersion, which leads to revivals and super revivals in the wave packet amplitudes. The classical motion and quantum revival structure of this free time evolution allows for the possibility of conveniently controlling the wave function at suitably chosen intervals of time.

We propose using short optical pulses to couple coherently the Rydberg manifold to a low-lying energy level in the atom. A transform-limited pulse that is short compared to the classical Kepler period of the system interacts with all the energy levels but couples only to the wave packet state that is nearest the atomic core, a phenomenon familiar in the excitation and photoionization of Rydberg states . The classical electron has the largest momentum, and hence absorbs energy most efficiently, when near the core. During the pulse, the wave packet state nearest the nucleus and the ground state form an effectively two-level system, and undergo Rabi oscillations.

Coherent control of wave packets based on spatial localization has been studied by Ahn et al. , who show that half-cycle pulses can be used to extract binary phase information stored in a radial wave packet. The use of optical pulses shorter than a Kepler period for quantum control of atomic wave packets in the strong-field regime has been studied analytically by Araujo et al. . A weak field analysis by Noel and Stroud  showed that a uniform train of $`d`$ pulses within one Kepler period can be used to create an arbitrary Rydberg state of $`d`$ levels whose amplitudes are related by a discrete Fourier transform to the pulse amplitudes. We propose to use such pulses for selective processing of quantum information in the wave packet basis.

Consider $`d`$ Rydberg energy levels with angular momentum $`l=1`$ in a hydrogen or alkali atom as a quantum information basis,

$`|j_\nu =|\overline{n}+j,1,0,`$ (1)
$`j=d/2+1,d/2+2,\mathrm{},d/2,`$ (2)
where $`\overline{n}`$ is the mean value for the principal quantum number in the wave packet superposition, and the subscript $`\nu `$ is used to denote a state in the energy basis. We have assumed above that the number of levels $`d`$ is an even number, but the arguments below are easily extended to odd $`d`$. These levels can be excited simultaneously from the ground state by a strong linearly polarized laser pulse with sufficient bandwidth to overlap all the levels. The pulsed excitation creates a wave packet which evolves in time as

$$|\psi (t)=\underset{n}{}b_j\mathrm{exp}(i\omega _jt)|j_\nu ,$$
(3)
where $`\mathrm{}\omega _j`$ is the energy of the $`j`$th Rydberg level and the amplitudes $`b_j`$ are determined by the spectrum of the exciting pulse. If we assume that the pulse spectrum, centered at $`\omega _0`$, is uniform over the $`d`$ computational levels and falls off rapidly outside this manifold, we have a uniform amplitude distribution, $`b_j=1/\sqrt{d}`$, which is time independent in the absence of external fields.

The radial wave packet was first studied by Parker and Stroud , who showed that it corresponds to a radially localized shell of probability distribution moving periodically between the classical turning points for short times. Eventually, the uneven spacing in the energies of the levels in Eq. (1) leads to dispersions and revivals of the wave packet, corresponding to the time-dependent phases in Eq. (3) interfering destructively or constructively. The rates at which the phases evolve with respect to each other depend on the frequencies $`\omega _j`$ relative to the average frequency $`\omega _0`$, which can be expanded in a Taylor series in $`j=n\overline{n}`$,

$`\omega _{j0}`$ $`=`$ $`\omega _j\omega _0={\displaystyle \frac{1}{2(\overline{n}+j)^2}}+{\displaystyle \frac{1}{2\overline{n}^2}}`$ (4)
$`=`$ $`2\pi \left[{\displaystyle \frac{j}{T_K}}{\displaystyle \frac{j^2}{T_{\mathrm{rev}}}}+{\displaystyle \frac{j^3}{T_{\mathrm{sr}}}}\mathrm{}\right],`$ (5)
where in atomic units, $`T_K=2\pi \overline{n}^3`$ is the classical Kepler period for round trip motion, $`T_{\mathrm{rev}}/2!=2\pi \overline{n}^4/3`$ is when the wave packet revives first, $`T_{\mathrm{sr}}/3!=\pi \overline{n}^5/6`$ is when the wave packet super-revives first, etc. At $`\overline{n}=180`$ for example, $`T_K0.89`$ ns, $`T_{\mathrm{rev}}106`$ ns, and $`T_{\mathrm{sr}}14`$ $`\mu `$s. For large $`\overline{n}`$, the different time scales are well separated in magnitude, and for a small number of levels $`d`$, the wave packet is more spread out in space and disperses more slowly in time. For $`d^2\overline{n}`$, we see from Eq. (4) that the revival dynamics of the wave packet can be neglected to good approximation during the first few Kepler periods, and Eq. (3) becomes

$`|\psi (t)`$ $``$ $`{\displaystyle \frac{1}{\sqrt{d}}}{\displaystyle \underset{j}{}}\mathrm{exp}\left(i2\pi jt/T_K\right)|j_\nu ,`$ (6)
$`=`$ $`|\psi _K(t)`$ (7)
which evolves periodically in time with period $`T_K`$.

The Fourier kernel in Eq. (6) involves $`j`$ and $`t`$, and can be thought of as underlying an uncertainty relation between energy and time, with more levels leading to greater spatio-temporal localization of the wave packet. Although the Rydberg energies are discrete, the time $`t`$ for the wave packet varies continuously around the orbit. However, if we consider $`d`$ discrete times during any Kepler period,

$`t_k=k(T_K/d),`$ (8)
$`k=d/2+1,d/2+2,\mathrm{},d/2,`$ (9)
we can construct an alternate discrete wave packet basis $`|k_\tau `$ for the system of energy levels in Eq. (1) based on the quantum Fourier transform (QFT),

$`|k_\tau `$ $`=`$ $`|\psi _K(t_k)`$ (10)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{d}}}{\displaystyle \underset{j}{}}\mathrm{exp}\left(i2\pi jk/d\right)|j_\nu .`$ (11)
This is an orthogonal basis in time (denoted by the subscript $`\tau `$), and the basis amplitudes evolve periodically with the classical period $`T_K`$, as illustrated in Fig. 1. Thus $`k=0`$ and $`k=d/2`$ correspond to wave packet states centered at the inner and outer turning points of the classical orbit respectively. The sign of $`k`$ indicates the direction of propagation of the wave packet. We can also define the basis for odd $`d`$, except that in this case $`k=\pm (d1)/2`$ correspond to wave packet states centered near (but not exactly at) the outer turning point, having nonzero and opposite momenta. The energy eigenstates can be recovered from $`|k_\tau `$ by an inverse QFT,

$$|j_\nu =\frac{1}{\sqrt{d}}\underset{k}{}\mathrm{exp}\left(i2\pi jk/d\right)|k_\tau .$$
(12)
In general, any excited state $`|\mathrm{\Psi }(t)`$ of the $`d`$ Rydberg levels can be expanded in either basis,

$`|\mathrm{\Psi }(t)`$ $`=`$ $`{\displaystyle \underset{j}{}}b_j(t)\mathrm{exp}(i\omega _jt)|j_\nu `$ (13)
$`=`$ $`\mathrm{exp}(i\omega _0t){\displaystyle \underset{k}{}}\stackrel{~}{b}_k(t)|k_\tau ,`$ (14)
where we have removed the center frequency $`\omega _0`$ from the wave packet amplitudes $`\stackrel{~}{b}_k`$. A classical discrete Fourier transform (DFT) relates the two sets of amplitudes,

$$\stackrel{~}{b}_k(t)=\frac{1}{\sqrt{d}}\underset{j}{}b_j(t)\mathrm{exp}(i\omega _{j0}t)\mathrm{exp}\left(i2\pi jk/d\right).$$
(15)
For free time evolution, $`b_j(t)=b_j(0)`$, and Eq. (15) can be written as

$$\stackrel{~}{b}_k(t)=\underset{k^{}}{}\stackrel{~}{b}_k^{}(0)u_{k^{}k}(t),$$
(16)
where $`u_{k^{}k}(t)`$ is the DFT of the time evolution kernel,

$$u_{k^{}k}(t)=\frac{1}{d}\underset{j}{}\mathrm{exp}\left[i\omega _{j0}ti2\pi j(k^{}k)/d\right].$$
(17)
In the Kepler regime of Eq. (4), the levels are equally spaced, $`\omega _{j0}2\pi j/T_K`$, and $`u_{k^{}k}`$ becomes a delta function in $`k^{}`$, resulting in a cyclic permutation of the wave packet amplitudes in Eq. (16),

$$\stackrel{~}{b}_k(nT_K/d)=\stackrel{~}{b}_{kn}(0),$$
(18)
where $`n`$ is an integer and $`kn`$ is taken modulo $`d`$. Equation (18) amounts to a SHIFT gate in the wave packet basis, the multivalued equivalent of the NOT gate for binary logic. More complex unitary transforms occur for longer times, when the revival and superrevival dynamics become relevant. Thus free time evolution alone becomes useful for quantum information processing in the wave packet basis. This is a consequence of the nonstationary character of the basis and the unequally spaced energy levels.

We now turn to the control of the multilevel atomic basis using laser fields. Universal multilevel transforms of the amplitudes $`b_j`$ can in principle be achieved using narrow-band lasers to couple neighboring transitions among the energy levels . However this requires multiple laser frequencies tuned to the $`d1`$ transitions. A simpler control scheme is possible for the wave packet amplitudes $`\stackrel{~}{b}_k`$, as shown below. This involves short broadband pulses and takes advantage of the cyclic evolution of the wave packet states described above.

A laser pulse with a spectrum overlapping the $`d`$ Rydberg levels couples most strongly to only that combination of the levels with a phase relation between them that localizes the electron near the core. This corresponds to the wave packet state $`|0_\tau `$, for which all the phases are equal since $`k=0`$. Consider a pulsed field with a carrier frequency $`\omega \omega _0`$ and pulse profile $`f(t)`$,

$$𝐄(t)=\mathit{ϵ}E_0f(t)\mathrm{exp}(i\omega t)+\text{c.c.},$$
(19)
where $`E_0`$ and $`\mathit{ϵ}`$ are the amplitude and polarization of the field. Assume that the pulse couples the $`d`$ Rydberg levels to the ground state $`|g`$ in the atom. In the energy level basis for the Rydberg states, the dipole interaction with the field yields the following coupled equations for time evolution in the rotating wave approximation,

$`\dot{b}_g`$ $`=`$ $`{\displaystyle \frac{i}{2}}f(t){\displaystyle \underset{j}{}}\mathrm{\Omega }_j\mathrm{exp}(i\mathrm{\Delta }_jt)b_j,`$ (20)
$`\dot{b}_j`$ $`=`$ $`{\displaystyle \frac{i}{2}}f(t)\mathrm{\Omega }_j\mathrm{exp}(i\mathrm{\Delta }_jt)b_g,`$ (21)
where $`\mathrm{\Omega }_j`$ and $`\mathrm{\Delta }_j`$ are the Rabi frequency and field detuning defined as

$`\mathrm{\Omega }_j`$ $`=`$ $`2g|\widehat{𝐝}|j_\nu \mathit{ϵ}E_0/\mathrm{},`$ (22)
$`\mathrm{\Delta }_j`$ $`=`$ $`\omega _j\omega ,`$ (23)
and we have taken $`\mathrm{\Omega }_j`$ to be real in Eqs. (20) and (21). Using the Fourier inverse of Eq. (15), we rewrite the equation of motion for $`b_g`$ in terms of the wave packet amplitudes,

$$\dot{b}_g=\frac{i}{2}f(t)\mathrm{exp}(i\mathrm{\Delta }_0t)\underset{k}{}\stackrel{~}{\mathrm{\Omega }}_k\stackrel{~}{b}_k,$$
(24)
where $`\mathrm{\Delta }_0`$ is the mean detuning of the field, and $`\stackrel{~}{\mathrm{\Omega }}_k`$ is the DFT of the Rabi frequency $`\mathrm{\Omega }_j`$,

$$\stackrel{~}{\mathrm{\Omega }}_k=\frac{1}{\sqrt{d}}\underset{j}{}\mathrm{\Omega }_j\mathrm{exp}\left(i2\pi jk/d\right).$$
(25)
To find the equation of motion for the wave packet amplitude $`\stackrel{~}{b}_k`$, we take the time derivative of Eq. (15) and

use Eq. (21) for $`\dot{b}_j`$. Assuming that the pulse width is short compared to a Kepler period, we can ignore the free atomic evolution during the pulse that led to Eq. (16). In the short pulse regime, we are left with

$$\dot{\stackrel{~}{b}}_k=\frac{i}{2}f(t)\stackrel{~}{\mathrm{\Omega }}_k^{}\mathrm{exp}(i\mathrm{\Delta }_0t)b_g.$$
(26)
Equations (24) and (26) describe Rabi oscillations between the ground state and the wave packet states, where the Rabi frequencies $`\stackrel{~}{\mathrm{\Omega }}_k`$ are Fourier transforms of their energy level counterparts $`\mathrm{\Omega }_j`$. For Rydberg levels, $`\mathrm{\Omega }_j`$ varies slowly with $`j`$,

$$\mathrm{\Omega }_j(\overline{n}+j)^{3/2},$$
(27)
making only the dc component $`\stackrel{~}{\mathrm{\Omega }}_0`$ significant for large $`\overline{n}`$. For $`\overline{n}=180`$ and $`d=8`$, we find $`\stackrel{~}{\mathrm{\Omega }}_{\pm 1}0.01\text{ }\stackrel{~}{\mathrm{\Omega }}_0`$. Thus to good approximation, we can ignore the Rabi frequencies $`\stackrel{~}{\mathrm{\Omega }}_{k0}`$ in Eqs. (24) and (26). These equations reduce to

$`\dot{b}_g`$ $``$ $`{\displaystyle \frac{i}{2}}f(t)\stackrel{~}{\mathrm{\Omega }}_0\mathrm{exp}(i\mathrm{\Delta }_0t)\stackrel{~}{b}_0,`$ (28)
$`\dot{\stackrel{~}{b}}_0`$ $``$ $`{\displaystyle \frac{i}{2}}f(t)\stackrel{~}{\mathrm{\Omega }}_0\mathrm{exp}(i\mathrm{\Delta }_0t)b_g.`$ (29)
These describe two-level Rabi oscillations between the ground state and the innermost wave packet state during the pulse, with a detuning $`\mathrm{\Delta }_0=\omega _0\omega `$ between the peak of the pulse spectrum and the average Rydberg frequency. The two-level approximation is valid only for short pulses. To affect only the $`k=0`$ wave packet state, we need the pulse width $`\tau _p`$ to be less than $`T_K/d`$, the time taken for this wave packet element to leave the atomic core. However, the bandwidth of the pulse cannot be much larger than the atomic frequencies spanning the $`d`$-level basis, to avoid coupling into levels outside this basis. Since the Rydberg frequency spacing is approximately $`1/T_K`$, this implies a bandwidth $`\mathrm{\Delta }\omega /2\pi d/T_K`$. To summarize, we require that

$$\tau _p<\frac{T_K}{d}\frac{2\pi }{\mathrm{\Delta }\omega },$$
(30)
or that the pulse be transform-limited. For $`\overline{n}=180`$ and $`d=8`$ levels, we use a Gaussian pulse of width $`\tau _p=0.5\mathrm{ln}2`$ $`T_K/d110`$ ps (FWHM amplitude), corresponding to a spectral width $`\mathrm{\Delta }\omega /2\pi =(8/2\pi )`$ $`d/T_K`$ for the electric field. Figure 2 shows an exact numerical solution for the wave packet amplitudes before, during and after the pulse. Assuming that the population is uniformly distributed initially, we find that a $`\pi `$ pulse transfers nearly all the population in the $`k=0`$ state to the ground state, creating a “dark” wave packet in the Rydberg basis that evolves cyclically in the Kepler evolution. About $`5\%`$ of the population is lost from neighboring states $`k=\pm 1`$ during the pulse. Thus a coherent coupling of individual wave packet states with the ground state can be modeled as a two-level system to good approximation.

Since the wave packet basis is nonstationary in time, the Rydberg wave function undergoes dispersion and revivals as the energy level phases evolve in the Schrödinger picture. This leads to aperiodicity in the free time evolution of the wave packet amplitudes $`\stackrel{~}{b}_k`$, going beyond the cyclic regime of Eq. (18). At $`\overline{n}=180`$ and $`d=8`$ for example, we find that the population in a wave packet state decays by about $`5\%`$ after a Kepler period, but recovers to within $`2`$-$`4\%`$ at a revival or super-revival. This scales as $`\overline{n}^2`$, making dispersion losses less significant for larger $`\overline{n}`$. By appropriate timing of successive laser pulses, we can thus take advantage of the revival structure of the atomic energy spectrum when addressing the wave packet amplitudes in the basis.

An arbitrary unitary transform of the $`d`$ wave packet states can be decomposed into a product of $`𝒪[d^2]`$ two-dimensional transforms , each of which can be implemented sequentially on two de-excited wave packet amplitudes at a time, as shown in Fig. 3. Two low lying energy levels $`g`$ and $`e`$ temporarily store the wave packet amplitudes in the form $`\stackrel{~}{b}_k|g+\stackrel{~}{b}_k^{}|e`$. Phase shifts or rotations of this state are accomplished by coupling the two levels with a pulse of appropriate length, Rabi frequency and detuning. The transformed state is then restored to the Rydberg basis. Thus multilevel transforms in the atom are possible with a timed sequence of laser pulses, without the need for multiple lasers to address the different Rydberg transitions. Moreover the fidelity of the gate operations compares favorably with that in an equivalent system of entangled two-level atoms .

For efficiently scalable quantum computing, we require conditional transforms on two entangled multilevel atoms. The method used in Ref.  is based on the linear ion trap scheme for binary quantum logic . This scheme provides a means to evolve one two-level atom in the trap array conditional on the state of another two-level atom. Since the wave packet approach allows each basis amplitude in the atom to be addressed sequentially through low-lying energy levels, conditional transforms on other atoms can be carried out using the same methods as in the binary scheme. That is, conditional on a given wave packet state (vis à vis energy level) in the control atom, a $`d`$-dimensional transform is performed on the wave packet basis of the target atom using the methods outlined above.

The time taken to couple two atoms in this approach can be on the order of the atomic Kepler period, which permits a fast non-adiabatic coupling scheme. Also, the mean energy and number of Rydberg levels used depend ultimately on the ionization threshold imposed by external fields in the environment. We anticipate that neutral atom schemes that allow fast coupling times, as have been suggested recently , will provide a useful setting for multilevel quantum information applications.

This work was supported by the Army Research Office through the MURI Center for Quantum Information.


format: text

text:
# Conductance fluctuations and boundary conditions

## Abstract

The conductance fluctuations for various types for two– and three–dimensional disordered systems with hard wall and periodic boundary conditions are studied, all the way from the ballistic (metallic) regime to the localized regime. It is shown that the universal conductance fluctuations (UCF) depend on the boundary conditions. The same holds for the metal to insulator transition. The conditions for observing the UCF are also given.

PACS numbers: 71.30.+h, 71.55.Jv

The influence of the boundary conditions (bc) on critical phenomena in disordered mesoscopic systems has been demonstrated by studies of the conductance distribution and energy level statistics. The ensemble average of the logarithm of the conductance, $`\mathrm{ln}(g)`$, is smaller for hard wall boundary conditions than for periodic boundary conditions. The variance $`\mathrm{ln}^2(g)\mathrm{ln}(g)^2`$ on the other hand is larger for hard wall than for periodic boundary conditions. The distribution of nearest neighbour energy level separations $`P(s)`$ becomes more Wigner–Dyson–like for periodic boundary conditions. Thus, systems with periodic boundary conditions exhibit a “more metallic” behaviour than those with hard wall boundary conditions. Different boundary conditions also lead to different values for the universal conductance fluctuations (UCF) in the diffusive metallic regime where the mean free path $`l`$ is much smaller and the localization length $`\xi `$ much bigger than the system size. The variation of ensemble fluctuations of the conductance as the disorder increases throughout the metallic regime has been studied, but the influence of the boundary conditions in the ballistic regime where $`l`$ exceeds the system size, and close to the localized regime where $`\xi `$ becomes comparable with the system size, has not been studied in detail. Also, the role of the correlation length in samples with a true metal–insulator–transition (MIT) — i.e. in systems where the localization length is “infinite” in the metallic regime — has not been discussed. We present here numerical studies of cubic systems and squares with spin–orbit scattering — all of which have a true MIT — for both hard wall and periodic boundary conditions in the direction(s) perpendicular to transport.

We are using the tight–binding model with the Hamiltonian

$$=\underset{n,\tau }{}|n\tau \epsilon _nn\tau \left|+\underset{n,\tau ,n^{},\tau ^{}}{}\right|n\tau V_{n,n^{}}n^{}\tau ^{}|$$
(1)
where $`n,n^{}`$ are nearest neighbour lattice sites on a square or cubic lattice. For systems with spin–orbit–interactions $`\tau `$ and $`\tau ^{}`$ take on values of $`+1`$ or $`1`$ and the hopping integrals $`V_{n,n^{}}`$ are $`2\times 2`$ matrices; without spin–orbit–interactions, the hopping integrals are scalar and the spin “variables” have only one value. We take the site energies $`\epsilon _n`$ (independent of $`\tau `$) to be random variables, chosen from an interval $`[W/2;W/2]`$ with a uniform probability distribution. The parameter $`W`$ serves thus as a measure of disorder strength. The conductance is calculated using the transfer matrix method and the Landauer formula.

We have used the analytical derivation of Lee et al. to calculate theoretical values for the UCF for both types of boundary conditions. In order to change the boundary conditions to periodic, one needs to make the following changes (references to equations are from the Appendix of Lee et al. (1987)): in the eigenfunctions $`Q_m`$ to the diffusion equations, the cosines in the transverse directions (Eq. (A9)) must be replaced by exponentials with a factor of $`2\pi `$ instead of $`\pi `$ in the argument; this will lead in effect to a factor $`4`$ in the $`m_x^2`$ and $`m_y^2`$ terms in the modified eigenvalues $`\stackrel{~}{\lambda }_m`$ (Eq. (A13)), and to a summation over all integers (including negative ones) for $`m_x`$ and $`m_y`$ in Eqs. (A15), (A16), (A24) and (A25). The results are presented in Table I. The values are only half those given by Lee et al. due to a factor $`2`$ in the definition of $`g`$. Also, our result for the three–dimensional case is slightly higher, probably due to our calculating the involved sums to a higher precision. The boundary conditions have of course no effect for the quasi–one–dimensional case. The values given here are those for the standard deviation $`\sigma _g`$ for the orthogonal universality class of Random Matrix Theory. The values for the other universality classes are obtained by dividing the variance, i.e. $`\sigma _g^2`$, by the universality class parameter $`\beta `$, where $`\beta =1,2,4`$ for the orthogonal, unitary, and symplectic universality classes respectively.

In Fig. 1 we show the standard deviation of the conductance in ensembles of 10,000 samples for different system sizes (squares with $`L\times L`$ lattice sites) and boundary conditions (open symbols: periodic boundary conditions; filled symbols: hard wall boundary conditions) as a function of the inverse of the average conductance. It is well known that in a two–dimensional disordered tight–binding model all the states are exponentially localized. The typical structure of the fluctuations with increasing disorder strength can be seen: after an initial strong increase in the ballistic regime (large $`g`$) it reaches a peak value which becomes more pronounced for larger systems; then the fluctuations drop back to the universal value and finally decrease again in the strongly localized regime. The boundary conditions have apparently no influence on the behaviour outside the region of UCF in this case. Notice that for the case of periodic boundary conditions for the large system size of $`L=128`$, $`\sigma _g`$ approaches the theoretical value of $`0.393`$ given by the lower horizontal line.

In Fig. 2 the same data is plotted for a two–dimensional disordered system of size $`L=64`$ with periodic boundary conditions. In the same plot the mean free path and the localisation length as a function of disorder strength $`W`$ are given. Both the mean free path and the localization length were obtained from the numerical results of Economou et al. The localization length was obtained by the transfer matrix method, while the mean free path was obtained by the coherent potential approximation (CPA). Notice that $`\xi `$ is always larger than $`l`$. So for a given system size ($`L=64`$ in this case), there is a finite region where $`lL\xi `$. Only in this region there is a plateau visible at the correct UCF–value. For $`W1`$, $`l`$ is larger than $`L`$, and we are in the ballistic regime where one observes a monotonic increase of $`\sigma _g`$ followed by the characteristic maximum as the system enters the crossover between the ballistic regime and the regime characterized by UCF.

Fig. 3 shows that the same overall behaviour is observed also for three–dimensional systems. As there is a MIT (indicated by vertical lines in Fig. 3) where the conductance distribution and therefore also its standard deviation become universal, i.e. independent of system size (though still depending on the boundary conditions), this value is approached after leaving the region of UCF. A direct comparison of the results shows again that the boundary conditions have only minimal effect outside that region. The additional peak noticeable in some of the periodic boundary conditions data are due to a near–degeneracy of eigen–energies for very small disorder.

In Fig. 4 we plot the data for one of the systems again as a function of $`W`$, together with the mean free path $`l`$ and the correlation length $`\xi `$. In the three–dimensional disordered case, the mean free path and the correlation length $`\xi `$ were again obtained by the CPA and the transfer matrix method. The size of the cube is $`L=16`$ and periodic boundary conditions were used. In the three–dimensional case, $`l`$ drops as the disorder strength $`W`$ increases, while $`\xi `$ increases as $`W`$ increases. The plateau in $`\sigma _g`$ is seen only when $`lL`$ and $`\xi L`$, so that we expect a wider plateau for larger systems. The fluctuations begin to approach the critical value as soon as $`L\xi `$. In the three–dimensional case too, when $`W2`$, $`l`$ is larger than $`L`$ and the ballistic regime is observed, followed by the maximum in the crossover regime.

Finally, Fig. 5 shows data for square systems with spin–orbit–interactions. We have chosen the Evangelou–Ziman model, where even in the absence of diagonal disorder there is disorder in the hopping matrices $`V_{n,n^{}}`$, which accounts for the fact that the fluctuations do not vanish for small diagonal disorder. Apparently, boundary conditions have a noticeable influence on the fluctuations even outside the region of UCF, but this is likely due to the peculiar overall structure in this case, most significantly the fact that the UCF value is much smaller than the critical value, causing another increase in the standard deviation as one approaches the MIT.

In conclusion, we have investigated the conductance fluctuations for various types of systems with both hard wall and periodic boundary conditions from the ballistic regime to the localized regime. The boundary conditions seem to have a relevant influence on the conductance fluctuation only in the region of UCF and at the critical point of the MIT. In true metallic systems, the fluctuations begin to deviate from the UCF–value and approach the critical value as soon as the correlation length approaches the system size.

Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. W–7405–Eng–82. This work was supported by the Director for Energy Research, Office of Basic Science. The authors thank Professor E. N. Economou for helpful discussion. P.M. would like to thank Ames Laboratory for their hospitality and support and the Slovak Grant Agency for financial support.


text:
# Untitled Document

Response to Technical Comments regarding the paper ”Direct Detection of Galactic Halo Dark Matter,” by B. R. Oppenheimer, N. C. Hambly, A. P. Digby, S. T. Hodgkin, D. Saumon, Science 292, 698 (2001).

Gibson and Flynn cite several arguments to support their claim that our study (1) greatly overestimated the space density of halo white dwarfs. The first stems from a recomputation of our $`1/V_{\mathrm{max}}`$ calculation. Their new value actually lies within the uncertainty in our estimate, however, so it is not clear to begin with that there is any discrepancy. Regardless of how the numbers were calculated, the uncertainties in the space densities are larger than the differences, and it is not meaningful to discuss them at the accuracies that Gibson and Flynn demand. Furthermore, Gibson and Flynn themselves do not uniformly apply their own criteria for what constitutes an ”important” level of difference. They clearly view as significant the change wrought by their initial recalculation, which implies a 14% decrease in the calculated white dwarf density, and the change wrought by applying a white dwarf mass estimate drawn from globular clusters, which implies an additional 17% decrease in the density (paragraph 3 of their comment). Yet they disregard as having ”little effect” a change of 11% in the opposite direction due to different assumptions for limiting apparent magnitude (paragraph 4 of their comment). Finally, the fact that we rounded our final number clearly indicates that we do not believe it is meaningful to distinguish between 1.6% and 2% based on current data. Deeper surveys are crucial for assessing the full extent of this population.

Gibson and Flynn also argue that our study used an average white dwarf mass, 0.6 M, that was too high by 20%. We do not regard this argument as valid. Although the average white dwarf mass in globular clusters is indeed nearer to 0.5 M than to 0.6 M, other research suggests that if our findings are correct and the population of stars represented by the white dwarfs in our study did emerge from a nonstandard initial mass function, the average white dwarf mass may actually be much higher. Chabrier (2), for example, found an average white dwarf mass of 0.7 to 0.8 M.

Gibson and Flynn claim that our results are inconsistent with the LHS catalog, pointing out that no halo white dwarfs are seen in the northern portion of the LHS. Because of the similarities between our survey and the LHS, they maintain, there can be no halo white dwarfs in our survey. That assertion is incorrect, in part because there are several crucial differences between our survey and LHS and in part because there clearly are halo white dwarfs in the LHS. We discuss these issues in greater detail below, in our response to the Graff comment.

Finally, Gibson and Flynn conclude that only 0.5% of the dark matter is explained by white dwarfs. The argument through which they reach that number, however, is not logical. They take the dimmest star, WD0351$``$564, and place the entire space density of all of the white dwarfs found in our survey into the bin of the luminosity function corresponding to that faintest star—ignoring the fact that it only accounts for 7% of the space density. Then they argue that our survey did not find enough stars in this faintest luminosity bin. This circular argument is used to claim that we have overestimated the space density. Further, even if one chooses to accept that flawed argument, one must still contend with the fact that the 0.5% figure is only a lower limit. The complex detection limits (in proper motion and magnitude) and the incomplete sampling of the velocity parameter space suggest that a substantially larger population exists and that we have only found the tip of the iceberg.

Graff’s arguments likewise do not withstand close scrutiny. He begins with the premise, also proposed by Gibson and Flynn, that our survey and the LHS catalog are very similar. As we have already noted, that premise is incorrect; our survey and the LHS are not directly comparable. Our survey reached more than a magnitude deeper and found many objects to which the LHS was not sensitive. A clear indication that we have not simply repeated the LHS work is contained in our reduced proper motion diagram \[figure 1 in (1)\]: The LHS catalog, as Gibson and Flynn mention, has relatively few objects with reduced proper motion $`H_\mathrm{R}>22`$; much of our sample, by contrast, is drawn from the stars with $`H_\mathrm{R}>22`$. The Liebert et al. (3) sample of LHS white dwarfs, on which Graff particularly focuses, likewise included stars with a smaller range of proper motions than we allowed in our survey. It is more logical to interpret these differences as indicating that our survey is deeper and more complete than to assume that the previous studies of white dwarfs in the LHS are complete and we have made a mistake.

In addition, one of the coolest white dwarfs known before our survey, LHS 3250, was first cataloged in the LHS and is notably absent from the previous studies of the white dwarf content of LHS. This star, and other cool white dwarfs that were previously cataloged, had never been studied spectroscopically before our study. Rather, the assumption was that their colors suggested that they were main sequence stars, not blue, cool white dwarfs with collision-induced absorption (4). Indeed, to this day, there remain many objects in the LHS that have not been measured spectroscopically.

These issues render comparison of our survey with the LHS, or the Liebert et al. sample of that survey, much more complicated than simple, back-of-the-envelope calculations based on the tables in Liebert et al., as proposed by Graff. In reality, a complete Monte Carlo simulation of the populations and the survey is necessary to assess the sensitivity of our survey in detail. Let us assume for the moment, however, that the Liebert et al. survey and our survey are indeed directly comparable, as Graff suggests. He claims that the stars we classify as halo stars would simply have been dubbed thick disk stars by Liebert et al. We point out, though, that the Liebert et al. sample contained stars that unambiguously are members of the galactic halo, along with others that may indeed be members of the halo. For example, the star LHS 542, which is in their survey, was clearly shown to be a halo member by Ibata et al. (5). The Liebert et al. (3) study is also missing stars: LHS 3250, which may or may not be a halo member, is conspicuously absent from the Liebert et al. study, partly because it was shown to be a white dwarf only recently (4), and Liebert et al. took known white dwarfs from LHS without a complete search through the LHS catalog for white dwarfs that had not been observed in detail. The halo star WD 0346+246 (6, 7) is missing from both Liebert et al. (3) and the LHS, even though it was within the photometric and proper motion detection limits of the LHS—which suggests that LHS is not complete at the R = 18.5 level for white dwarfs. The latter two stars both have peculiar spectral energy distributions.

In short, the LHS survey certainly does contain dark halo white dwarfs that have not previously been identified as such. Indeed, 11 of the 38 stars that we listed \[table 1 in (1)\] were in the LHS or LP catalog. However, the LHS catalog contains only a small number of halo white dwarfs, and it is certainly not complete at the detection levels necessary to reveal a convincing fraction of the halo white dwarf population.

Graff continues by claiming that because the local density of halo main sequence stars is 600 times smaller than the density of disk main sequence stars (8), the same should be true for the white dwarfs. That claim has no basis in our current understanding of these two different populations of stars. First of all, if the halo is composed of substantially older stars that formed roughly coevally, as is generally believed, one would expect substantially different ratios of white dwarfs and main sequence stars in the halo and disk. The disk is believed to be a population of stars that have been continuously forming since the disk formed. The comparison is thus moot, and the assertion that 1 in 30 of the stars in our sample may be halo white dwarfs is rendered incorrect. The first examples of halo white dwarfs were discovered convincingly only in the past few years, and the construction of relative numbers of these stars is impossible if one disregards the results that we published in (1).

According to Graff, because our survey is proper-motion limited, we necessarily have included more of the disk stars than we thought. Effectively, he claims, our 2-$`\sigma `$ exclusion is relegated to a 1-$`\sigma `$ exclusion of the disk stars. That statement clearly does not hold in all cases and, most important, in this case. Halo and even thick-disk stars should have average proper motions higher than those of the disk; their kinematics are necessarily different from those of the Sun. Furthermore, the 94 km s-1 number—which is actually centered at the point $`(V,U)=(35,0)`$, not $`(V,U)=(0,0)`$ as Graff seems to have assumed—comes from the survey by Chiba and Beers (9), which examined the velocity distributions of stars that were not kinematically selected. Therefore, there is no question that 94 km s-1 is a 2-$`\sigma `$ value.

To respond to Graff’s final point, we have not yet assessed the sensitivity of the survey in (1) as a function of proper motion with any accuracy. \[We did point out in (1) that there was a less than 10% chance that we would find any stars with 3 arcseconds of motion per year or greater.\] To assess that sensitivity—and, more important the sensitivity of our survey in the $`VU`$ parameter space that we plotted \[figure 3 in (1)\]—will require detailed modeling of the survey and the various galactic populations.

B. R. Oppenheimer

Astronomy Department, University of California, Berkeley,

Berkeley, CA 94720-3411, USA, E-mail: bro@astron.berkeley.edu

N. C. Hambly, A. P. Digby

Institute for Astronomy, University of Edinburgh, Royal Observatory,

Blackford Hill, Edinburgh, EH9 3HJ, UK

S. T. Hodgkin

Institute of Astronomy, Cambridge University, Madingley Road,

Cambridge, CB3 0HA, UK

D. Saumon

Department of Physics and Astronomy, Vanderbilt University,

Nashville, TN 37235, USA

References and Notes

1. B. R. Oppenheimer, N. C. Hambly, A. P. Digby, S. T. Hodgkin, D. Saumon, Science 292, 698 (2001).

2. G. Chabrier, Astrophys. J. Lett. 513, L103 (1999).

3. J. Liebert, C. Dahn, D. G. Monet, Astrophys. J. 332, 891 (1988).

4. H. C. Harris et al., Astrophys. J. 524, 1000 (1999).

5. R. A. Ibata et al., Astrophys. J. Lett. 532, L41 (2000).

6. N. C. Hambly et al., Astrophys. J. Lett. 489, L157 (1997).

7. S. T. Hodgkin et al., Nature 403, 57 (2000).

8. This point is actually a major source of contention among experts in this field and the number varies from Graff’s 600 to well below 250 \[(9); also (10) and references therein\].

9. M. Chiba and T. C. Beers, Astron. J. 119, 2843 (2000).

10. Gould, C. Flynn, J. N. Bahcall, Astrophys. J. 503, 798 (1998).

4 May 2001; accepted 18 May 2001


text:
# 1 Introduction

Magnetic Susceptibility Study of the Heavy Rare Earth Stannate Pyrochlores

V. Bondah-Jagalu1 and S. T. Bramwell2

1 Department of Chemistry, The University of Edinburgh, The King’s Buildings, West Mains Road, Edinburgh, EH9 3JJ, United Kingdom
2 Department of Chemistry, University College London, 20 Gordon Street, London, WC1H 0AJ, United Kingdom

Corresponding authors: e-mail: vin@ed.ac.uk and s.t.bramwell@ucl.ac.uk

The series of magnetic rare earth pyrochlore stannates R2Sn2O7 (R = rare earth, except Ce and Pm) have been investigated by powder susceptibility measurements down to T $`=1.8`$ K. The results are compared to results for the analogous titanate series, which are well-known frustrated magnets. Unlike the titanates, the whole series can be formed in the cubic pyrochlore structure. Possible experimental advantages of studying the stannates are discussed.

## 1 Introduction

The magnetic properties of the pyrochlore compounds can be very unusual on account of the high degree of frustration between antiferromagnetically coupled magnetic ions . Cubic pyrochlore oxides A2B2O6O’ are described by the face centred space group $`Fd3m`$ where the trivalent A atoms are situated on the 16d sites, the tetravalent B atoms are on the 16c sites and the oxygen atoms are in the 48f(O) and 8b(O’) sites. The A and the B atoms form identical interpenetrating sublattices that are spatially displaced from each other by $`\frac{1}{2},\frac{1}{2},\frac{1}{2}`$. Each lattice can be described as a three dimensional array of corner linked tetrahedra and it is this basic tetrahedral motif that gives rise to the geometrical frustration of antiferromagnetic interactions. The A site is coordinated by two oxygen atoms (O’) at the centre of the tetrahedra and by six oxygen atoms arranged in a puckered hexagonal ring around the metal ions. The structural details of pyrochlore compounds are well described in a very comprehensive review of the oxide pyrochlores . The tetrahedral lattice is a feature of many other materials, examples including

ternary fluorides (e.g. NH4CoAlF6 ) and spinels (e.g. MgCr2O4 ).
Investigation of the magnetic properties of the pyrochlores has been pioneered mainly by J .E. Greedan and collaborators. Of particular interest is the discovery  that several compounds of the form R2Mo2O7 (R = rare earth or Y) show frozen in magnetic disorder in the apparent absence of chemical disorder. For these compounds the molybdenum(IV) ions are magnetic (d2, S = $`1`$) and the R ion can be chosen to be nonmagnetic. Depending on the rare earth, the molybdates may be metallic (R = Sm, Nd, Gd) or semiconducting (R = Y, Tb-Yb)  and it is not entirely clear to what extend local moment models of frustration would apply to these materials. It therefore remains of interest to characterise new families of magnetic pyrochlore materials with localised moments. Here we describe a susceptibility study of a family of rare earth stannates R2Sn2O7 in which the Sn(IV) ion is nonmagnetic and the rare-earth ions are chosen to be magnetic. This study follows on from an investigation of the isostructural titanate materials R2Ti2O7  that in turn developed the early work of Cashion  and Blöte .

While investigating Ho2Ti2O7 it was realised  that in the presence of local $`\mathrm{1\hspace{0.17em}1\hspace{0.17em}1}`$ easy axis anisotropy, frustration arose from a net ferromagnetic coupling between spins which has subsequently been established to be largely dipolar in origin . The simple $`\mathrm{1\hspace{0.17em}1\hspace{0.17em}1}`$ Ising ferromagnet was termed the “spin ice” model on account of its precise mapping onto the problem of proton disorder in ice. It gives a good first approximation to the properties of Ho2Ti2O7  and Dy2Ti2O7 . Following on from the titanates, we decided to look at the magnetic properties of the rare earth stannates, whose structural properties have been well studied . The compounds are cubic, with a linear relationship between the radius of the lanthanide and the lattice. Non-cubic pyrochlore stannates do exist, such as Bi2Sn2O7 .

## 2 Experimental

Our samples were prepared from SnO2 and R2O3 (R = La, Nd, Sm, Eu, Gd, Dy, Ho, Er, Tm, Yb, Lu and Y) ordered from Strem and Aldrich, all of a minimum of 99.9% purity. For the praseodymium and terbium stannate, Pr6O11 and Tb4O7 were used respectively and it was assumed that the excess oxygen would escape as a gas during annealing. It was not possible to prepare cerium stannate in this way, since the cerium remains as Ce+4. However, because of this fact, it would be possible to make La2Ce2O7 , although this material adopts a disordered fluorite structure. The oxides were ground, pressed into pellets and heated in air to $`1400^o`$C for 4 hours, this length of time being found sufficient for all of the compounds. The X-ray diffraction refinement of room temperature data gave lattice constants consistent with the literature values .

Magnetisation was measured in the range $`1.8`$ K to $`300`$ K at an applied field of $`10`$ Oe ($`T<20`$ K) and $`100`$ Oe (T $`<20`$ K) using a Quantum Design SQUID magnetometer, and divided by field and amount of substance to give the susceptibility per mole of lanthanide. The low field susceptibility was scaled to fit the higher field data, as the field positioning at $`10`$ Oe is not very accurate. Zero- field cooled versus field cooled susceptibility was determined in an applied field of $`10`$ Oe between T = $`1.84`$ K and T = $`20`$ K, and magnetisation versus field isotherms were measured at T = $`1.8`$ K in applied fields of up to $`7`$ T. Estimated Curie-Weiss temperatures and magnetic moments are listed in table 1. Individual curves of effective moment $`\mu _{eff}=\sqrt{8\chi T}`$ versus temperature are shown in figure 1.

## 3 Results and Discussion

$`R_2Sn_2O_7`$ (R = Nd, Tb, Dy, Ho, Er, Yb) almost certainly have doublet ground states with strong single ion anisotropy. The latter is evidenced by the fact that the magnetisation approaches saturation at roughly half the free ion value. Amongst this group all except Yb2Sn2O7 have Curie constants that are consistent with a ground state for $`R^{3+}`$ consisting mainly of MJ = $`\pm `$ J for the free ion ground term. For R = Nd, Gd, Dy, Ho, Er, Tm, Yb the magnetic behaviour is extremely close to that of the analogous titanate. Thus one expects “dipolar spin ice” behaviour for R = Ho, Dy and ordering at T $`<1.8`$ K for R = Nd, Er, Yb, while R = Tm clearly has a singlet ground state derived from crystal field splitting of the free ion term . Both Er2Sn2O7 and Er2Ti2O7 exhibit a pronounced FC-ZFC splitting in the susceptibility at about $`3.4`$ K. This feature needs to be investigated further. For Tb2Sn2O7, the Curie-Weiss temperature $`\theta _{CW}`$ = $`11.1`$ K, is apparently much less than in Tb2Ti2O7, $`\theta _{CW}`$ = $`19`$ K . It is possible that our sample of Tb2Sn2O7 is slightly non-stoichiometric. Consistent with this, we observe a slight FC-ZFC splitting in the susceptibility below $`4.5`$ K, which is not reported for the titanate. The latter has been described as a co-operative paramagnet with no sign of magnetic order down to T = $`0.07`$ K .

There have not been detailed reports of the magnetism of titanates analogous to R2Sn2O7 with R = Pr, Sm, Eu, so we consider their behaviour in more detail.

$`Pr_2Sn_2O_7`$: The saturation magnetisation and effective moment approach the free ion value, while the Curie-Weiss constant is small. This suggests that crystal field effects and magnetic coupling are small in this material. This behaviour can be contrasted with Pr2O3  which has a large $`\theta _{CW}`$ $``$ -60 K arising from the crystal field splitting of the H6 free ion term of Pr3+ to give a ground state singlet and low lying doublet.

$`Sm_2Sn_2O_7`$: The effective moment versus temperature appears to have a term in (temperature)$`^{\frac{1}{2}}`$ dominating at high temperature, suggestive of a Van Vleck paramagnetism. It is typical for the susceptibility of Sm3+ salts to be dominated by the Van Vleck term near room temperature, and to show increasingly Curie-Weiss behaviour at low temperature . Ascribing the value $`0.00143`$ erg Oe-2 mol-1 to the Van Vleck term gives a Curie-Weiss temperature of $`1.36`$ K and moment derived from the Curie constant of $`0.30`$, which is less than the moment of the 6H$`_{\frac{5}{2}}`$ free ion term, $`0.84\mu _B`$ atom-1. This reduction in moment is consistent with the ground term in Sm2Sn2O7 being a doublet dominated by free ion MJ levels with MJ significantly less than the maximum, $`\frac{5}{2}`$. A similar effect is observed in Sm2O3, which shows a similar susceptibility curve .

$`\mathrm{𝐸𝑢}_\mathit{2}\mathrm{𝑆𝑛}_\mathit{2}O_\mathit{7}`$: The susceptibility is roughly constant below 100 K, behaviour typical for Eu3+, in which the free ion 7F0 ground term is non-magnetic, leaving only the Van Vleck term. Similar behaviour is observed in Eu2O3 . The observed magnitude of the Van Vleck term ($`0.008`$ erg Oe-2 mol-1) in Eu2Sn2O7 is close to that of Eu2O3, which has been shown to be enhanced by anisotropic exchange between the Eu ions .

## 4 Concluding Remarks

The rare earth stannates can be prepared in the cubic form across the whole lanthanide series and make an ideal system for the comparative study of magnetic properties. Experimental investigation of the series to lower temperatures is in progress. The presence of Sn in the structure will facilitate experimental investigation of these materials by Mossbauer 119Sn spectroscopy. For example, using this technique, Sn2B2O7 (B = Nb, Ta) have been shown to be actually of the form Sn$`{}_{}{}^{2+}{}_{22u}{}^{}`$(Sn$`{}_{}{}^{4+}{}_{2u}{}^{}`$B$`{}_{}{}^{5+}{}_{2u}{}^{}`$)O7-(5/2)u ; $`u`$ = $`0.180.20`$) . Solid solutions of the type RxR’1-xSn2O7 (R = rare earth, R’ = Ln, Y) can also be formed. In this context, 119Sn MAS NMR spectroscopy has been shown to be a particularly powerful probe of local cation geometry on account of the very large 119Sn chemical shifts . We finally note that, to complete the series, cerium stannate pyrochlore (Ce2Sn2O7) has recently been prepared for the first time .


text:
# Classical molecular dynamics simulations of amorphous silica surfaces

## I introduction

Most studies on surfaces have been focused on their electronic aspect. The goal was usually to analyze materials used in electronic devices which characteristics depend on the surface properties. As a consequence, crystalline surfaces are currently well known and a lot of attention has been devoted to the doping and transport properties of these surfaces as well as to the surface reconstruction. In the last years, many studies, theoretical  or experimental  have also been devoted to the analysis of liquid metal surfaces. In particular, the existence of an atomic layering in the density profile has been found. A great interest in understanding amorphous surfaces exists , because of their participation in the glass formation, and therefore several numerical studies are now available on this subject . Pioneering work in using classical potentials to describe amorphous silica surfaces is due to Garofalini et al. who studied mainly the presence of defects (non bridging oxygens (NBOs), small-sized rings..) at the surface . These defect sites are the most chemically reactive ones but experimental studies have, so far, not clearly shown their existence . An alternative way to learn more about the surface is to study the chemical reactions taking place on it  and in particular the interaction between water and amorphous SiO2 surfaces has been the topic of both experimental  and numerical studies . In most of the numerical studies available in the literature, the total number of atoms did not exceed 800 atoms and only a fraction of these atoms were considered as “mobile” (i.e. participating in the formation of the surface). (It should be noted that very recently (in fact during the process of writing the present article) a molecular-dynamics study of silica surfaces has been proposed involving several thousands of atoms ). Moreover the long-range Coulomb forces were either cut-off  or treated in a pseudo 2D geometry , not taking into account exactly the reduced dimensionality of free silica surfaces. A recent study of the influence of long-range forces on the structure and dynamics of a model silica glass has shown that over-damped Coulomb interactions can alter the quality of the results .

In the present work the molecular dynamics force calculation scheme has been adapted to the specific 2D geometry of samples with free surfaces. Our aim was mainly to study in detail the structural properties of amorphous silica surfaces with respect to what is already known for the corresponding bulk samples and to test the quality of the widely used “BKS” potential proposed by van Beest et al. . This classical pair potential describes quite well the structural , and vibrational , as well as relaxational  and thermal  properties of both, supercooled viscous liquid and glassy bulk silica samples: it is therefore justified to address its effectiveness in the case of “bi-dimensional” silica systems even more so since in a recent study based on ab initio simulations, Ceresoli et al. have put into question the structural description of the surfaces obtained with the BKS potential . In addition to the work already proposed by Roder et al.  who used the same potential but aimed their study mainly on silica clusters, we show that the effect of the surface is firstly reflected in the density profile which exhibits stratification, explained as a tendency to have more short range order at the surface. We then analyze in detail the ring size distribution as well as the orientation of the rings with respect to the surface, and show that the number of two-fold rings (whose signature is visible in the radial pair distribution functions and the angle distributions) is coherent with the experimental estimates and that they are positioned perpendicularly to the surface. By analyzing the influence of the annealing time on this effect, we show that it is a true surface reconstruction and not an artifact due to out-of equilibrium phenomena. Finally, in order to elucidate the diffusion process at the surface, we present the mean square displacement of the bridging and non-bridging oxygen atoms. We find no signature of a liquid like behavior at the surface but the dynamics of the NBOs seems to be anisotropic.

## II simulations

We have performed classical molecular dynamics simulations starting with cubic samples of edge length L=35.80 Å, containing 3000 particles. These samples were used previously to study the properties of bulk silica therefore their mass density is $`\rho 2.18`$ g/cm3, which is very close to the experimental value of 2.20 g/cm3 , and the periodic boundary conditions were assumed in the three space directions. The potential used to describe the interactions between the atoms is the two-body “BKS” potential. After the use of this potential to study densified silica samples , we apply it here to study silica surfaces assuming this potential is still valid in the case of SiO2 free surfaces. Of course this description can never be as accurate as an ab initio study, nevertheless we can simulate much larger samples and typical structural features should not depend on the details of the potential. The classical equations of motion are resolved by using the Velocity-Verlet algorithm with a timestep $`\mathrm{\Delta }t=`$ 0.7 fs. The initial configuration is given by a bulk sample of amorphous silica at $`T=0`$ K. This configuration has been obtained by quenching a well equilibrated liquid sample around 7000 K with a quench rate of $`2.3\times 10^{14}`$ K/s. The free surfaces are created by breaking the periodic boundary conditions along the $`z`$-direction, normal to the surface, thus creating two free surfaces located at $`L/2`$ and $`L/2`$. Doing so, the Ewald summation for the calculation of the long range forces has to be modified to take into account the loss of periodicity in the $`z`$-direction. There exists no standard method since several possibilities based on either a modified Ewald summation  or multipole expansions  are proposed in the literature. A purely two-dimensional method exists  which nevertheless necessitates the adjunction of charged plates. Here we have chosen a strictly two-dimensional approximate technique which has the advantage to be very simple . In the three dimensional Ewald summation technique, one considers all the point charges surrounded by a spherical distribution of charges that has same magnitude but opposite sign . Then, the potential is separated into two terms. The first term $`V_c`$ corresponds to the potential due to the point charge distribution minus that of the screening spherical charge distribution and is quite rapidly convergent in real space. The second one $`V_\gamma `$ is the potential created by the spherical distribution and is rapidly convergent in the reciprocal space. Since the real space term is a short range potential with a spherical symmetry, there is no need to change it for the two-dimensional summation. But, of course the reciprocal space potential has to be modified. While for the 3d Ewald summation the expression for the reciprocal space term of the potential is:

$`V_\gamma ={\displaystyle \frac{2\pi }{L^3}}{\displaystyle \underset{\stackrel{}{k}0}{}}{\displaystyle \frac{1}{k^2}}e^{\frac{k^2}{4\kappa ^2}}\left|{\displaystyle \underset{i}{}}q_ie^{i\stackrel{}{k}.\stackrel{}{r_i}}\right|^2`$ (1)
with

$`\stackrel{}{k}=n_x{\displaystyle \frac{2\pi }{L}}\stackrel{}{e_x}+n_y{\displaystyle \frac{2\pi }{L}}\stackrel{}{e_y}+n_z{\displaystyle \frac{2\pi }{L}}\stackrel{}{e_z},`$ (2)
In our approximate two-dimensional method it becomes:

$`V_\gamma ={\displaystyle \frac{\pi }{L^2}}{\displaystyle \underset{\stackrel{}{k}0}{}}{\displaystyle \frac{e^{\frac{k^2}{4\kappa ^2}}}{k}}\left|{\displaystyle \underset{i}{}}q_ie^{i\stackrel{}{k}.\stackrel{}{r_i}}\right|^2`$ (3)
with

$`\stackrel{}{k}=n_x{\displaystyle \frac{2\pi }{L}}\stackrel{}{e_x}+n_y{\displaystyle \frac{2\pi }{L}}\stackrel{}{e_y}.`$ (4)
Here $`\stackrel{}{k}`$ is a two-dimensional vector and the summation runs over the two integers $`n_x`$ and $`n_y`$. Although formula (3) looks like a straightforward two-dimensional transcription of formula (1) it is only approximate. The expression would have been exact if, instead of a Coulomb potential with a spherical symmetry one would have considered a cylindrical potential. In particular, one practical consequence is that the k-contribution to the $`z`$-component of the force is zero, since the scalar product $`\stackrel{}{k}.\stackrel{}{r_i}`$ is independent of the third coordinate $`z_i`$ of $`\stackrel{}{r_i}`$. We have used the characteristic constant $`\kappa =5.0/L`$, where $`L`$ is the cubic box size, and considered 49 k-vectors in reciprocal space to insure a relative error smaller than $`2\times 10^5`$ for the potential energy. The removal of the periodic boundary conditions in the $`z`$-direction breaks the bonds between the atoms at the “bottom” and those at the “top” of the film. This will lead to bond rearrangements near each free surface and of course increase the temperature of the system. It is also worth noticing that we let the position of the atoms adjust freely in the $`z`$-direction and therefore we are not stricto sensu in a microcanonical ensemble. To justify a posteriori the approximations and check the results we compared this method with the one in which the box length in the $`z`$-direction, Lz, is artificially increased in order to simulate a pseudo 2D system. We find that the two techniques give similar results when Lz is five times larger than the box length in the other 2 directions which is in contrast with the study of Bakaev et al.  in which Lz was taken only twice as large as in the $`x`$ and $`y`$ directions. In any case the method used here is computationally much more efficient since the number of k-points included in the Ewald summation has not to be increased. It is also faster than the more frequently used technique proposed in Ref.  which is not a strictly 2D technique and which is computationally very costly.
The heating of the slab mentioned earlier will tend to evaporate some oxygen atoms located at the surface, and to avoid this vaporization, it is necessary to control the temperature of the slab. This has been done at 1000K by rescaling the velocities during 30000 time steps. This temperature is below the simulated glass transition temperature as will be shown in a forthcoming study. Then the system was allowed to relax 30000 supplemental steps and finally we performed simulations during 60000 supplemental steps to collect the results, the system temperature remaining almost constant around 1000 K in this time range. Such a simulation represents 45 days of computer time on the latest IBM-SP2 processor and in order to improve the statistics of the data we averaged our results over ten independent liquid samples.

## III results

In figure 1 we have reported the averaged density profile of our samples in the $`z`$-direction. Each point (filled circle) of the figure corresponds to the density $`\rho (z)`$ calculated within a slab parallel to the surface of area $`L\times L`$ and of width $`\mathrm{\Delta }z=0.4475`$ Å and reported as a function of the mean slab position $`z`$. We have averaged the results over the two slabs located at $`z`$ and $`z`$. For comparison, we have reported the profile in the $`x`$-direction calculated within a slab of same width $`\mathrm{\Delta }x=0.4475`$ Å but of area $`L\times L_z`$ with $`L_z=26.85`$ Å, taken sufficiently smaller than $`L`$ to avoid edge effects due to the presence of the free surface in the $`z`$direction. In both cases, we get an average bulk density of about $`2.3`$ g/cm3, significantly larger than the density $`\rho =2.18`$ g/cm3 of the original cubic box with periodic boundary conditions. This increase of the density is due to the relaxation of the system (the system contracts along the $`z`$ direction) when creating the free surfaces since it has been shown that the equilibrium density of the BKS potential is larger than the experimental one  (this means that the initial sample with PBC had a negative pressure). This shows that creating free surfaces is a physical way to obtain the equilibrium density and thus permits also to validate a posteriori our modus operandi. For the profile in the $`z`$-direction, one can use a hyperbolic tangent fit derived from the Van der Waals theory of surface tension . Although this formula is generally used to study the thermodynamic properties of simple liquids, some characteristics of the density profile of our amorphous samples can be determined. This fit uses three parameters, $`\rho _0`$ the mass density in the underlying bulk, $`z_0`$ the position of the surface and $`d`$ the surface thickness:

$`\rho (z)={\displaystyle \frac{\rho _0}{2}}\left(1\mathrm{tanh}\left({\displaystyle \frac{2(zz_0)}{d}}\right)\right).`$ (5)
In fact, this formula fits very well the vapor-liquid interfaces which exhibit symmetric density profiles . But here it does not take into account the asymmetric density maximum near the free surfaces (see figure 1) and therefore the density profile is not very well fitted near this maximum. As a consequence, if the parameter $`\rho _0`$ ($`\rho _0=2.36g/cm^3`$) gives a good quantitative value for the bulk density, the parameter $`d`$ ($`d=1.18`$ Å ) gives only a rough (underestimated) value of the actual surface thickness. As observed for liquid surfaces , we have an asymmetric peak with a broad tail in the density profile at the surface and a layering behavior. Our goal is to identify the participation of the stratification in the oscillations observed in the density profile along the $`z`$-direction. However when comparing the density profiles in the $`z`$ and $`x`$ directions in figure 1, one can hardly see a difference in the oscillations around the plateau. One way to address this point is to compute the standard deviation of these data over different samples. If the oscillations are only due to statistical fluctuations, the standard deviation $`\sigma `$ of the density should behave like $`N_S^{1/2}`$, where $`N_S`$ is the number of statistical independent samples over which the density profiles are averaged. Along the $`x`$ and $`y`$ axis, we find that the standard deviations behave like $`N_S^{1/2}`$ while for the $`z`$-direction they are quite independent on $`N_S`$. Therefore the oscillations along the $`z`$-direction are most likely a signature of a layer formation whereas the oscillations observed in the $`x`$-direction are only due to statistical density fluctuations. This stratification is likely to be enhanced at the free surfaces and to extend in the underlying bulk. Computing the radial pair distribution for atoms within a sphere of radius $`R=L/4`$, whose center is the center of the the initial cubic box, we found that it is similar to that observed for bulk amorphous silica samples. This is an evidence that the layering behavior is not due to crystallization. The inset of Fig. 1, exhibits the contribution of each type of atom to the mass density. Since the stoichiometry is respected, the quantity $`\mathrm{\Delta }n`$ = 2\[Si\] - \[O\] is oscillating around a zero mean within the underlying bulk. Near the surface, the apparition of two peaks reflects some stoichiometry breaking. The silicon atom excess under the surface and the oxygen atom excess at the surface are respectively represented by the positive and the negative peak. This surface segregation appears not to be model dependent since it has been found also by Roder et al.  as well as by Litton et al.  with another potential.
In figure 2 we present the radial pair distribution functions (RPDFs) for the different pairs of atoms in order to investigate the structure of our samples. The space has been divided into two regions, the “surface” layer which extends 6 Å below the $`z`$ position of the outermost atom, and the “interior” which represents the rest of the sample. The width of the surface layer has been chosen to encompass the broad peak in the density profile visible in figure 1 and is coherent with previous studies . Also shown in figure 2 are the RPDFs obtained in a bulk silica sample. As a general trend, the RPDFs in the interior are very similar to the ones obtained in the bulk. Nevertheless differences can be seen especially in the O-O RPDF, where a slight shift towards smaller distances can be seen. This is a consequence of the relaxation towards a higher equilibrium density. The major difference between the surface and the interior can be seen in the Si-Si RPDF, since a small shoulder appears around 2.5 Å . This shoulder which can be attributed to 2-membered rings as we will see below, is absent in the interior or bulk curve since this structural unit is not present in the bulk . In figure 3 we show the distribution functions of the angle between an oxygen (silicon) atom and its 2 silicon (oxygen ) neighbors (two atoms are neighbors if their separation is smaller than the location of the first minimum in the corresponding RPDF) where we distinguish again the surface, the interior and the results obtained for bulk silica (this is a further test of the validity of the BKS potential which does not contain 3-body terms). Concerning the O-Si-O angle distribution, there is almost no difference between the interior and the bulk. The distribution at the surface is slightly broader and shows a second peak around 80 which is absent in the bulk. Concerning the Si-O-Si angle distribution we firstly note a slight shift of the whole distribution towards smaller angles compared to the bulk distribution which is again a consequence of the density relaxation. Next we note a strong shift of the surface distribution towards smaller angles by more than 10 which is coherent with the peak in the density profile observed in figure 1 and an important second peak around 100 which is absent in the bulk and interior curves. The existence of the additional peaks in the surface angle distributions is in agreement with previous studies . The location of these peaks is different from the angles O-Si-O and Si-O-Si expected in 2 membered rings : 94 and 86 respectively. Here we find 80 and 100 which indicates a weakness of the pair potential on this particular point since in a recent ab initio study on silica clusters Lopez et al. found 91 and 89 respectively .
So far we only have the signature of the presence of 2-membered rings at the surface therefore it is justified to address the influence of the free surfaces on the variation of the ring size distribution. Such quantity is of interest since it can be extracted from infrared and Raman spectroscopy. In particular the highly strained 2-membered rings result in infrared-active Si-O stretching modes at 888 and 908 cm-1 . The ring size distribution has been determined using the algorithm described in , where a ring is defined as the smallest loop starting from one oxygen atom O1 nearest neighbor of a given Si atom and ending on O2 another of its nearest neighbor. Then the size of a ring is given by the number of its constitutive Si-O segments. To get the distribution of each ring size along the $`z`$ axis in the 2 regions (interior and surface) we compute the probability $`P_n`$ for a given Si atom, whose coordinate along the normal direction to the surface is in a given region, to be a member of a $`n`$-fold ring. These results are reported in the top part of figure 4 and compared to the bulk distribution. We calculated the ring size distribution in the surface layer on the initial 2D samples (not relaxed) and on the 2D samples annealed for 42ps (relaxed) in order to see the influence of the annealing time on the results. Firstly it should be noted that the distributions in the interior region (relaxed or not relaxed) are very similar to the bulk distribution (which is a broad Gaussian centered around six-fold rings ) and are therefore not reported in figure 4a. The situation is different at the surfaces where the global trend after the annealing period is a decrease of the probability for a Si atom to be a member of a large ring (as expected, since a lot of bonds are broken) together with an increase of the probability to be a member of a small ring, more precisely, two-fold, three-fold and four-fold rings. A particular case of this behavior is the inversion of the proportion of five-fold and six-fold rings at the annealed surface. Energetically, there is no significant difference between these two ring sizes, therefore if six-fold rings are broken at the surface (because they are relatively large) a fraction of these rings will become five-fold and five-fold rings will become predominant: in a sense this corresponds to the reconstruction of an amorphous surface. An other significant feature is the relatively large proportion of two-fold rings at the annealed surface whereas it stays very weak in the bulk and at the freshly created surface. More precisely we find an average density of 0.25 rings per nm2 which is in the range of the experimental values which vary between 0.1 and 0.4 rings per nm2 . This proportion is smaller than the one found in ref  (0.6 rings per nm2), where the statistics were done at 3000K and for silica clusters at equilibrium (which is certainly not the case here).
Together with the ring size distribution we have also investigated the ring orientations by computing $`<\mathrm{cos}^2\theta >`$, for a given ring size $`n`$, within a given region of the sample, where $`\theta `$ is the angle between the normal to the surface and to the ring. The results for the surface region are reported in the bottom of figure 4 as a function of the annealing time, and compared with those obtained in a 3D sample (which are similar to those in the interior region). In the bulk, all the $`n`$-fold rings are oriented in an isotropic way since the value of $`<\mathrm{cos}^2\theta >`$ is close to $`\frac{1}{3}`$. In the freshly created surface region, the $`2,3`$ and $`4`$-fold have still an isotropic orientation, while the larger rings are oriented rather parallel to the surface since $`<\mathrm{cos}^2\theta >`$ is slightly larger than $`\frac{1}{3}`$. In the annealed surface layer we note that the rings larger than the three-fold ones have a strong tendency to orient parallel to the surface. On the contrary the two-fold rings are now perpendicular to the surface since $`<\mathrm{cos}^2\theta >`$ is much smaller than $`\frac{1}{3}`$. This is a result of the annealing and can therefore be extrapolated to a film in equilibrium in which two-fold rings are certainly oriented perpendicular to the surface. This is coherent with the ab initio study of Ceresoli et al.  even if they had only one two-fold ring at the surface.
After the description of the structure let us now turn in the final part of this work to the diffusion mechanism inside our thin silica films. The motivation of this study is connected to the idea that similarly to crystalline faces which develop a microscopic liquid film  the surface layers of an amorphous sample show an enhanced diffusivity compared to the bulk . The simplest way to address this question is to compute the mean-squared displacement (MSD) $`r^2(t)`$ = $`|𝐫_i(t)𝐫_i(0)|^2`$ for a given particle type $`i`$. This is what is represented in figure 5 for the oxygen atoms (the behavior of the silicon atoms is qualitatively similar as shown in Ref. ). The interior and surface regions have been determined with respect to the initial position of the atoms. In figure 5 we have also distinguished the behavior of the non-bridging oxygens (NBOs) from that of the bridging oxygens (BOs). At the surface we find approximately 15 % of NBOs while in the interior only 0.5 % of the oxygen atoms are single-coordinated, which is in good agreement with a recent MD study of nanoporous silica , especially for the number of NBOs in the surface region. We have also calculated the components of the MSD that are parallel ( $`1.5\times x^2+y^2`$) and perpendicular ($`3\times z^2`$) to the surface, and these components are also represented in figure 5. In order to study the influence of the annealing time, the MSD has been calculated for 42ps at the beginning of the relaxation process (thin lines) and after an annealing time of 105ps (bold lines). For the BO there is no difference in the MSD between the interior and the surface region. After the annealing period, the displacement of the atoms is smaller and more isotropic. Hence for a fully equilibrated film the dynamics of the particles inside the sample is expected to be almost isotropic.
For the NBOs the situation is different. Firstly as expected, the NBOs diffuse more than the BOs. Next, after annealing, the dynamics appear to be more anisotropic in both regions (interior and surface). In the interior region the perpendicular MSD is more important than the parallel one as if the NBOs were attracted towards the surface. On the contrary in the surface layer the NBOs move parallel to the surface rather than perpendicular to it. This can be understood since most of the Si-O bonds point outward of the surface and since it is easier for an NBO to have a bond “bending” motion rather than a bond “stretching” motion. Therefore concerning the NBOs, we find that their dynamics is anisotropic both in the interior and the surface region. Since this anisotropy is enhanced after annealing one can expect that this behavior is also true in an equilibrated sample.

## IV conclusion

We have used classical molecular dynamics to study the structural properties and the mean squared displacement of the particles inside an amorphous silica film. The first step in this work was to adapt our computational method to samples without periodic conditions in one spatial direction. This has been done using an original method which is fast and relatively accurate. With this technique we have calculated the density profile along the normal direction $`z`$ to the surface. We have found that the system is contracting since the equilibrium density of the potential is greater than the mass density of the initial amorphous silica sample. This relaxation is rather fast and creating a surface seems therefore to be a physical and efficient way to obtain the equilibrium density of a given system. This shows also that the BKS potential describes properly the physics inherent to the creation of free surfaces. Like in liquid metals, the density profile exhibits a layering behavior that seems to be enhanced at the surface. This feature was interpreted as a sign of an increase of the short range order near the surface. The asymmetric peak near the surface is a sign of a stoichiometry breaking: while the silicon atoms are in excess in the inner part of the surface the oxygen atoms are in excess in its outer part as already mentioned in previous studies . We have calculated also the radial pair distribution functions and the angle distributions and we found the signature of the presence of 2-fold rings at the surface similarly to previous studies . Therefore we have also focused our attention on the ring size distribution and ring orientation along the $`z`$-direction. We see that at the free surfaces, the most energetically unfavorable configurations (small rings) are enhanced and correlatively, the large rings show a tendency to disappear. A particular case of this feature is the inversion of the proportion of five-fold and six-fold rings. Also the orientation of the small-sized rings with respect to the surface (2-fold rings are perpendicular to the surface) agrees well with a recent ab initio study . Concerning the MSD of the oxygen atoms, we find that the dynamics of the BOs is unaffected by the presence of the surface. On the contrary the dynamics of the NBOs is clearly anisotropic: in the interior (far from the surface) they move rather perpendicular to the surface while in the surface region they move parallel to it.


text:
# Ultrahigh reflection from a medium with ultraslow group velocity

\[

## Abstract

We show that an incident wavepacket at the boundary to a medium with extremely slow group velocity, experiences enhanced reflection and a substantial spatial and temporal distortion of the transmitted wave packet. In the limit of vanishing group velocity, light cannot be transferred into the medium due to its perfect reflectivity.

\]

Currently there is a great deal of interest in slowing down light in dispersive media. It is triggered by a continuous experimental progress using electromagnetically induced transparency (EIT) in gaseous media. The attained group velocities are as small as 17 m/s for the medium being a Bose–Einstein condensate, 90 m/s in a hot rubidium gas, and even 8 m/s in a rubidium vapor at room temperature . EIT provides here a large enhancement of the dispersion of medium while absorption is largely suppressed .

A theoretical description of optical radiation propagating with group velocities $`v_g10^8c`$ (where $`c`$ is the speed of light in vacuum) or even slower, may require a careful reconsideration of some approximations routinely employed in optics. For example, the reflection of an optical beam at the boundary of the medium is usually considered to be determined by the refractive index of the medium. This is perfectly valid for a monochromatic beam, and serves as a perfect approximation for a wavepacket in a typical dielectric with weak dispersion, where each spectral component experiences essentially the same refractive index. However, when the dispersion steepens, a substantial variation of the refractive index over the spectrum of a wavepacket may be present. Then the usual approximation fails and new effects emerge that are dominated by the spectrally dependent refraction of the transmitted light. Besides a severe spatio-temporal distortion of the transmitted beam, one might expect difficulties in actually transferring the incident optical wavepacket into a medium of such strong dispersion.

In this Letter we show the significance of these effects for optical experiments with extremely slow group velocity. We point out a fundamental obstacle preventing one from infinitely decreasing the speed of light inside a medium: The slowing down of light inside the medium is accompanied by an increase of the reflectivity of the boundary of the medium. Near the limit of stopping light, almost no energy can penetrate from outside (as well as leak from inside) and the boundary of the medium will behave as a perfect mirror.

Let us consider the reflection and transmission of an optical wavepacket at a plane boundary between vacuum and a linear dielectric with frequency dependent dielectric function $`\epsilon (\omega )`$. Since a relatively small absorption will have little effect on our conclusions we take $`\epsilon (\omega )`$ to be real valued . Expanding the dielectric function of the medium around the center frequency $`\omega _0`$ of the incident wavepacket and keeping only the leading term, we arrive at

$$\epsilon (\omega )=n_0^2+\frac{2\alpha }{\omega _0}(\omega \omega _0)+\mathrm{},$$
(1)
where $`n_0`$ is the (absolute) refractive index of the medium at $`\omega =\omega _0`$ and $`\alpha `$ determines the steepness of the dispersion . In the medium, the magnitude of the wave vector and the dielectric function are related by $`k(\omega )=\omega /c\sqrt{\epsilon (\omega )}`$, from which we obtain for the group velocity inside the medium,

$$v_g=\frac{c}{n_0+\alpha /n_0}.$$
(2)
The transmittance at the boundary of the medium can be illustrated in more detail by considering the intensity reflection and transmission coefficients for the case of normal incidence :

$$R(\omega )=\left|\frac{1\sqrt{\epsilon (\omega )}}{1+\sqrt{\epsilon (\omega )}}\right|^2,T(\omega )=\left|\frac{2}{1+\sqrt{\epsilon (\omega )}}\right|^2,$$
(3)
As can be seen from Fig. 1, slowing down the group velocity brings forth two distinct features in the reflectivity: (a) For frequencies below the increasing cutoff frequency,

$$\omega _c=\omega _0\left(1\frac{n_0^2}{2\alpha }\right),$$
(4)
the low-frequency wing experiences total reflection, and (b) an overall increase of the reflectivity for the high-frequency wing. Usually these effects play no role in optical experiments. However, for a small ratio $`v_g/c`$, the cutoff frequency may move into the spectrum of the incident optical wavepacket. One clearly observes from Fig. 1, that whereas for $`c/v_g=5\times 10^7`$ and a wavepacket of 1$`\mu \mathrm{s}`$ duration the reflectivity is still small, for lower group velocities the cutoff frequency approaches $`\omega _0`$, while the reflectivity for the high-frequency tail also continuously increases. In the extreme limit $`v_g/c0`$ we then obtain overall reflection of the incident optical wavepacket for all frequencies apart from the center frequency.

The parameter $`r`$ defined as the ratio of the characteristic spectral width of the incident wavepacket $`\mathrm{\Delta }\omega `$ and the passband,

$$r=\frac{\mathrm{\Delta }\omega /2}{\omega _0\omega _c}=\frac{\alpha }{n_0^2}\frac{\mathrm{\Delta }\omega }{\omega _0}\frac{c/n_0}{v_g}\frac{\mathrm{\Delta }\omega }{\omega _0},$$
(5)
can be used as a measure for the relevance of the boundary effects. When $`r>1`$ the cutoff frequency considerably moves into the spectrum of the incident wavepacket, so that the variation of the refractive index over the spectrum leads already to an enhanced reflection and substantial spectral (temporal) distortion of the wave packet. The spectral modifications are asymmetric with respect to the center frequency $`\omega _0`$. Since the low-frequency tail experiences a larger reflection compared to the higher frequencies, this results in an overall blue shift of the transmitted wavepacket.

For the ideal case of EIT with purely radiative transitions of rate $`\gamma `$, negligible Doppler effect, and long-living coherence, one finds $`n_01`$ and $`r3/(8\pi ^2)N\lambda ^3\gamma \mathrm{\Delta }\omega /\mathrm{\Omega }^2`$ with $`N`$ and $`\mathrm{\Omega }`$ being the atomic density and the Rabi frequency of the driving field, respectively. The relevant spectral width of the wavepacket may be taken as the EIT transparency window $`\mathrm{\Delta }\omega \mathrm{\Omega }^2/\gamma `$, and we immediately obtain the estimation $`r3N\lambda ^3/(8\pi ^2)`$. In conclusion, for optical fields, say of wavelength $`\lambda =800\mathrm{nm}`$, $`r`$ may already approach unity for rather modest densities of $`N5\times 10^{13}\mathrm{cm}^3`$.

In general, the importance of the boundary effects can be demonstrated by the energy transfer into the medium, characterized by the ratio: $`T_E=E_t/E_i`$. Here $`E_t`$ and $`E_i`$ are the total energies of the transmitted and incident wavepackets, respectively, defined as the normal component of the Poynting vector in and outside the medium, integrated over time and the transverse area of the beam. Fig. 2 shows a typical example of how the energy-transfer ratio $`T_E`$ changes with group velocity for an incident Gaussian wavepacket. The continuous decrease of $`T_E`$ with decreasing group velocity indicates a vanishing energy flow through the boundary for $`v_g/c0`$.

The above result remains conceptually valid for a dielectric function of a more general form than Eq. (1). Given that any additional term appearing in Eq. (1) is small compared to the first-order term, we expect only a minor modification of the cutoff frequency and of the overall index of refraction.

Analogous effects arise when the sharp boundary is replaced by a layer of finite width $`L`$, where the dielectric function varies continuously along the direction of propagation. Such a regime is readily available in EIT with the continuous variation of $`\epsilon `$ being provided by the externally controllable driving field. Despite arbitrary intermediate modifications, the reflection properties depend only on the difference between the initial and final values of the dielectric function, $`\epsilon (\omega ,z=0)`$ and $`\epsilon (\omega ,z=L)`$, respectively . Then our results are recovered upon replacement of $`\epsilon (\omega )`$ by $`\epsilon (\omega ,z=L)`$ in Eq. (1).

For the more general situation where the incident beam has some spread, we expect along with temporal also spatial distortions. Typically the spread of the transmitted beam will increase with increasing refractive index. Approaching the extreme case of very large dispersion, when the cutoff frequency appears inside the spectrum, the lower frequency components are totally reflected, whereas the other components are refracted over a wide range of angles.

The enhanced reflectivity from the highly dispersive medium draws attention also to the importance of backscattering effects in the field dynamics in inhomogeneous/nonlinear media. Whenever the parameter $`r`$ (or $`1/T_E`$) gains appreciable values, the widely used approximation of unidirectional propagation may fail. Instead, the full wave equation must be used for a proper inclusion of the counter-propagating waves.

The occurrence of a cutoff frequency where wave propagation is inhibited is a well known phenomenon in the ionosphere and waveguides . However, compared to these cases we deal here with a dielectric function of a conceptually different form. To illustrate this, we compare in Fig. 3 the frequency dependence of the group velocity for these cases. Whereas, in the case of the ionosphere or waveguides the group velocity is approaching $`c`$ for increasing frequency, in our case a maximum can be observed, after which the group velocity is continously falling to zero. That is, in our case a suppression of wave propagation is present also for higher frequency.

In conclusion, an incident wavepacket at the boundary to a highly dispersive medium with extremely slow group velocity experiences enhanced reflection and a substantial spatio-temporal distortion of the transmitted wavepacket. In the case of EIT, the frequency selective reflectivity may be utilized for generating pulses (in reflection) revealing spectral shapes sensitively controllable by the external driving field. One may also expect that spectrally resolved polarization measurements of the reflected beam can provide an efficient tool for probing the dispersive properties of such media. On the other hand, the transmission losses due to the enhanced reflectivity pose a serious restriction on possible applications of slow light for quantum-optical purposes. Finally, we have shown, that in the extreme limit of vanishing group velocity, light cannot penetrate the boundary due to the perfect reflectivity of the medium.

S.W. acknowledges support from the Studienstiftung des deutschen Volkes by a BASF research fellowship.