We present the moving picture of a molecular bond, in phase-space, in real-time, at resolution limited by quantum uncertainty. The images are tomographically reconstructed Wigner distribution functions (WDF) obtained from four-wave mixing measurements on Br2-doped ice. The WDF completely characterizes the dissipative quantum evolution of the system, which despite coupling to the environment retains quantum coherence, as evidenced by its persistent negative Wigner hole. The spectral decomposition of the WDF allows a direct visualization of wavefunctions and spatiotemporal coherences of the system and the system-bath interaction. The measurements vividly illustrate nonclassical wave mechanics in a many-body system, in ordinary condensed matter.

Dynamical systems are most completely described by their evolution in phase space, which can be represented and visualized by Wigner's quasi-probability distribution function (WDF) along conjugate momentum p and position q coordinates, W(p, q; t).1–3 The quasi-probability designation of the WDF underscores that unlike classical distributions it can locally assume negative values, and the negative volume of the WDF serves as a measure of non-classicality.4 Given a continuous record of the marginal (shadow5) of the WDF, such as the observable probability density along coordinate space ρ(q; t) = ∫W(q, p; t)dp, the WDF can be reconstructed through the inverse Radon transform.6 This tomographic procedure is well established in quantum optics,7,8 and has been implemented to reconstruct the quantum state of a photon,9–12 of ions in an ion trap,13 of atoms in an atom beam,14 and of diatomic sodium molecules in the gas phase.15 Strictly positive WDFs arise in classically coherent collective dynamics.16 Superpositions that interfere in phase space exhibit negative Wigner holes are manifestly non-classical; and as such, are broadly designated as Schrodinger “cat” states.4 W-holes have been demonstrated in Bose condensates of photons11,12 and atoms.14,17 Here, we present the tomographic reconstruction of the dissipative motion of a molecule in ordinary matter: bromine in ice. The essential vibronic dynamics of this system has been previously characterized through electronically resonant transient grating (TG) measurements.18 The spectrally resolved extension of the same method (SRTG), which we implement here, allows imaging of the marginal ρ(q, t) = ψ*(q, t)ψ(q, t) of the molecular bond after setting it into motion on its electronically excited state. Although measurement effects cannot be entirely eliminated, a faithful reconstruction of the WDF is possible and its negative hole can be demonstrated. We capture and present the complete wave mechanical motion of a bond in real time, with spatial resolution limited by quantum uncertainty, in moving images that vividly illustrate wave-particle duality in molecular matter.

The observations are remarkable because they are made in ordinary condensed matter, in which quantum states are expected to rapidly collapse through decoherence,19 effected by energy transfer (dissipation) and loss of deterministic phase evolution (dephasing) as a system interacts with its environment. The observed quantum coherent motion of the Br2 bond develops after electronic dephasing and after extensive vibrational dissipation. Despite coupling to the surrounding, the environmentally induced quantum coherence persists at the measurement temperature of T = 120 K, as evidenced by the persistence of its Wigner hole. Since the WDF provides a complete description of the state of the system, it can be spectrally decomposed to directly observe wavefunctions and the spatiotemporal coherence of the system-bath interaction. We visualize the entanglement |s, b⟩ + |s′, b′⟩ that appears as the inseparable superposition of the Br2 bond and the breathing motion of the cage of water molecules within which it enclathrates.

The experiments are conducted on a thin solid film of Br2/H2O obtained by condensing a premixed vapor (1/500 mole ratio) on a sapphire plate held at 120 K. Under these conditions, the solid is expected to be low-density amorphous ice.20 Indeed, the visible absorption spectrum of Br2/ice is indistinguishable from that of Br2 dissolved in water, as shown in Figure 1. In both cases, the transitions that give bromine its brown color are shifted and broadened. As is typical to ordinary condensed matter, the combination of structural disorder and fast electronic dephasing render the spectra featureless. There is no evidence of a discrete level structure, which is a prerequisite for quantum evolution.

FIG. 1.

Absorption spectra of Br2 in the gas phase (black), in water or ice (white), and of Br2 clathrate hydrate crystal (red). The structure of the Br2 hydrate lattice viewed along the z-axis is overlaid on the image of the single crystal. In the clathrate, the lone pairs on oxygen are hydrogen bonded, and are not available to interact with the guest molecule.

FIG. 1.

Absorption spectra of Br2 in the gas phase (black), in water or ice (white), and of Br2 clathrate hydrate crystal (red). The structure of the Br2 hydrate lattice viewed along the z-axis is overlaid on the image of the single crystal. In the clathrate, the lone pairs on oxygen are hydrogen bonded, and are not available to interact with the guest molecule.

Close modal

Bromine is known to form clathrate hydrates — a ubiquitous class of crystalline solids in which the guest molecule occupies regular polyhedra made of hydrogen-bonded water molecules.21 In crystalline hydrates, the Br2 absorption bands undergo moderate shifts22 and long Raman progressions are observed, which are characteristic to the various cages the molecule is known to occupy.23 In water, the Raman spectrum is quenched.24 These differences are explained by local structure.25 In water and disordered ice, H2O may bind to Br2 via its lone pair of electrons, as in the H2O-Br2 complex.26 This geometry is not accessible in the clathrate cage where O atoms are fully H-bonded (Fig. 1). That a sub-ensemble of molecules enclathrates in ice was inferred through nonlinear time-resolved spectroscopy.18 Independent of initial excitation energy, at wavelengths ranging from 495 nm–565 nm, a coherent vibrational superposition consisting of v = 0–2 was observed to form on the electronic B-state, after the loss of 3000–4000 cm−1 in energy (Fig. 2(b)).18 Measurements in single clathrate crystals show that the wavepacket is prepared during rebound, after the stretching molecule collides with its immediate surrounding, and that the curve crossing located near v = 3 limits the highest vibrational level that can be trapped.27 It is this indirect, environmentally induced vibrational coherence that we characterize here.

As a spectrally resolved, electronically resonant, coherent four-wave mixing process, SRTG probes the phase-space.28 The scheme is illustrated in Fig. 2. Two identical non-collinear laser pulses with wavevectors k1 and k2 are used in coincidence to prepare a superposition of two vibrational wavepackets |B, k1⟩ + |B′, k2⟩ on the electronic B-state, and a delayed pulse resonant with the T ← B charge-transfer transition is used to induce coherent scattering between them, P(3) = ⟨B′, k2|μ(k4)|T⟩⟨T|μE(k3)|B, k1⟩. This third order polarization over a volume of ∼103 μm3 forms the signal which propagates along k4 = k1− k2+ k3. The measured material response is understood with the aid of the time-circuit diagram (Fig. 2(c)). During t1, the excitation evolves as a vibrational density |$\rho _{B,B^\prime } (t_1)$|ρB,B(t1). The probe pulse creates the electronic coherence, ρB,T(t2; t1), which evolves during t2 to close the measurement loop through radiation P(3)(t) = Tr[μρB,T(t2; t1)]. The two-time (t2; t1)-correlation that determines P(3) maps on the time-frequency domain of conjugate variables (ω, t1):

\begin{equation}S(\omega,t_1) = \left| {\int\limits_{ - \infty }^\infty {dt_2 e^{i\omega t_2 } } P^{(3)} \left( {t_2 ;t_1 } \right)} \right|^2.\end{equation}
S(ω,t1)=dt2eiωt2P(3)t2;t12.
(1)

The measurement consists of recording the power spectrum of this time correlation along t2 at a given delay t1 between writing and reading of the grating.

FIG. 2.

(a) Experimental arrangement: pulses along k1 and k2 prepare the grating; the delayed pulse along k3 serves as probe, the polarization along k4 is recorded in TG, while the dispersed spectrum is recorded in SRTG. (b) Potential energy diagram: the preparation is resonant with the molecular B ← X transition; the probe is resonant with the T ↔ B charge-transfer transition. The vibrational coherence appears after dissipation, below the curve crossing near v = 3. (c) Time circuit diagram of the molecular response: during t1 the vibrational coherence |B⟩ + |B′⟩ evolves, the probe prepares the electronic superposition |B⟩ +|T⟩, which evolves during t2. The radiation closes the measurement loop.

FIG. 2.

(a) Experimental arrangement: pulses along k1 and k2 prepare the grating; the delayed pulse along k3 serves as probe, the polarization along k4 is recorded in TG, while the dispersed spectrum is recorded in SRTG. (b) Potential energy diagram: the preparation is resonant with the molecular B ← X transition; the probe is resonant with the T ↔ B charge-transfer transition. The vibrational coherence appears after dissipation, below the curve crossing near v = 3. (c) Time circuit diagram of the molecular response: during t1 the vibrational coherence |B⟩ + |B′⟩ evolves, the probe prepares the electronic superposition |B⟩ +|T⟩, which evolves during t2. The radiation closes the measurement loop.

Close modal

The SRTG data are shown in Fig. 3(a). The oscillatory signal is unmistakably the wavepacket of a damped oscillator, albeit mapped along the spectral coordinate. Slices of the signal at early time show spectral sub-structure due to recursion in the signal. This can be seen in the (t2, t1)-correlation map of Fig. 3(b), which is obtained by Fourier transforming the signal along the spectral coordinate:

\begin{equation}S(t_2,t_1) = \int {S(\omega,t_1)\,e^{i\omega t} } d\omega.\end{equation}
S(t2,t1)=S(ω,t1)eiωtdω.
(2)

The slices along the t2-axis provide a detailed dissection of the electronic dephasing process of the probe coherence. They establish that the |B⟩⟨T| coherence is limited to inertial motion (Gaussian decay), with dephasing time equal to the probe pulse width of 60 fs. The recursion at t2 = 180 fs characterizes the curvature of the T-state by its harmonic frequency: ωe(T) = 185 cm−1. Frequency analysis along t1 characterizes the B-state potential by its harmonic frequency ωe(B) = 173.2 cm−1 and anharmonicity ωexe(B) = 1.75 cm−1. These values are modified in ice in comparison to those of the free molecule: ωe = 167.6 cm−1 and ωexe = 1.08 cm−1.29 The similarity of frequencies in B and T states indicates that the molecule retains its bond-order during excitation, consistent with the assignment of the probe resonance to Br2:H2O(s)+(T)  ←  Br2(B) charge-transfer transition. The similarity of curvatures in Br2(B) and Br2 potential minima30 implies that the difference potential ΔV(q) = VTVB is a linear function of the bond length, q, and the Franck-Condon factors are nearly diagonal |$|\langle v({\rm B})|v^\prime (T)\rangle |^2 \sim \delta_{v,v^\prime }$||v(B)|v(T)|2δv,v. Given the linear relation ω = ΔV(q) and the signal is trivially transformed from (ω, t1)-plane to (q, t1)-plane as implemented in Fig. 3(c). As long as ∂ΔV(q)/∂q ≠ 0, the color of the radiation defines the instantaneous bond length. Both turning points of the packet are seen, therefore the full amplitude of motion Δq remains within the probe bandwidth Δω = 350 cm−1 given by the convolution of pump and probe laser spectra. The requirement for capturing the complete marginal of the evolving density, (∂ΔV(q)/∂qq ⩽ Δω, is fulfilled. Finally, the slope of the difference can be directly read-off the terminal width of the signal, which is given by the thermal width of the nearly stationary state.

FIG. 3.

(a) The SRTG signal with overlaid trajectory of a damped, anharmonic oscillator. The spectral modulation, most prominent in the horizontal section at 0.41 ps, is due to a weak recurrence in the polarization. (b) Fourier transform of the signal along ω, in log10 scale, generates the (t2, t1)-correlation map. The recursion (∼1%) at t2 = 180 fs identifies the vibrational period of the T-state. The time slices track the dephasing of the probe coherence. (c) The back-transform of S(t2,t1) truncated at t2 = 100 fs forces the instantaneous scattering limit and eliminates the spectral modulation. The abscissa is transformed to the molecular coordinate using the linear difference potential between B- and T-states: ΔV(q) = 24 700 + 4520(q − 2.65) in cm−1, Å units.

FIG. 3.

(a) The SRTG signal with overlaid trajectory of a damped, anharmonic oscillator. The spectral modulation, most prominent in the horizontal section at 0.41 ps, is due to a weak recurrence in the polarization. (b) Fourier transform of the signal along ω, in log10 scale, generates the (t2, t1)-correlation map. The recursion (∼1%) at t2 = 180 fs identifies the vibrational period of the T-state. The time slices track the dephasing of the probe coherence. (c) The back-transform of S(t2,t1) truncated at t2 = 100 fs forces the instantaneous scattering limit and eliminates the spectral modulation. The abscissa is transformed to the molecular coordinate using the linear difference potential between B- and T-states: ΔV(q) = 24 700 + 4520(q − 2.65) in cm−1, Å units.

Close modal

To eliminate the interference due to evolution past the probe pulse, we truncate t2 at 100 fs and back-transform the data:

\begin{equation}S(\omega (q),t_1) = \int_0^{t = 100fs} {S(t_2,t_1)\,e^{i\omega t_2 } dt_2 }.\end{equation}
S(ω(q),t1)=0t=100fsS(t2,t1)eiωt2dt2.
(3)

The filtered signal, plotted against the coordinate axis in Fig. 3(c), shows the motion of the bond with amplitude that decays as it relaxes to form the stationary Gaussian v = 0 state. Inspection shows that the oscillating distribution is skewed, and this can be directly traced to the finite time of measurement (see below). The signal in this imposed instantaneous scattering limit is treated as the experimental marginal after correcting for bandwidth, by dividing it with the spectrum at t = 0:31 

\begin{equation}\rho ^e (q,t) = S(q,t_1 - t_0)/S(q,t_0).\end{equation}
ρe(q,t)=S(q,t1t0)/S(q,t0).
(4)

The WDF is constructed through the inverse Radon transform of ρe(q, t), using the filtered back-transform algorithm over continuous 100-fs time segments that span a half-period of motion.32 

Snapshots of the continuously normalized WDF are shown in Fig. 4, and the complete animation is provided in Video S1. As highlighted by the zero-contours, we clearly see a persistent negative hole. A variety of metrics have been defined to quantify non-classicality, of which a convenient measure is the negative volume fraction:4 

\begin{equation}\delta {\rm \_(t)} = \frac{{\int\!\!{\int {dpdq[ {| {W(p,q;t)}| - W(p,q;t)} ]} } }}{{\int\!\!{\int {dpdqW(p,q;t)} } }},\end{equation}
δ_(t)=dpdq[|W(p,q;t)|W(p,q;t)]dpdqW(p,q;t),
(5)

which is plotted in Fig. 4. The measured W-hole is a periodic function of time, with amplitude that decays from an initial value of 0.3 to 0. Contrary to isolated systems, the interrogated density does not undergo unitary time evolution. It evolves and decays due to contact with the environment. The decay reflects relaxation of the system toward v = 0 — a state that does not sustain a hole: |$\delta \_(v = 0) = 0$|δ_(v=0)=0. The periodic modulation is an artifact of the measurement. We illustrate the measurement effect in Fig. 5 by considering pure states of the isolated molecule for which exact calculations are possible. The origin of the modulation is clearest when considering the measurement in W-representation, where the signal is formed by the overlap of the WDF of the B-state and its copy on the T-state. Evolution under the pulse corresponds to rotation of the respective WDFs with relative phase velocity determined by detuning Δ = ℏω − Te. As such, the measured distribution is subject to angular convolution under the probe pulse of width τ,

\begin{equation}W^e (\bar \vartheta (t)) = \int_\vartheta ^{\vartheta + \Delta \tau } \langle W(\vartheta)| W(\vartheta - \vartheta ^{\prime}) \rangle \,d\vartheta ^{\prime}.\end{equation}
We(ϑ¯(t))=ϑϑ+ΔτW(ϑ)|W(ϑϑ)dϑ.
(6)

Thus, the finite time of measurement leads to modulation, smoothing, and addition of wings to the WDF (see Fig. 5). The extent of distortion depends on details of potentials, laser bandwidth and detuning. However, independent of specifics, the visibility of the hole can only be reduced upon measurement since the full convolution of two W-functions leads to a strictly positive distribution.2 In effect, the measured δ is a lower bound of the true negative volume fraction. We may conclude that the evolving density is non-classical. It retains this character throughout the 4 ps of evolution. Moreover, a weak coherence transfer from the environment can be noticed in the slow growth of the hole near 2.5 ps.

FIG. 4.

Snapshots of the reconstructed WDF at indicated times (in fs): the W-hole is highlighted by the white zero-contour line. Reference cat states: (a) Gaussian running cat (Eq. (7) of text) matches the initial state; (b) three-state superposition of quantum harmonic oscillator at t = 50 fs.; (c) superposition of v = 0 and v = 1 matches frames near 1 ps; (d) mostly v = 0, which describes frames during the last ps; (e) graph of the negative volume fraction. For complete animation, see Video S1 (enhanced online). (Video S1) [URL: http://dx.doi.org/10.1063/1.4813437.1]

FIG. 4.

Snapshots of the reconstructed WDF at indicated times (in fs): the W-hole is highlighted by the white zero-contour line. Reference cat states: (a) Gaussian running cat (Eq. (7) of text) matches the initial state; (b) three-state superposition of quantum harmonic oscillator at t = 50 fs.; (c) superposition of v = 0 and v = 1 matches frames near 1 ps; (d) mostly v = 0, which describes frames during the last ps; (e) graph of the negative volume fraction. For complete animation, see Video S1 (enhanced online). (Video S1) [URL: http://dx.doi.org/10.1063/1.4813437.1]

Close modal
FIG. 5.

The measurement effect: (top row) WDF of pure states |$\sum\nolimits_v {c_v |v\rangle }$|vcv|v; (second row) explicit simulation of their measurement. The state amplitudes are: (a) c0 = 0, c1 = 0.8, c2 = 1, (b) c0 = 0.1, c1 = 1, c2 = 0.7, (c) c0 = 1, c1 = 0.5, c2 = 0. The graph is the time-dependent negative volume fractions of the measurement-simulated WDFs (d)–(f), with actual values of the three states indicated in the graph as blue line (d), green line (e), and red line (f), respectively.

FIG. 5.

The measurement effect: (top row) WDF of pure states |$\sum\nolimits_v {c_v |v\rangle }$|vcv|v; (second row) explicit simulation of their measurement. The state amplitudes are: (a) c0 = 0, c1 = 0.8, c2 = 1, (b) c0 = 0.1, c1 = 1, c2 = 0.7, (c) c0 = 1, c1 = 0.5, c2 = 0. The graph is the time-dependent negative volume fractions of the measurement-simulated WDFs (d)–(f), with actual values of the three states indicated in the graph as blue line (d), green line (e), and red line (f), respectively.

Close modal

The nonclassical nature of the motion of the bond is clearly seen. Consider for example the bimodal density along the outer ring in the first frame of Fig. 4. It represents a bond that is simultaneously stretched and compressed, with negative probability amplitude for being at the classical equilibrium geometry. The distribution imitates the Schrodinger cat state of Fig. 4(a), which is constructed from a superposition of two displaced Gaussians:4 

\begin{equation}W_{cat} (q,p) = W_ + + W\_ + \gamma W_{{\rm int}}.\end{equation}
Wcat(q,p)=W++W_+γW int .
(7)

The comparison is to a “running cat,” with non-zero momentum: q0 = 0.8, p0 = 1.8 and reduced visibility γ = 0.5.33 This is the initial superposition state in which the bond is prepared with two laser pulses |B(k1t)⟩ + |B′(k2t)⟩, after rebounding from the cage of water molecules.27 The bimodal preparation of the density may be understood to arise from a trajectory ensemble that splits along two paths, e.g., one-bounce and two-bounce paths, to fall below the curve crossing shown in Fig. 2. In the process, the surrounding is driven into motion, and system and bath are entangled |s, b⟩ + |s′, b′⟩.34 To the extent that the process is sudden, classical coherence subject to phase noise would be expected. This is familiar in both solid35 and liquid36 phase pump-probe measurements where wavepacket motion is observed, which for the most part can be reproduced using classical trajectory ensembles.37 A bimodal classical distribution centered at q+ and q would be given by the first two terms in (7). The quantum interference between the two paths Wint, which gives the negative hole, is the nontrivial observation. As in any two-slit experiment, the interference identifies wave mechanics. Here, the interference of the entangled superposition implies phase coherent wave mechanics in a many-body environment. Since the measurement is limited to the Br-Br coordinate, the observable is the wave-mechanical motion of the Br2 bond, which can be recognized by comparison to the reference states of the quantum harmonic oscillator provided in Fig. 4(b)–4(d). The early evolution is in a superposition of v = 0, 1, 2, and relative amplitudes can be selected to emulate individual frames. As v = 2 decays, the quantum equivalent of an oscillating bond appears as the c0|0⟩ + c1|1⟩ superposition shown in Fig. 4(c). As v = 1 decays, so does the hole (Fig. 4(d)). The spiral distributions seen in Fig. 4 (e.g., at 825 fs and 1375 fs) illustrate the distortion due to angular convolution (6), as in the simulations in Fig. 5. Note that the reference distributions are for pure states. In effect, after a period of evolution during which the system loses ∼0.5 eV of energy, in what must be regarded as classical evolution since it occurs in the continuum of states seen in the absorption spectrum of Fig. 1, a quantum coherent superposition of the system is prepared. Remarkably, while the prepared state dissipates, its quantum interference persists even though at 120 K both system and bath degrees of freedom are thermally accessible. The spectral decomposition, which we consider next, allows a dissection of the underlying dynamical contributions.

Motion in a quantum state is contained in the off-diagonal elements of the density matrix. Given the small number of contributing vibrational states, it is possible to decompose the spatiotemporal coherence into its principle spectral components, as a sum of system and system-bath contributions:

\begin{eqnarray}\rho (t) &=& \rho _s (t) + \rho _{sb} (t)\nonumber\\&=& \sum\nolimits_{v,v' = 0}^2 {c_{v,v'} (t) | v \rangle \langle {v'} |e^{ - i\omega _{v,v'} t} + \rho _{sb} (t)} .\end{eqnarray}
ρ(t)=ρs(t)+ρsb(t)=v,v=02cv,v(t)|vv|eiωv,vt+ρsb(t).
(8)

The task is simplified by the nearly harmonic level structure and the prevailing weak damping limit. Accordingly, the spectral components can be retrieved by adding and subtracting a time-shifted copy of ρe(q, t) from itself, with time-shift taken at half-period of the characteristic frequencies:38 |$\omega _{v,v^\prime } = \omega _s$|ωv,v=ωs, 2ωs, ωb, i.e., at fundamental, overtone, and bath frequencies. The residue yields the non-oscillatory diagonal density |$\sum\nolimits_v {p_v (t)\left| v \right\rangle \left\langle v \right|}$|vpv(t)vv, with slowly varying time-dependent population weights, pv(t). We present the stepwise dissection of the spatiotemporal coherence and the WDF of components.

The spatiotemporal coherence at the fundamental frequency is isolated by using the fundamental period, Ts = 194 fs, in the shift-subtract operation:

\begin{equation}\Delta ( - T_S /2)\rho ^e \equiv \rho ^e (q,t_1) - \rho ^e (q,t_1 + T_S /2) = 2\rho _{v,v + 1}^e.\end{equation}
Δ(TS/2)ρeρe(q,t1)ρe(q,t1+TS/2)=2ρv,v+1e.
(9)

The extracted map is shown in Figure 6(a), with the abscissa in cm−1 units of the raw data. We will refer to this axis as the spatial coordinate, since its linear ω → q transformation (which was implemented in Fig. 3(c)) is apparent by inspection. Its spatial sections (Fig. 6(b)) can be recognized as somewhat distorted products of familiar harmonic oscillator wavefunctions φ0(q1(q) and φ1(q2(q) of similar appearance. At early time the spatial profile along a horizontal cut is distorted due to evolution under the pulse, and this effect can be removed by considering the profile on a moving frame, i.e., along a skewed line (see Fig. 6(b)). The required slope, ∼100 fs, identifies the slippage to be due to motion of the overtone during the measurement. Consistent with this, the effect disappears as v = 2 decays (Fig. 6(b)). The temporal profile (Fig. 6(c)) is modulation at the fundamental beat, as clarified by the power spectrum (Fig. 6(d)). The spectrum establishes the purity of the isolated component. The fundamental beat decays with a time constant of t = 1.17 ps.

FIG. 6.

Fundamental coherence: (a) spatiotemporal profile, (b) spatial cuts taken along the indicated lines: the horizontal cut (red) at t = 110 fs shows the distorted profile of a|0⟩⟨1| + b|1⟩⟨2|, which is restored along the skewed line cut (blue). The skewing caused by the measurement corresponds to ∼90 fs delay across the width of the coherence. The skewing disappears at later time (as v = 2 decays), as seen in the horizontal cut (green) at 1425 fs. (c) The fundamental modulation decays with a time constant of 1165 fs. (d) Power spectrum of C identifies the purity of the decomposed coherence at the fundamental frequency ωs = 173 cm−1.

FIG. 6.

Fundamental coherence: (a) spatiotemporal profile, (b) spatial cuts taken along the indicated lines: the horizontal cut (red) at t = 110 fs shows the distorted profile of a|0⟩⟨1| + b|1⟩⟨2|, which is restored along the skewed line cut (blue). The skewing caused by the measurement corresponds to ∼90 fs delay across the width of the coherence. The skewing disappears at later time (as v = 2 decays), as seen in the horizontal cut (green) at 1425 fs. (c) The fundamental modulation decays with a time constant of 1165 fs. (d) Power spectrum of C identifies the purity of the decomposed coherence at the fundamental frequency ωs = 173 cm−1.

Close modal

The shift-sum operation approximates the diagonal density, reinforces the overtone coherence, and adds a phase shift to components that evolve more slowly:

\begin{eqnarray}\Delta ( + T_S /2)\rho ^e & \equiv & \rho ^e (q,t_1) + \rho ^e (q,t_1 + T_S /2)\nonumber\\&=& \rho _{v,v}^e + 2\rho _{v,v + 2}^e + A\rho _{SB}.\end{eqnarray}
Δ(+TS/2)ρeρe(q,t1)+ρe(q,t1+TS/2)=ρv,ve+2ρv,v+2e+AρSB.
(10)

The spatiotemporal coherence at the overtone frequency can be isolated from (9) using Ts/4 time for the shift-subtract operation:

\begin{equation}\Delta ( - T_S /4)\Delta ( + T_S /2)\rho ^e = 4\rho _{v,v + 2}^e.\end{equation}
Δ(TS/4)Δ(+TS/2)ρe=4ρv,v+2e.
(11)

The result is shown in Figure 7. The modulation in time occurs at the center of the spatial distribution, as would be expected for the φ0(q2(q) cross-term, identifiable by the spatial cuts in Fig. 7(b). The modulation decays with a time constant of 0.57 ps, at nearly half the fundamental lifetime of 1.17 ps. The power spectrum of the time slice (Fig. 7(d)) shows the overtone beat, 2ωs = 330 cm−1, flanked by weak sidebands that identify system-bath combination beats. Note that a 2:1 ratio of decay rates between v = 2 and v = 1 is to be expected for a harmonic oscillator subject to linear coupling to the environment, since γv, v − 1∝|⟨v′|q|v⟩|2 = |⟨v′|a + a†|v⟩|2 = vγ1, 0.

FIG. 7.

(a) Overtone coherence. (b) Spatial profiles along the horizontal show the same distortion as the fundamental: the horizontal cut at t = 110 fs is distorted, while along the skewed line we recover the spatial profile of the |0⟩⟨2| coherence. (c) The modulated time profile decays with a time constant of 567 fs. (d) The power spectrum of C peaks at the overtone frequency of ∼ 330 cm−1, flanked by weak sidebands due to the bath, at −ωb = 40 cm−1 and +ωb = 80 cm−1.

FIG. 7.

(a) Overtone coherence. (b) Spatial profiles along the horizontal show the same distortion as the fundamental: the horizontal cut at t = 110 fs is distorted, while along the skewed line we recover the spatial profile of the |0⟩⟨2| coherence. (c) The modulated time profile decays with a time constant of 567 fs. (d) The power spectrum of C peaks at the overtone frequency of ∼ 330 cm−1, flanked by weak sidebands due to the bath, at −ωb = 40 cm−1 and +ωb = 80 cm−1.

Close modal

The diagonal density is obtained by the repeated shift sum operation:

\begin{equation}\Delta ( + T_S /4)\Delta ( + T_S /2)\rho ^e = \rho _{i,i}^e + A\rho _{SB}.\end{equation}
Δ(+TS/4)Δ(+TS/2)ρe=ρi,ie+AρSB.
(12)

The resulting map is shown in Figure 8. It is void of any modulation at Br2 vibrational frequencies; however, a faint, slow modulation can be seen in its time profile (Fig. 9(c)). The population decays with a time constant of t = 1.9 ps. Its bimodal spatial profile at early time (Fig. 8(b)) identifies it as being dominated by |$\varphi _1^2 (q) + \varphi _2^2 (q)$|φ12(q)+φ22(q). At t = 1.4 ps, the system can be seen to have decayed into |$\varphi _1^2 (q) + \varphi _0^2 (q)$|φ12(q)+φ02(q). We see the vibrational dissipation through the wavefunctions.

FIG. 8.

(a) Map of the diagonal density. (b) The bimodal spatial profile at t = 110 fs (red) can be recognized as the sum of three states |$\sum\nolimits_{v = 0}^2 {a_v |v\rangle \langle v|}$|v=02av|vv| of a harmonic oscillator, with dominant contributions from v = 1 and 2. In the cut at t = 1425 fs (green), v = 2 is a minor component. (c) The weak modulation on the time profile is due to the bath. The population decays with a time constant of 1.9 ps. (d) The power spectrum of C establishes that the system coherences have been effectively filtered out.

FIG. 8.

(a) Map of the diagonal density. (b) The bimodal spatial profile at t = 110 fs (red) can be recognized as the sum of three states |$\sum\nolimits_{v = 0}^2 {a_v |v\rangle \langle v|}$|v=02av|vv| of a harmonic oscillator, with dominant contributions from v = 1 and 2. In the cut at t = 1425 fs (green), v = 2 is a minor component. (c) The weak modulation on the time profile is due to the bath. The population decays with a time constant of 1.9 ps. (d) The power spectrum of C establishes that the system coherences have been effectively filtered out.

Close modal

The system-bath coherence is isolated by repeated shift-sum operations, by varying Tb and demanding maximum visibility:

\begin{equation}\Delta ( + T_b /2)\Delta ( + T_S /4)\Delta ( + T_S /2)\rho ^e = A\rho _{SB}^e.\end{equation}
Δ(+Tb/2)Δ(+TS/4)Δ(+TS/2)ρe=AρSBe.
(13)

The result is shown in Figure 9. The time-course of this coherence shows that after an initial swing, it rings for the duration of the population, with a characteristic frequency of ωb = 38 cm−1. A second frequency near the overtone of the bath, at 2ωb ∼ 80 cm−1, is apparent in the time profile (Fig. 9(b)), and can be seen as the combination band in both system-bath and overtone coherence (Figs. 9(d) and 7(d)). Such persistent lattice motion is usually reserved to zone-boundary phonons.39 Although the system is embedded in ice, contrary to phonons of amorphous ice,40 the cage motion appears as a nondispersive resonant mode. The observed frequencies are characteristic to breathing modes of the clathrate cage observed in neutron scattering.41–43 The observed dynamical localization of the cage modes distinguishes clathrate hydrates from ice and the large disparity in their thermal conductivities. Although an extended clathrate is not expected in this 1:500 guest:host composition, the persistence of the system-bath coherence suggests that we are interrogating a sub-ensemble of enclathrated molecules. We recognize that the local cage is effective in mechanically isolating the molecule from the rest of the environment. Indirectly, the molecule also appears to be electrically screened as if in a Faraday cage. The latter is inferred by noting that the measurements would not be possible if the probe transition were inhomogeneously broadened. Given the large dipole of the charge-transfer transition, its decoupling from dipolar disorder of amorphous ice requires effective screening; presumably, through the rapid mechanical reaction of the hydrogen-bonded cage and clathrate network to dielectric fluctuations in the medium. In effect, the primary bath coordinate that is directly coupled to the system can be identified as the cage-breathing mode, which in turn is poorly coupled to the rest of the solid.

FIG. 9.

Entangled system-bath coherence. (a) Spatiotemporal coherence. (b) Spatial profiles along horizontal cuts show that at t = 326 fs (red) and 775 fs (blue) the coherence is nearly entirely negative. Much later, at t = 2354 fs, it oscillates about zero. (c) Time profiles along the indicated lines in (a) show the modulation of the system by the bath, and the negativity of this coherence at early time. (d) Power spectrum of the black curve in C yields the principle oscillation frequency, ωb = 38 cm−1, which is assigned to the breathing mode of the local clathrate cage.

FIG. 9.

Entangled system-bath coherence. (a) Spatiotemporal coherence. (b) Spatial profiles along horizontal cuts show that at t = 326 fs (red) and 775 fs (blue) the coherence is nearly entirely negative. Much later, at t = 2354 fs, it oscillates about zero. (c) Time profiles along the indicated lines in (a) show the modulation of the system by the bath, and the negativity of this coherence at early time. (d) Power spectrum of the black curve in C yields the principle oscillation frequency, ωb = 38 cm−1, which is assigned to the breathing mode of the local clathrate cage.

Close modal

The WDFs of the spectral components are illustrated by phase-space snapshots in Fig. 10 along with reference states of the harmonic oscillator. The WDF of the diagonal density is simply the weighted sum of wavefunctions in W-representation (see Figs. 10(a)–10(c) and Video S2)—once the WDF is reconstructed tomographically, or interferometrically,44 wavefunctions become observables.45 The comparison with the reference states highlights the fidelity with which the system wavefunctions are observed. The coherences can be seen as the time dependent quantum interference between members of the superposition. The coherence at the fundamental frequency (2nd column in Fig. 10) starts as highly non-classical c0,1|0⟩⟨1| + c1,2|1⟩⟨2|, and only after dissipation of v = 2 does it approach the odd distribution |0⟩⟨1| (see Figs. 10(d) and 10(e), and Video S3), which conforms to the classical notion of a bond that alternatively stretches and compresses. The |2⟩⟨0| overtone coherence (3rd column in Fig. 10), which lasts for the duration of the v = 2 lifetime, does not have a classical equivalent — the bond simultaneously materializes at its stretched and compressed extremes (Fig. 6(f) and Video S4). Since this spectral component does not have any interferences, its evolution is a simple rotation of the |2⟩⟨0| phase space distribution (see Video S4). While we clearly see dissipation by the change in composition of the states, there is no evidence of dephasing or decoherence along the system coordinate. Clearly, the system is not isolated. However, the system coherence is isolated through the projective measurement |$\hat P_S$|P̂S along the Br2 coordinate, enforced by the resonant interrogation of the density, which has the effect of tracing out the bath, Tr[|$\hat P_S$|P̂Sρ] = Trb[ρ]. To the extent that the total density |si, bk⟩⟨sj, bl| can be decomposed as product states of system and bath, the coherence of the system can be isolated Trb[ρ] = ρij δk,l in diagonal bath states and therefore time evolution strictly given by the system. Identifying the system states as the vibrations of Br2, we obtain the first term in Eq. (8).

FIG. 10.

Spectral decomposition of the WDF (upper matrix) and harmonic oscillator states (lower matrix). The stationary (ω = 0) component is the population-weighted sum of wavefunctions in W-representation: |$\sum\nolimits_v {p_v (t)W(\psi _v)}$|vpv(t)W(ψv). The three frames are compared to states with vibrational populations in v = 0, 1, 2 of (a) 0.075, 0.55, 0.37, (b) 0.37, 0.37, 0.26, and (c) 0.625, 0.375, 0. The fundamental coherence (ωe) starts as the sum of two cross-terms 0.55|0⟩⟨1|+0.44|1⟩⟨2| and evolves into |0⟩⟨1|, (d) and (e). The overtone (2ωe) coherence |0⟩⟨2| exists at early time. Its comparison with (f) clearly shows the measurement distortion due to the angular convolution. The system-bath coherence (ωb) does not have a reference state. Its time evolution at late time shows periodic momentum kicks the system receives from the local cage, which appear as the pure momentum state seen in the last frame. The line graphs are profiles along the indicated cuts on the experimental frames. The reference line profiles are along p = 0, except for (f), which is taken along a curved path to account for the angular convolution. Complete animation provided in Videos S2–S5 (enhanced online). (Video S2) (Video S3) (Video S4) (Video S5) [URL: http://dx.doi.org/10.1063/1.4813437.2] [URL: http://dx.doi.org/10.1063/1.4813437.3] [URL: http://dx.doi.org/10.1063/1.4813437.4] [URL: http://dx.doi.org/10.1063/1.4813437.5]

FIG. 10.

Spectral decomposition of the WDF (upper matrix) and harmonic oscillator states (lower matrix). The stationary (ω = 0) component is the population-weighted sum of wavefunctions in W-representation: |$\sum\nolimits_v {p_v (t)W(\psi _v)}$|vpv(t)W(ψv). The three frames are compared to states with vibrational populations in v = 0, 1, 2 of (a) 0.075, 0.55, 0.37, (b) 0.37, 0.37, 0.26, and (c) 0.625, 0.375, 0. The fundamental coherence (ωe) starts as the sum of two cross-terms 0.55|0⟩⟨1|+0.44|1⟩⟨2| and evolves into |0⟩⟨1|, (d) and (e). The overtone (2ωe) coherence |0⟩⟨2| exists at early time. Its comparison with (f) clearly shows the measurement distortion due to the angular convolution. The system-bath coherence (ωb) does not have a reference state. Its time evolution at late time shows periodic momentum kicks the system receives from the local cage, which appear as the pure momentum state seen in the last frame. The line graphs are profiles along the indicated cuts on the experimental frames. The reference line profiles are along p = 0, except for (f), which is taken along a curved path to account for the angular convolution. Complete animation provided in Videos S2–S5 (enhanced online). (Video S2) (Video S3) (Video S4) (Video S5) [URL: http://dx.doi.org/10.1063/1.4813437.2] [URL: http://dx.doi.org/10.1063/1.4813437.3] [URL: http://dx.doi.org/10.1063/1.4813437.4] [URL: http://dx.doi.org/10.1063/1.4813437.5]

Close modal

The inseparable part of the system-bath superposition, |s, b⟩ + |s′, b′⟩, namely the entanglement, is seen as the system-bath coherence in phase space (see last column in Fig. 10 and Video S5). After an initial period where a negative trough is observed, the entanglement appears as the phase portrait of a pure momentum state (as in the last frame of Fig. 10) that alternates sign at every cycle (see Video S5). The interaction appears as periodic momentum exchanges that last throughout the observation time of 4 ps. The coupled motion completely preserves phase coherence. Given the separation of time scales between system and bath, ωe/ωb ∼ 5, the inseparable motion can be regarded as evolution of the base states of the system, i.e., system states that adiabatically follow a cage that rings throughout the B-state lifetime. This in turn relegates the requirement for preservation of the W-hole to phase coherent evolution of the cage-breathing mode. The observed entanglement |s, b⟩ + |s′, b′⟩ is reduced to identifying b with vibrations of the immediate cage, which in turn are coupled to the rest of the bath. The extent to which the cage mode induces coherence in the lattice is interesting to consider, since this would be the mechanism for nucleating a “cat”-state that may be macroscopically distinguishable,46 e.g., through thermal conductivity. That through entanglement the driven system may induce a more extended coherence in the bath did not materialize in recent semi-classical simulations of a closely related system: I2 in a Kr cluster.47 This can be shown to occur in the more special case of a superfluid helium bath, where adiabatic evolution of rotational states of a guest molecule can be observed to follow the circulation they induce in the droplet.48 Here, it may be argued that the guest and its immediate cage are well isolated from the rest of the lattice, and that it is therefore not too surprising to see quantum coherent evolution of an isolated system of two coupled modes. Their reduction to the model of two bilinearly coupled qubits in a quantum harmonic bath allows exact solutions in the spin-boson framework.49 For a sudden disturbance away from thermal equilibrium, the concurrence due to bath coherence, which is closely related to entanglement, very nearly approximates the time profile of the coherence we see in Fig. 9(c).50 Rather than quantifying the extent of such system-bath entanglement, here we are satisfied by characterizing its phase-space portrait.

We presented the complete quantum mechanical moving picture (see Video S1) of a molecular bond, captured at a rate of 1 frame every 5 fs, and at resolution limited by quantum uncertainty. The non-classical motion, as certified by the negative volume in its Wigner distribution function, was recorded in ordinary matter through a parametric four-wave mixing process in which a volume coherence on a scale of ∼103 μm3 is prepared and interrogated. The quantum coherent response of the ensemble allows it to be viewed as identical copies of one state — a pure state limited in span to the phase space of three system vibrations. Its tomographic reconstruction yields the complete phase-space description of a system undergoing coherent dissipation. Through the spectrally decomposed Wigner distribution functions, we observe the wavefunctions in which the system develops and the spatiotemporal coherence of the system-bath entanglement. The generality of the observation of arrested decoherence is underscored by its realization in rather ordinary matter — impure ice. We show that hidden behind the dephasing window, buried deep in the thermally accessible density of states of the medium, quantum coherent mechanics survives and can be observed with properly tailored measurements. Indeed, such coherence is implied in realizations of interferometric control in condensed matter, the most remarkable of which are demonstrations in the liquid phase.51 Our measurements demonstrate them through the direct observation of wavefunctions.

This research was completed under support of the NSF Center for Chemistry at the Space-Time Limit (CHE-082913). This work has benefited from fruitful discussions with J. A. Cina and N. Schwentner over the years. We thank G. Kerenskaya for preparing Figures 1–3.

1.
E.
Wigner
,
Phys. Rev.
40
,
749
(
1932
).
2.
M.
Hillery
,
R.
O’Connell
,
M.
Scully
, and
E.
Wigner
,
Phys. Rep.
106
,
121
(
1984
).
3.
W. B.
Case
,
Am. J. Phys.
76
,
937
(
2008
).
4.
A.
Kenfack
and
K.
Zyczkowski
,
J. Opt. B: Quantum Semiclassical Opt.
6
,
396
(
2004
).
5.
D.
Leibfried
,
T.
Pfau
, and
C.
Monroe
,
Phys. Today
51
,
22
(
1998
).
6.
K.
Vogel
and
H.
Risken
,
Phys. Rev. A
40
,
2847
(
1989
).
7.
W. P.
Schleich
,
Quantum Optics in Phase Space
(
Wiley-VCH
,
Berlin
,
2001
).
8.
W. P.
Schleich
and
M. G.
Raymer
,
J. Mod. Opt.
44
,
2021
(
1997
).
9.
D. T.
Smithey
,
M.
Beck
,
M. G.
Raymer
, and
A.
Faridani
,
Phys. Rev. Lett.
70
,
1244
(
1993
).
10.
A. I.
Lvovsky
,
H.
Hansen
,
T.
Alchele
,
O.
Benson
,
J.
Mlynek
, and
S.
Schiller
,
Phys. Rev. Lett.
87
,
050402
(
2001
).
11.
A. I.
Lvovsky
and
J.
Mlynek
,
Phys. Rev. Lett.
88
,
250401
(
2002
).
12.
P.
Bertet
,
A.
Auffeves
,
P.
Maioli
,
S.
Osnaghi
,
T.
Meunier
,
M.
Brune
,
J. M.
Raimond
, and
S.
Haroche
,
Phys. Rev. Lett.
89
,
200402
(
2002
).
13.
D.
Leibfried
,
D. M.
Meekhof
,
B. E.
King
,
C.
Monroe
,
W. M.
Itano
, and
D. J.
Wineland
,
Phys. Rev. Lett.
77
,
4281
(
1996
).
14.
C.
Kurtsiefer
,
T.
Pfau
, and
J.
Mlynek
,
Nature (London)
386
,
150
(
1997
).
15.
T. J.
Dunn
,
I. A.
Walmsley
, and
S.
Mukamel
,
Phys. Rev. Lett.
74
,
884
(
1995
).
16.
R. L.
Hudson
,
Rep. Math. Phys.
6
,
249
(
1974
).
17.
C.
Monroe
,
D. M.
Meekhof
,
B. E.
King
, and
D. J.
Wineland
,
Science
272
,
1131
(
1996
).
18.
I. U.
Goldschleger
,
V.
Senekerimian
,
M. S.
Krage
,
H.
Seferyan
,
K. C.
Janda
, and
V. A.
Apkarian
,
J. Chem. Phys.
124
,
204507
(
2006
).
19.
W. H.
Zurek
,
Rev. Mod. Phys.
75
,
715
(
2003
).
20.
C. A.
Angell
,
Science
267
,
1924
(
1995
).
21.
E. D.
Sloan
,
Clathrate Hydrates of Natural Gases
(
Marcel Dekker
,
New York
,
1998
).
22.
G.
Kerenskaya
,
I. U.
Goldschleger
,
V. A.
Apkarian
, and
K. C.
Janda
,
J. Phys. Chem. A
110
,
13792
(
2006
).
23.
I. U.
Goldshleger
,
G.
Kerenskaya
,
K. C.
Janda
, and
V. A.
Apkarian
,
J. Phys. Chem. A
112
,
787
(
2008
).
24.
E. T.
Branigan
,
N.
Halberstadt
, and
V. A.
Apkarian
,
J. Chem. Phys.
134
,
174503
(
2011
).
25.
M. I.
Bernal Uruchurtu
,
G.
Kerenskaya
, and
K. C.
Janda
,
Int. Rev. Phys. Chem.
28
,
223
(
2009
).
26.
A. C.
Legon
,
J. M. A.
Thumwood
, and
E. R.
Waclawik
,
Chem.- Eur. J.
8
,
940
(
2002
).
27.
I. U.
Goldshleger
,
G.
Kerenskaya
,
V.
Senekerimian
, and
K. C.
Janda
,
Phys. Chem. Chem. Phys.
10
,
7226
(
2008
).
28.
S.
Segale
and
V. A.
Apkarian
,
J. Chem. Phys.
135
,
024203
(
2011
).
29.
G.
Herzberg
,
Molecular Spectra and Molecular Structure: I. Constants of Diatomic Molecules
(
Van Nostrand Reinhold Company
,
New York
,
1950
).
30.
E. C. M.
Chen
and
W. E.
Wentworth
,
J. Phys. Chem.
89
,
4099
(
1985
).
31.
A ramp filter is applied to eliminate fringes that arise by the sharp cutoff in (3), which otherwise are amplified by the division in (4).
32.
U.
Janicke
and
M.
Wilkens
,
J. Mod. Opt.
42
,
2183
(
1995
).
33.
The running cat is the W-transform of the superposition of two 2D Gaussians, a|G+⟩ + b|G⟩, where |G+⟩ = G(q+q0,p-p0) and |G⟩ = G(q-q0, p-p0) with p0 ≠ 0.
34.
D.
Segale
,
M.
Karavitis
,
E.
Fredj
, and
V. A.
Apkarian
,
J. Chem. Phys.
122
,
111104
(
2005
).
35.
Z.
Li
,
R.
Zadoyan
, and
V. A.
Apkarian
,
J. Phys. Chem.
99
,
7453
(
1995
).
36.
N.
Pugliano
,
D. K.
Palit
,
A. Z.
Szarka
, and
R. M.
Hochstrasser
,
J. Chem. Phys.
99
,
7273
(
1993
).
37.
Z.
Li
,
J-Y.
Fang
, and
C. C.
Martens
,
J. Chem. Phys.
104
,
6919
(
1996
).
38.
The time shift operation takes advantage of the identities: cos(ωt) + cos(ωt+π) = 0, cos(ωt) − cos(ωt+π) = 2cos(ωt), and cos(ωt) ± cos(ωt+ϕ) = A cos(ωt+ϕ/2).
39.
M.
Guhr
and
N.
Schwentner
,
Phys. Chem. Chem. Phys.
7
,
760
(
2005
).
40.
M. M.
Koza
,
Phys. Rev. B
78
,
064303
(
2008
).
41.
J. S.
Tse
,
V. P.
Shpakov
,
V. R.
Belosludov
,
F.
Trouw
,
Y. P.
Handa
, and
W.
Press
,
Europhys. Lett.
54
,
354
(
2001
).
42.
J. S.
Tse
and
Z. Li
Uehara
,
Europhys. Lett.
56
,
261
(
2001
).
43.
M.
Celli
,
D.
Colognesi
,
L.
Ulivi
,
M.
Zoppi
, and
A. J.
Ramirez-Cuesta
,
J. Phys.: Conf. Ser.
340
,
012051
(
2012
).
44.
J. A.
Cina
,
Annu. Rev. Phys. Chem.
59
,
319
(
2008
).
45.
M. G.
Raymer
,
Contemp. Phys.
38
,
343
(
1997
).
46.
A. J.
Leggett
,
J. Phys. Condens. Matter
14
,
R415
(
2002
).
47.
M.
Buchholz
,
C-M.
Goletz
,
F.
Grossmann
,
B.
Schmidt
,
J.
Heyda
, and
P.
Jungwirth
,
J. Phys. Chem.
116
,
11199
(
2012
).
48.
M. N.
van Staveren
and
V. A.
Apkarian
,
J. Chem. Phys.
133
,
054506
(
2010
).
49.
A. J.
Leggett
,
S.
Chakravarty
,
A. T.
Dorsey
,
M. P. A.
Fisher
,
A.
Garg
, and
W.
Zwerger
,
Rev. Mod. Phys.
59
,
1
(
1987
).
50.
A. G.
Dijkstra
and
Y.
Tanimura
,
Phys. Rev. Lett.
104
,
250401
(
2010
).
51.
T.
Brixner
and
G.
Gerber
,
ChemPhysChem
4
,
418
(
2003
).