text:
# Dynamic Fano Resonance of Quasienergy Excitons in Superlattices

## Abstract

The dynamic Fano resonance (DFR) between discrete quasienergy excitons and sidebands of their ionization continua is predicted and investigated in dc- and ac-driven semiconductor superlattices. This DFR, well controlled by the ac field, delocalizes the excitons and opens an intrinsic decay channel in nonlinear four-wave mixing signals.

The Fano resonance (FR) results from quantum coupling between a discrete state and a degenerate continuum of states and manifests itself in optical spectra as asymmetric lineshape. FR was observed in a variety of atomic and molecular systems. Recently, this phenomenon has been reported in semiconductor quantum wells , in biased semiconductor superlattices (SSL’s) , and in bulk GaAs in the presence of magnetic field  where discrete excitons couple to continua of the lower transitions through Coulomb interaction . In addition to the frequency-domain experiments, transient four-wave mixing (FWM) in magneto-confined bulk semiconductors and in SSL’s has demonstrated that the FR is a fundamental type of dephasing mechanism in irreversible quantum dynamical processes.

An intense ac field, e.g. the free electron laser, can induce the dynamical localization (DL)  of electrons in SSL’s. The DL is directly related to the collapse of the quasienergy band of Floquet states , the temporal analogue to the Bloch states in spatially periodic potential . Research taking into account Coulomb interaction has shown that the DL may cause the dimension crossover of excitons  and the enhancement of exciton binding and opacity in linear  and nonlinear spectra . Meanwhile the dynamical delocalization (DDL) effect of an ac field is also investigated on localized Wannier-Stark (WS) states  in biased SSL’s , which is hindered by Coulomb binding .

In this Letter, we predict the dynamic Fano resonance (DFR) between quasienergy exciton states in biased SSL’s driven by an intense THz-field. In particular, the novelty of this DFR consists in the coupling between the discrete quasienergy exciton and the neighbor sideband of its ionization continuum. This DFR will also lead to the dynamical ionization of bound excitons in SSL’s, which itself is a new effect of the ac field.

The one-dimensional tight-binding model  is adopted in the present investigation. With excluding the realistic three-dimensional excitonic motion, the results are hardly compared to experimental data quantitatively, but we believe this simple model does present qualitatively correct results, in view of the fact that the ground WS exciton-state usually possesses much larger oscillator strength than its ionization continuum of in-plane motion, and the coupling induced by the ac field between the in-plane motion associated with different WS states is negligible because of their approximate orthogonality. The Coulomb interaction is treated on the basis of Hartree-Fock approximation. The neglected genuine many-body correlation  may lead to almost instantaneous decay in time-integrated (TI) FWM signals, which is irrelevant to the intrinsic dephasing . This effect is crucial for the magneto-confined FR  due to the magnetic field-induced enhancement of many-body correlation , but is less important for the decay of TI-FWM signals in biased SSL’s, which is, in essence, an intrinsic dephasing process .

The semiconductor Bloch equations in the presence of electric fields  for the exciton amplitude $`p_ka_{\mathrm{h}k}a_{\mathrm{e}k}`$ and electron-hole (e-h) pair density $`f_ka_{\mathrm{e}k}^{}a_{\mathrm{e}k}=a_{\mathrm{h}k}^{}a_{\mathrm{h}k}`$ read

$`_tp_k=`$ $`i\left[E_0+\epsilon _k(t)\right]p_k+\left[\omega _{\mathrm{BO}}+\omega _1\mathrm{cos}(\omega t)\right]_kp_k`$ (2)
$`+i\chi ^R(t)(12f_k)\gamma _2p_k,`$
$`_tf_k=`$ $`\left[\omega _{\mathrm{BO}}+\omega _1\mathrm{cos}(\omega t)\right]_kf_k2\mathrm{}\left\{\chi ^R(t)p_k^{}\right\},`$ (3)
where $`a_{\mathrm{e}k}`$ ($`a_{\mathrm{h}k}`$) is the electron (hole) annihilation operator for the quasi-momentum state $`|\mathrm{e}(\mathrm{h}),k`$, $`E_0`$ is the separation between centers of the conduction and the valence minibands, $`\epsilon _k(t)=\epsilon _k2_k^{}V_{k,k^{}}f_k^{}`$ is the renormalized combined e-h miniband dispersion with $`\epsilon _k=\frac{\mathrm{\Delta }}{2}\mathrm{cos}k`$, $`\omega _{\mathrm{BO}}`$ is the Bloch oscillation frequency due to the dc-field , $`\omega _1`$ is the strength of the driving ac field with frequency $`\omega `$, $`\chi ^R(t)=\chi (t)+_k^{}V_{k,k^{}}p_k^{}`$ is the renormalized strength of a near infrared optical excitation (here the dipole element is assumed $`k`$-independent), and $`V_{k,k^{}}`$ comes from the Coulomb potential.

The Schrödinger equation for the exciton wave function is just the linear part of Eq. (2) with the excitation and dephasing terms removed. The quasienergy is obtained by numerically diagonalizing the propagator

$$U(T+t,t)\widehat{\mathrm{T}}\mathrm{exp}\left(i_t^{T+t}H(t)𝑑t\right)$$

\[$`H(t)`$ is the Hamiltonian, $`T2\pi /\omega `$\], according to the secular equation

$$U(T+t,t)u_q(k,t)=\mathrm{exp}\left(i\epsilon _qT\right)u_q(k,t),$$

where $`u_q(k,t)`$ is the eigenstate with quasienergy $`\epsilon _q`$. Obviously, for any integer number $`m`$, $`u_{q,m}u_q\mathrm{exp}(im\omega t)`$ is also a solution with quasienergy $`\epsilon _q+m\omega `$, and, in analogy to the conventional Brillouin zone in quasi-momentum space, so is defined the Brillouin zone of quasienergy $`\left\{[m\omega ,(m+1)\omega ]\right\}`$.

In the ac-driven system, the absorption spectrum $`\alpha (\mathrm{\Omega })\mathrm{}\left[P^{(1)}(\mathrm{\Omega })/\chi (\mathrm{\Omega })\right]`$, where $`P^{(1)}(\mathrm{\Omega })`$ is the linear term of the Fourier transformation of the optical response $`P(t)_kp_k`$, and $`\chi (\mathrm{\Omega })`$ is the excitation profile . In this work, only the continuous wave absorption is considered, whereby the interband susceptibility and oscillator strength are of the same form as in the static case, except that the transition matrix element is replaced by its time-average and the transition energy by the quasienergy . The FWM signals are calculated following Ref. REFERENCES, in which two degenerate Gaussian-shaped pulses \[$`\chi (t)=\chi \mathrm{exp}(i\mathrm{\Omega }_0tG^2t^2/2)`$\] are separated by a delay time $`\tau `$.

For the sake of simplicity, let us focus on a specific case: The ac field is resonant with the Bloch oscillation, and the Coulomb potential is of on-site type, i.e. $`V_{k,k^{}}=V_0/N`$ ($`N`$ is the size of the basis set $`\{|k\}`$). Parameters for our calculations are taken as follows. $`V_0=10`$ meV, $`\mathrm{\Delta }=4V_0`$, $`\omega _{\mathrm{BO}}=\omega =1.5V_0`$, $`\gamma _2=0.1V_0`$, $`\mathrm{\Omega }_0=E_0V_0`$, and $`G=\omega _{\mathrm{BO}}`$. $`N`$ is chosen to be 40 for calculating quasienergy, and 80 for optical spectra. Our calculation shows that a larger basis produces no significant difference.

As shown in Fig. 1, the quasienergy spectrum consists of continuous minibands and a few well-separated discrete states below. As the ac field strengthens, owing to the DDL, the miniband broadens and the discrete excited states  merge into the continuum one by one. Meanwhile, the discrete ground state is repelled by its continuum towards lower energy, and eventually dips from above into the ionization continuum of its neighbor sideband, causing a series of anticrossing in the spectrum. Further enhancement of the ac field may suppress the quasienergy miniband, then the discrete exciton states can be released from the miniband one by one. At certain strength, namely $`\omega _1=\omega _{\mathrm{DL1}}`$, the miniband collapses and DL takes place. After that the evolution described above will repeat with increasing the ac field.

For an intuitive understanding, let us briefly review the picture of dressed WS ladders introduced in Ref. REFERENCES. In weak ac field limit, only one-photon-assisted hopping between neighbor WS states $`|n_\mathrm{x}`$ takes effect , thus the time-periodic states $`\{\mathrm{exp}(in\omega t)|n_\mathrm{x}\}`$ form the Wannier basis set in a tight-binding ‘lattice’ with the nearest-neighbor hopping coefficient of $`\omega _1\mathrm{\Delta }/(8\omega _{\mathrm{BO}})`$ . Those Wannier states with small $`n`$ can be viewed as ‘impurities’ in a crystal since their ‘on-site energy’ $`E_n^\mathrm{x}n\omega `$ deviates from zero significantly due to the Coulomb interaction . The discrete states in the quasienergy spectrum are concentrated on these ‘impurities’; on the other hand, the remote WS states are almost equally spaced and the resonant photon-assisted hopping results in the formation of the continuum as wide as without the Coulomb interaction. As the DDL is enhanced by the ac field, the discrete excited states may be ionized, and merge into the continuum. Calculation with this picture agrees well with the exact results for $`\omega _1<\omega /2`$ (not shown). However, when the ac field is strong enough, the inter-sideband interaction with multi-photon processes has to be invoked to account for the mixing between the discrete ground state and the continuum, the band suppression, and the DL.

The absorption spectrum in the absence of an ac field presents just the excitonic WS ladder (Fig. 2, $`\omega _1/\omega _{\mathrm{DL1}}=0`$). When the ac field is turned on, the WS states are coupled through photon-assisted hopping, so their sidebands start to gain in oscillator strengths , which looks as if the WS peaks were split (Fig. 2, $`\omega _1/\omega _{\mathrm{DL1}}=1/16`$). As the excited discrete states merge into the continuum, the intra-sideband DFR occurs, and the oscillator strength is shared among the continuum states which otherwise have negligible opacity. This manifests itself in the spectrum as broadening of the original discrete line (Fig. 2, $`\omega _1/\omega _{\mathrm{DL1}}=1/16`$ and $`1/8`$). When the ground exciton state couples to energetically degenerate continuum associated with its neighbor sideband, (for $`\omega _1/\omega _{\mathrm{DL1}}=1/8`$, 1/4, 1/2, and 3/4 in Fig. 2), the inter-sideband DFR happens, as characterized in the absorption spectra by asymmetric peaks with FWHM larger than $`\gamma _2`$ . As shown in the inset of Fig. 2, these peaks have perfect Fano lineshape

$$\alpha (\mathrm{\Omega })=\alpha _0+\alpha _\mathrm{c}\frac{(q\gamma +\mathrm{\Omega }E_\mathrm{x})^2}{\gamma ^2+(\mathrm{\Omega }E_\mathrm{x})^2},$$
(4)
where $`\alpha _0`$ is the background constant, $`\alpha _\mathrm{c}`$ represents the continuum absorption without the inter-sideband coupling, $`q`$ is the lineshape parameter, $`\gamma `$ is related to the resonance broadening, and $`E_\mathrm{x}`$ denotes the position of the discrete state. When further enhancement of the ac field reduces the quasienergy-band width and the discrete states emerge out of the continuum, the absorption spectrum evolves from the DFR to Lorentzian-shaped discrete lines.

The lineshape of the absorption peaks vary drastically when the ac field strength passes through the field at which the DFR occurs. As shown in Fig. 3, when the exciton ground state meets the continuous sideband, $`1/q`$ becomes nonzero, demonstrating the asymmetrical non-Lorentzian lineshape, and the Fano coupling parameter, $`\gamma `$, is larger than the static dephasing rate of excitons, $`\gamma _2`$, indicating that the peak is further broadened by the DFR. According to Ref. REFERENCES,

$$\gamma =\pi |V_E|^2$$

(in the present calculation, $`\gamma `$ also includes the contribution from $`\gamma _2`$), and

$$1/q^2=\left|\pi V_E\left(\mu _E/\mu _X\right)\right|^2,$$

where $`V_E`$ is the coupling matrix element between the discrete and the continuum states, and $`\mu _E`$ and $`\mu _X`$ are the optical transition matrix element of the continuum and the discrete states, respectively. $`V_E`$ is determined mainly by the overlap integral of the bound and extended excitons, which can be modulated by the ac field through changing the localization property of the excitons. Besides, the ac field redistributes the oscillator strength among the quasienergy states and their sidebands, changing the relative opacity of the continuum. Thus both $`1/q`$ and $`\gamma `$, depending on the sideband index, are readily tuned by the ac field \[see Fig. 3 (a) and (b)\].

The TI-FWM traces in Fig. 4 display the exponential decay superimposed with quantum beats, which result from the interference between the quasienergy excitons with different transition energy . The effective dephasing rate $`\gamma _2^{\mathrm{eff}}`$, extracted from Fig. 4 and plotted in Fig. 3 (c), follows approximately the Fano coupling parameter $`\gamma `$ \[Fig. 3 (b)\], indicating that the DFR plays a role in intrinsic phase-breaking process as the usual FR does in biased SSL’s.

To study the localization property of excitons in the dc- and ac-driven SSL’s, by integrating the Schrödinger equation in the Wannier representation, we have calculated the mean-square root of the exciton radius $`W\widehat{r}^2^{1/2}`$ and the probability $`R_0`$ of finding the electron and hole at the same site, at various delay time after a weak $`\delta `$-pulse excitation \[see Fig. 3 (d) and (e)\]. A $`\delta `$-pulse excites the e-h pair at the same site, so at $`t=0^+`$, $`W=0`$ and $`R_0=1`$. $`W`$ increases unboundedly with time as a result of the resonant photon-assisted hoping between the evenly spaced remote WS states, whenever the quasienergy miniband has finite width. This DDL effect is handicapped somewhat by the Coulomb binding: The $`R_0`$ may saturate with time whenever there remains a discrete state outside the continua. When the inter-sideband DFR occurs, however, the excitons are dynamically ionized with no such saturation behavior for $`R_0`$.

In summary, the dynamic Fano resonance between a discrete exciton and the sideband of ionization continuum is predicted in dc- and ac-driven SSL’s, which, well controlled by the ac field strength, manifests itself as broadened asymmetric lineshape in absorption spectra and intrinsic decay in FWM signals, and leads to the dynamic ionization for bound excitons. All these effects stem essentially from the unequal spacing of the excitonic WS ladder together with the DDL effect of the ac field. This DFR effect was absent in the absorption spectra calculated in Ref. REFERENCES, because the sidebands of quasienergy miniband are not effectively overlapped for the chosen parameters . It should be pointed out that the DFR is quite different in nature from the FR between WS states and the degenerate in-plane continua associated with low-lying WS states in dc-biased SSL’s . The DFR is expected to be experimentally distinguishable from the in-plane FR in its ac field dependence. Moreover, it’s possible to adjust the dc-field strength to minimize the in-plane FR effect  so as to emphasize the dynamic interference.

This work was supported by the National Science Foundation of China.


format: text

text:
# Extraction of work from a single thermal bath in the quantum regime

## Abstract

The stationary state of a quantum particle strongly coupled to a quantum thermal bath is known to be non-gibbsian, due to entanglement with the bath. For harmonic potentials, where the system can be described by effective temperatures, thermodynamic relations are shown to take a generalized Gibbsian form, that may violate the Clausius inequality. For the weakly-anharmonic case a Fokker-Planck type description is constructed. It is shown that then work can be extracted from the bath by cyclic variation of a parameter. These apparent violations of the second law are the consequence of quantum coherence in the presence of the slightly off-equilibrium nature of the bath.

The laws of thermodynamics are at the basis of our understanding of nature, so we all expect them to govern also systems coupled to a bath in the quantum regime. However, recently thought-provoking claims were made about a violation of Thomson’s formulation of the second law (the impossibility to do work periodically without loosing heat) and even about a perpetuum mobile acting in an inhomogeneous superconducting ring .

Most of our thermodynamic understanding is based on the Gibbs distribution. The laws of equilibrium thermodynamics apply equally well to closed classical and quantum systems, as to open classical subsystems . The setting for the classical case is well known: Under general statistical conditions  one derives a Langevin equation. The corresponding probability distribution is described by the Fokker-Planck equation, and it converges in time to the Gibbs distribution.

Much less is known about the quantum Langevin equation . The stationary distribution has only been obtained for the harmonic potential. It depends explicitly on the damping constant and becomes Gibbsian only in the limit of weak damping , thus preventing the applicability of equilibrium thermodynamics. Entanglement is the very reason of this crucial difference, as subsystems are necessarily in a mixed state.

In the present paper we examine the standard model for quantum Brownian motion, the so-called Caldeira-Leggett model , see eq. (1). Hereto we employ methods developed recently for glasses . For a particle in harmonic potential we define effective temperatures, and put the thermodynamic relations in a generalized Gibbsian form. For weakly anharmonic confining potentials Fokker-Planck equations will be constructed, which allow to obtain the stationary distribution and elucidate important aspects of nonstationary properties.

We shall provide a nontrivial thermodynamic interpretation for the relaxation towards the steady non-gibbsian state and for the slow change of a system parameter. Our main results are rather dramatic, apparently contradicting the second law: We show that the Clausius inequality $`\mathrm{¯}\mathrm{d}QT\mathrm{d}S`$ can be violated, and that it is even possible to extract work from the bath by cyclic variations of a parameter (“perpetuum mobile”). The physical cause for this appalling behavior will be traced back to quantum coherence in the presence of the near-equilibrium bath.

The Quantum Langevin equation is derived from the exact Hamiltonian description of a particle and a thermal bath, when tracing out the degrees of freedom of the bath. For $`t<0`$ the particle and bath do not interact, and the whole system is described by a density matrix $`\rho _0=\rho _b\rho _s`$, where $`\rho _b`$ is the Gibbs distribution for the bath, and $`\rho _s`$ describes the state of the particle. At $`t=0`$ a linear coupling is switched on instantaneously, and the total Hamiltonian reads for $`t>0`$

$`={\displaystyle \frac{p^2}{2m}}+V(x)+{\displaystyle \underset{i}{}}\left[{\displaystyle \frac{p_i^2}{2m_i}}+{\displaystyle \frac{m_i\omega _i^2}{2}}(x_i{\displaystyle \frac{c_ix}{m_i\omega _i^2}})^2\right]`$ (1)
where $`p`$, $`p_i`$, $`x`$, $`x_i`$, and $`m`$, $`m_i`$ are, respectively, the momentum and coordinate operators and the masses of the particle and the modes of the bath. The latter have a constant frequency gap $`\mathrm{\Delta }`$, so $`\omega _i=i\mathrm{\Delta }`$. For the couplings we choose the Drude-Ullersma spectrum

$$c_i=\sqrt{\frac{2\gamma m_i\omega _i^2\mathrm{\Delta }}{\pi }\frac{\mathrm{\Gamma }^2}{\omega _i^2+\mathrm{\Gamma }^2}},$$
(2)
where $`\gamma `$ is the damping constant, quantifying the stength of interaction, and $`\mathrm{\Gamma }`$ is the cut-off frequency of the interaction with the bath. The thermodynamic limit is taken for the bath by sending $`\mathrm{\Delta }0`$. Notice that this creates an infinite timescale $`1/\mathrm{\Delta }`$, implying that in the remaining approach the limit of “large times” always means the state where time is still much less than $`1/\mathrm{\Delta }`$.

For fast switching on the interaction, an amount of work

$$𝒲_0=\mathrm{tr}\{\rho _0[(0^+)(0^{})]\}=\frac{\gamma \mathrm{\Gamma }}{2}\mathrm{tr}(\rho _sx^2)$$
(3)
has to be supplied to the system. The large internal energy of the whole system $`U=p^2_0/2m+V_0+\pi ^2T^2/(6\mathrm{}\mathrm{\Delta })`$$``$$`1/\mathrm{\Delta }`$ is increased by the finite amount $`𝒲_0`$. This brings the system slightly out of equilibrium, since the von Neumann entropy $`\mathrm{tr}(\rho _0\mathrm{ln}\rho _0)`$ is exactly conserved under the switching. Notice that the switching energy $`𝒲_0=\frac{1}{2}\gamma \mathrm{\Gamma }x^2_0`$ is purely a classical effect.

The resulting quantum Langevin equation reads

$`\dot{p}+{\displaystyle \frac{\gamma \mathrm{\Gamma }}{m}}{\displaystyle _0^t}dt^{}e^{\mathrm{\Gamma }(tt^{})}p(t^{})+V^{}(x)=\eta (t)`$ (4)
with Gaussian noise having $`\eta =0`$ and anticommutator

$`K(t)={\displaystyle \frac{\{\eta (t),\eta (0)\}}{2}}={\displaystyle \frac{\mathrm{d}\omega }{2\pi }\frac{\gamma \mathrm{}\omega \mathrm{coth}(\frac{1}{2}\beta \mathrm{}\omega )e^{i\omega t}}{1+(\omega /\mathrm{\Gamma })^2}}`$ (5)
The connection between properties of the noise and the friction kernel is the consequence of quantum fluctuation-dissipation theorem . Eq. (4) with physically suitable forms of the potential and friction describes a rich variety of physical phenomena, such as Josephson junctions , processes in plasma and condensed matter , interaction of atoms with black-body radiation , and the Lamb shift of an electron .

We shall restrict ourselves to the quasi-Ohmic case, where $`\mathrm{\Gamma }`$ is much larger than other characteristic times. For $`t1/\mathrm{\Gamma }`$ we may expand the memory kernel in Eq. (4)

$`\mathrm{\Gamma }{\displaystyle _0^t}dt^{}e^{\mathrm{\Gamma }(tt^{})}p(t^{})=p(t){\displaystyle \frac{\dot{p}(t)}{\mathrm{\Gamma }}}`$ (6)
Compared to the classical white noise case, one can thus adjust only the noise, see eq. (5), while keeping the friction instantaneous. The same conclusion was reached in refs.  for the exactly solvable harmonic potential.

Work and heat. Let us explain the ingredients for the thermodynamic description of the Brownian particle having Hamiltonian $`H=p^2/2m+V(x)`$ and Wigner function $`W(p,x,t)`$. A change with time of the mean energy $`U=HdpdxWH`$ is considered when slowly varying a system parameter $`\alpha `$, such as $`m`$ or $`V^{\prime \prime }(0)`$, according to a prescribed trajectory $`\alpha (t)`$:

$$\mathrm{d}U\mathrm{d}H=H\mathrm{d}\mathrm{ln}W+\mathrm{d}H$$
(7)
The last term is the averaged mechanical work $`\mathrm{¯}\mathrm{d}𝒲`$ produced by external sources . For a variation of an intrinsic parameter (for which $`\mathrm{d}H=\mathrm{d}`$), $`\mathrm{¯}\mathrm{d}𝒲`$ is the change of the total (particle$`+`$bath) mean energy. The first term in r.h.s. is due to the statistical redistribution of the phase space. We shall identify it with the change in heat $`\mathrm{¯}\mathrm{d}𝒬=H\mathrm{d}\mathrm{ln}W`$, so Eq. (7) is just the first law.

In the harmonic case $`V(x)=\frac{1}{2}ax^2`$, Eqs. (4,6) can be integrated directly, see e.g. . In particular, the stationary Wigner distribution reads

$`W_s(p,x)={\displaystyle \frac{1}{2\pi }}\sqrt{{\displaystyle \frac{a}{mT_pT_x}}}\mathrm{exp}[{\displaystyle \frac{p^2}{2mT_p}}{\displaystyle \frac{ax^2}{2T_x}}].`$ (8)
The effective temperatures $`T_p`$ and $`T_x`$ approach $`T`$ in the classical limit, and have for large damping, $`\gamma ^2am`$ the following values at $`T0`$:

$`T_p={\displaystyle \frac{\mathrm{}\gamma }{\pi m}}\mathrm{ln}{\displaystyle \frac{\mathrm{\Gamma }m}{\gamma }}+{\displaystyle \frac{\mathrm{}a}{\pi \gamma }},T_x={\displaystyle \frac{\mathrm{}a}{\pi \gamma }}\mathrm{ln}{\displaystyle \frac{\gamma ^2}{am}}`$ (9)
The brownian particle has semiclassical behavior due to its interaction with the bath, and entropy $`S=S_p+S_x`$ with $`S_p=\frac{1}{2}\mathrm{ln}emT_p/\mathrm{}`$, $`S_x=\frac{1}{2}\mathrm{ln}eT_x/a\mathrm{}`$.

For an adiabatic variation of a parameter $`\alpha `$ the situation is still described by Eq. (8) with $`\alpha =\alpha (t)`$. One can check for $`U`$ and the free energy $`F=UT_pS_pT_xS_x`$

$`\mathrm{d}U=\mathrm{¯}\mathrm{d}Q_{ad}+\mathrm{¯}\mathrm{d}𝒲_{ad}=T_p\mathrm{d}S_p+T_x\mathrm{d}S_x+\mathrm{¯}\mathrm{d}𝒲_{ad},`$ (10)
$`\mathrm{d}F=S_x\mathrm{d}T_xS_p\mathrm{d}T_p+\mathrm{¯}\mathrm{d}𝒲_{ad}.`$ (11)
These generalized thermodynamical relations are in close analogy with those proposed for glassy systems .

To demonstrate a violation of the fundamental Clausius inequality $`\mathrm{d}QT\mathrm{d}S`$ we consider a slow variation of the mass $`m`$ at $`T0`$. Using Eqs. (8, 9, 10) one gets

$`\mathrm{¯}\mathrm{d}𝒲_{ad}=T_p{\displaystyle \frac{\mathrm{d}m}{2m}},\mathrm{¯}\mathrm{d}𝒬_{ad}=\mathrm{d}U\mathrm{¯}\mathrm{d}𝒲_{ad}={\displaystyle \frac{\mathrm{}\gamma }{m^2\pi }}{\displaystyle \frac{\mathrm{d}m}{2}}`$ (12)
Thus, there is the transfer of heat even for $`T=0`$, and the Clausius inequality is violated when $`\mathrm{¯}\mathrm{d}Q>0`$, i.e. for $`\mathrm{d}m>0`$. The latter also happens when varying $`a`$.

The Fokker-Planck equation. Except for the solvable harmonic potential, Eq. (4) is hardly tractable. Another approach having started from (4), goes further to the weak-coupling (small $`\gamma `$) limit described by a Markovian master-equation . In the opposite, strong-coupling limit one can consider the terms in (4) as c-numbers, but with the Gaussian quantum noise . The correspondence with the underlying quantum problem is established through the Wigner function, i.e., $`\delta (p(t)p)\delta (x(t)x)=W(p,x,t)`$. Since this approach is still exact for the harmonic case, one of the conditions of its validity involves a characteristic scale $`L`$ where the non-linearity remains small: $`L\sqrt{\mathrm{}/\gamma }`$ . Based on this condition we have derived a closed equation for the Wigner function. Here we only give the final result, while details will be presented elsewhere ,

$`{\displaystyle \frac{W(x,p,t)}{t}}={\displaystyle \frac{p}{m}}{\displaystyle \frac{W}{x}}+{\displaystyle \frac{}{p}}\left([{\displaystyle \frac{\gamma }{m}}p+V^{}(x)]W\right)`$ (13)
$`+{\displaystyle \frac{^2}{px}}\left(D_{xp}(x,t)W\right)+\gamma D_{pp}(x,t){\displaystyle \frac{^2W}{p^2}},`$ (14)
For our purposes it is enough to indicate the stationary values of the diffusion coefficients $`D_{xp}`$ and $`D_{pp}`$,

$`D_{pp}(x)={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\mathrm{d}\omega }{\pi m}}{\displaystyle \frac{\overline{K}(\omega )\omega ^2}{(\omega ^2+\omega _1^2)(\omega ^2+\omega _2^2)}}`$ (15)
$`D_{xp}(x)={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\mathrm{d}\omega }{\pi m^2}}{\displaystyle \frac{\overline{K}(\omega )(V^{\prime \prime }(x)m\omega ^2)}{(\omega ^2+\omega _1^2)(\omega ^2+\omega _2^2)}}`$ (16)
where $`\overline{K}(\omega )`$ is the spectrum of $`K(t)`$ in Eq. (5) and $`\omega _{1,2}=[\gamma \pm \sqrt{\gamma ^24mV^{\prime \prime }(x)}]/(2m)`$. Equation (13) was derived under the same assumptions as the semiclassical Langevin equation itself. For ensuring the convergence of the diffusion coefficients ($`\mathrm{Re}(\omega _{1,2})0`$) we will demand $`V^{\prime \prime }0`$ (local stability). In the classical limit: $`D_{xp}0`$, $`D_{pp}T`$, and Eq. (13) tends to usual Fokker-Planck equation . For the harmonic case $`D_{xp}`$ and $`D_{pp}`$ become space-independent, and Eq. (13) is in agreement with previous results , obtained for a more general type of the environment.

In the physically interesting case of overdamped motion, where the characteristic times of the momenta $`m/\gamma `$ and coordinate $`\gamma /V^{\prime \prime }(x)`$ are widely separated $`\gamma /V^{\prime \prime }(x)m/\gamma `$, a solution of the Fokker-Planck equation for $`tm/\gamma `$ can be presented as

$$W(p,x,t)=\frac{\mathrm{exp}[p^2/(2mD_{pp}(x,t))]}{\sqrt{2\pi mD_{pp}(x,t)}}W(x,t)$$
(17)
Here $`W(x,t)`$ is the solution of a reduced equation

$`_tW(x,t)+{\displaystyle \frac{J(x,t)}{x}}=0,`$ (18)
$`J={\displaystyle \frac{1}{\gamma }}V^{}(x)W(x,t){\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{}{x}}[D(x,t)W(x,t)],`$ (19)
$`D(x,t)=D_{xp}(x,t)+D_{pp}(x,t).`$ (20)
It still contains the inhomogeneous, time-dependent diffusion coefficient. As we discussed, large, but finite $`\mathrm{\Gamma }\omega _{1,2}`$ is necessary only for the statistics of the momenta, see Eq. (9). Taking $`\mathrm{\Gamma }\mathrm{}`$ in Eq. (20) we obtain

$`D(x)={\displaystyle \frac{\mathrm{}\gamma V^{\prime \prime }(x)}{\pi m^2}}{\displaystyle _0^{\mathrm{}}}d\omega {\displaystyle \frac{\omega \mathrm{coth}\left(\frac{1}{2}\beta \mathrm{}\omega \right)}{(\omega ^2+\omega _1^2)(\omega ^2+\omega _2^2)}}`$ (21)
In the stationary state one has $`J(x)=0`$, which implies for the corresponding distribution

$`W_s(x)={\displaystyle \frac{e^{\beta V_e(x)}}{Z}};V_e(x)=T{\displaystyle _0^x}dy{\displaystyle \frac{V^{}(y)+D^{}(y)}{D(y)}}`$ (22)
where $`V_e`$ is an effective potential $`V_e`$. Notice that $`W_s(x)`$ is non-gibbsian. For the harmonic potential Eq. (8) is recovered from this expression while the classical Gibbs distribution appears in the limit $`\mathrm{}\beta 0`$. Eqs. (17, 22) show that the statistics of momenta is influenced by the coordinate, yielding its non-Maxwellian form.

In classical case universal (thermodynamical) properties of the relaxation are described by an $``$-theorem, which is intimately connected with a formulation of the second law . This theorem can be generalized in our case, at least for times $`tm/\gamma `$, where the momenta already came to the local equilibrium, and the relevant variable is $`x`$. The $``$-function is defined as

$$=dxW_1(x,t)\mathrm{ln}\frac{W_1(x,t)}{W_2(x,t)}0,$$
(23)
where $`W_{1,2}(x,t)`$ are solutions of Eq. (18) corresponding to different initial conditions. Calculating $`\dot{}`$ from Eq. (18) one performs partial integrations to obtain

$`\dot{}={\displaystyle \frac{1}{\gamma }}{\displaystyle dxW_1(x,t)D(x,t)\left[\frac{}{x}\mathrm{ln}\frac{W_1(x,t)}{W_2(x,t)}\right]^2}`$ (24)
Since $`D(x,t)>0`$, $``$ is a monotonically decreasing function, and attains its minimum in the stationary state . On the other hand, $``$ is limited from below by zero, so we conclude that all solutions of the Fokker-Planck equation (18) converge with time to the stationary solution. Let us define the current-less state $`W_s(x,t)`$, obtained from $`W_s(x)`$ by the substitution $`D(x)D(x,t)`$. Starting from Eq. (24) one can show that the entropy of the coordinate sector obeys

$$\frac{\mathrm{d}S_x}{\mathrm{d}t}=dx\dot{W}(x,t)\mathrm{ln}W(x,t)=\frac{\mathrm{d}_eS_x}{\mathrm{d}t}+\frac{\mathrm{d}_iS_x}{\mathrm{d}t}.$$
(25)
Here the quantity

$`{\displaystyle \frac{\mathrm{d}_eS_x}{\mathrm{d}t}}={\displaystyle dx\dot{W}(x,t)\mathrm{ln}\frac{1}{W_s(x,t)}}={\displaystyle dxJ(x)\beta V_e^{}(x,t)}`$ (26)
can be interpreted as the flux of entropy. Indeed, $`J(x)`$ is the probability current and $`\beta V_e^{}`$ represents a force divided by temperature. The second term in (25),

$`{\displaystyle \frac{\mathrm{d}_iS_x}{\mathrm{d}t}}={\displaystyle \frac{1}{\gamma }}{\displaystyle dxW(x,t)D(x,t)\left[\frac{}{x}\mathrm{ln}\frac{W(x,t)}{W_s(x,t)}\right]^2}`$ (27)
is the entropy production during the relaxation. This quantity is strictly positive out of the steady state, and becomes zero in the long-time limit. In the classical case Eq. (25) leads to the well-known relation $`T\mathrm{d}S=\mathrm{d}𝒬+T\mathrm{d}_iS`$, since $`W_s(x,t)`$ becomes time-independent and coincides with the Gibbs distribution. However, in the general quantum case the relation between heat ($`\mathrm{¯}\mathrm{d}𝒬`$) and the flux of the statistical entropy ($`\mathrm{¯}\mathrm{d}𝒬=T\mathrm{d}_eS`$) appears to be broken. A similar relation can, though, still be recovered in the harmonic case, where the diffusion constant $`D`$ does not depend on the coordinate. Applying Eq. (25) we find $`T_p\mathrm{d}S_p+T_x\mathrm{d}S_x=\mathrm{d}U+\mathrm{¯}\mathrm{d}\mathrm{\Pi }`$, where $`T_p(t)=D_{pp}(t)`$ and $`T_x(t)=D(t)`$ are the time-dependent effective temperatures, and $`\mathrm{¯}\mathrm{d}\mathrm{\Pi }=T_x\mathrm{d}_iS_x0`$ is the energy dissipated during the relaxation.

Returning to the case of varying a system parameter, we are now interested in the first non-adiabatic correction to the stationary distribution, arising when the time of the variation is large, but finite. Since the dynamics is in the overdamped regime, the main correction comes from the deviation of $`W(x,t)`$ (and not the full $`W(p,x,t)`$) from its stationary form. The variation starts at $`t=t_i`$ (the particle already reached its stationary state), and ends at $`t=t_f`$. We shall assume a “smooth” start of the variation, i.e., $`\dot{\alpha }(t_i)=0`$. It can be checked directly from Eq. (18) that $`W`$ now reads for small $`\dot{\alpha }`$

$`W(x,t)=W_s(x,\alpha )\left[1+{\displaystyle \frac{\gamma }{T}}\dot{\alpha }\left(B(x,t)B\right)\right],`$ (28)
$`B(x,t)={\displaystyle _{\mathrm{}}^x}{\displaystyle \frac{\mathrm{d}y}{D(y,\alpha )W_s(y,\alpha )}}𝒜(y,{\displaystyle \frac{V_e}{\alpha }}),`$ (29)
$`𝒜(y,f)={\displaystyle _y^{\mathrm{}}}dzW_s(z,\alpha )(f(z)f),`$ (30)
where $`\alpha =\alpha (t)`$ and the average is taken w.r.t. $`W_s(x,\alpha )`$. The work finally becomes $`\mathrm{¯}\mathrm{d}𝒲=\mathrm{¯}\mathrm{d}𝒲_{ad}+\mathrm{¯}\mathrm{d}\mathrm{\Pi },`$ with

$`\mathrm{¯}\mathrm{d}𝒲_{ad}=\dot{\alpha }\mathrm{d}t{\displaystyle dxdpW_s(p,x,\alpha )\frac{H(p,x,\alpha )}{\alpha }},`$ (31)
$`\mathrm{¯}\mathrm{d}\mathrm{\Pi }={\displaystyle \frac{\gamma }{T}}\dot{\alpha }^2\mathrm{d}t{\displaystyle \frac{\mathrm{d}x𝒜(x,_\alpha V)𝒜(x,_\alpha V_e)}{D(x,\alpha )W_s(x,\alpha )}}`$ (32)
In the Gibbsian case at large $`T`$, where $`V_eV`$, one recovers $`𝒲_{ad}=UTS`$ and the known result for $`\mathrm{¯}\mathrm{d}\mathrm{\Pi }0`$ . For low $`T`$ $`\mathrm{¯}\mathrm{d}\mathrm{\Pi }`$ remains finite due to the explicit factor $`T`$ in $`V_e`$, see (22). We stress that $`\mathrm{¯}\mathrm{d}\mathrm{\Pi }`$ is always relevant for cyclic processes $`\alpha (t_i)=\alpha (t_f)`$, where $`\mathrm{\Delta }𝒲_{ad}=0`$.

In the harmonic case $`\alpha `$ can stand for the spring constant $`a`$. It is then straightforward to see that $`\mathrm{¯}\mathrm{d}\mathrm{\Pi }(a/T_x)/a`$. As expected, this is always positive.

For the oscillator with $`V(x)=ax^2/2+gx^4/12`$ with small $`g`$ the anharmonicity is displayed at the scale $`L=\sqrt{a/g}|x|`$. We shall now investigate the non-adiabatic correction to the work caused by temporal variation of $`L`$. For $`T0`$ we additionally take the limit of large $`\gamma `$ in Eqs. (22, 21), which yields

$`\beta V_e(x){\displaystyle \frac{\gamma \pi }{2\mathrm{}\mathrm{ln}\gamma }}\left[{\displaystyle \frac{x^2}{6}}+{\displaystyle \frac{L^2}{3}}\mathrm{ln}(1+{\displaystyle \frac{x^2}{L^2}})\right].`$ (33)
Although $`_LV(x)=[a/(6L^3)]x^4`$ is negative, $`_LV_e(x)+[\gamma \pi T/(6\mathrm{}L^3\mathrm{ln}\gamma )]x^4`$ is positive. Therefore $`\mathrm{¯}\mathrm{d}\mathrm{\Pi }`$ is negative. For small $`g`$ we may calculate the integral in Eq. (32) in the harmonic approximation, to get

$$\mathrm{\Delta }\mathrm{\Pi }=\frac{7}{3}\gamma x^2^3_{t_i}^{t_f}dt\frac{\dot{L}^2(t)}{L^6(t)},$$
(34)
where $`x^2=(2\mathrm{}\mathrm{ln}\gamma )/(\pi \gamma )`$ is the dispersion in this approximation. The possibility to extract this energy from the bath is due to its non-equilibrium state, which is ensured by the energy (3) supplied in the switching.

In conclusion, we have considered a brownian quantum particle strongly interacting with a quantum thermal bath. The non-gibbsian statistics of the particle is completely described by Fokker-Planck equations (13, 18). An $``$-theorem is formulated in eqs. (24), (25). For the harmonic potential generalized Gibbsian relations can be constructed in terms of effective temperatures (8), (9), as happens also in glassy systems . Two formulations of the second law, namely the Clausius inequality and the impossibility to extract work during cyclical variations, can be apparently violated at low temperatures. One could thus speak of a “perpetuum mobile of the second kind”. We should mention, however, that the number of cycles can be large, but not arbitrarily large. As a result, the total amount of extractable work is modest .

These violations of the second law are due to quantum coherence in the presence of the slightly off-equilibrium nature of the bath. We call them apparent violations, since, the standard requirements for a thermal bath not being fulfilled, thermodynamics just does not apply. Let us stress that also in the classical regime the harmonic oscillator bath is not in full equilibrium, but there the Gibbs distribution saves the day and thermodynamics does apply. Our results thus make clear that the characterization of the heat bath should be given with care. If it thermalizes on a short timescale, standard thermodynamics always applies. If it does not thermalize within the observation time, it acts as a “mechanical” part of the system, and thermodynamics need not have a say. Although it still applies in the classical case, we have discussed that it fails to do so in the quantum case.

Let us close by noting that in the harmonic case the unequal effective temperatures do not cause heat currents that equalize them. This situation is reminiscent of the classical paradox that atoms should radiate, but, being in the quantum regime, they do not.

We thank H. Dekker for useful discussions. A.E. A is grateful for support by FOM (The Netherlands).


text:
# On the existence of internal modes of sine-Gordon kinks

## I Introduction

Solitons, which were originally a concept arising in the study of integrable systems , have rapidly become a key paradigm in perturbed, non-integrable systems describing very many physical applications . A very succesful picture of soliton dynamics under different perturbations has resulted from the use of collective coordinate (CC) techniques : In many instances, solitons behave basically like point-like particles, and therefore their time evolution can be predicted by equations of motion which are ordinary differential equations. Nevertheless, along the last decade it has been progressively realized that internal degrees of freedom of the solitons play a crucial role in a number of problems: They govern resonant soliton (or solitary wave) collisions , give rise to nontrivial soliton-impurity interactions , and can be excited both by ac forces  and by thermal noise . Therefore, it is very important to assess whether or not physically relevant solitons possess internal modes.

In this Rapid Communication we study solitons of the sine-Gordon (sG) equation, which arises in a diverse range of areas of physics, covering from crystal dislocation theory to Josephson junctions . In spite of the fact that linear stability analysis leads to only one eigenvalue in the discrete spectrum, corresponding to a zero frequency Goldstone mode , Boesch and Willis  claimed that they found an internal “quasimode” above but close to the lower phonon band edge. They described this “quasimode” as a long-lived oscillation of the width of the sG kink. The possibility of such a mode had been suggested earlier by Rice , by means of a variational approach that reproduced with great accuracy the internal mode frequency of the $`\varphi ^4`$ model  and predicted the existence of a similar mode for the sG problem. To our knowledge, Ref.  is the only paper reporting the observation of this mode, although hints of its existence are scattered all over the soliton literature . A few years later, it was theoretically found , within the Rice approach, that the internal mode of the sG and the $`\varphi ^4`$ kinks could be excited by a constant external force or by an initial distortion; however, no numerical evidence was presented to confirm these results. Therefore, we decided to revisit the issue of the existence of the internal mode of the sG equation in order to give a definitive answer to it. Our approach, based in part on our earlier findings , overcomes the difficulty that the internal “quasimode” would lie within the phonon band by allowing us to look at driving frequencies well in the gap. As we will show below, this permits us to establish clearly that this putative internal mode does not exist. Furthermore, we do observe the excitation of a discreteness induced internal mode , whose frequency is in excellent agreement with recent predictions by Kevrekidis and Jones . We substantiate these claims by briefly recalling our main analytical results in Sec. II, collecting our numerical results in Sec. III, and discussing our conclusions in Sec. IV.

## II Resonances due to ac forces

We begin by considering the perturbed sG equation

$$\begin{array}{c}\varphi _{tt}\varphi _{xx}=\mathrm{sin}\varphi +ϵ\mathrm{sin}(\delta t+\delta _0).\end{array}$$
(1)
where $`ϵ`$, $`\delta `$, and $`\delta _0`$ are the amplitude, frequency, and phase of the external periodic force. In  we analyzed analytically and numerically the above problem in the more general context of nonlinear (possibly damped) Klein-Gordon systems, which includes Eq. (1) as well as the $`\varphi ^4`$ model. Particularizing our results, obtained by means of the generalized travelling wave Ansatz  combined with the Rice Ansatz , we find that the momentum $`P(t)`$ and the width of the kink $`l(t)`$ (of a general nonlinear Klein-Gordon equation; specific constants for the sG case will be given below) obey the equations

$`{\displaystyle \frac{dP}{dt}}`$ $`=`$ $`qϵ\mathrm{sin}(\delta t+\delta _0),P(t){\displaystyle \frac{M_0l_0\dot{X}}{l(t)}},`$ (2)
$`\alpha \left(\dot{l}^22l\ddot{l}\right)`$ $`=`$ $`{\displaystyle \frac{l^2}{l_0^2}}\left(1+{\displaystyle \frac{P^2}{M_0^2}}\right)1,`$ (3)
where X(t) is the kink position. For sG, $`q=2\pi `$, $`M_0=8`$, $`l_0=1`$ and $`\alpha =\pi ^2/12`$. The equation for the momentum is linear and can be solved exactly, yielding

$$P(t)=\left[P(0)\frac{qϵ\mathrm{cos}(\delta _0)}{\delta }\right]+\frac{qϵ}{\delta }\mathrm{cos}(\delta t+\delta _0)$$
(4)
We thus see that the term $`P(t)^2`$ in Eq. (3) involves two frequencies, $`\delta `$ and $`2\delta `$. As shown in , for the $`\varphi ^4`$ kink this leads to dramatic resonance effects when the ac driving frequency $`\delta `$ is half the frequency $`\mathrm{\Omega }_i`$ of the internal mode. Without resonant excitation of kink modes, the kink oscillates around its initial position as a pure point-like particle, as shown analytically and numerically in . However, as shown in  for the above parameters, in case an internal mode is excited as in the $`\varphi ^4`$ system, the system energy grows and eventually leads to chaotic motion of the kink. In addition, the resonance at $`\mathrm{\Omega }_i/2`$ turns out to be much stronger than that at $`\mathrm{\Omega }_i`$, which makes it a convenient way to probe those modes: In case the sG kink would possess such a degree of freedom, we should observe similar phenomena. In , the frequency of the numerically observed “quasimode” was $`\omega _s=1.004\pm 0.001`$ in our units, whereas Rice’s theoretical analysis leads to  $`\mathrm{\Omega }_R=1.103`$. In both cases, the corresponding half frequencies lie well within the gap; we can thus force the system at the predicted resonances without exciting much radiation which could mask any internal mode effects.

## III Numerical results

In this Section, we present the results of our numerical search for internal modes or “quasimodes” of the sG equation. We have computed the numerical solution of the perturbed sG equation (1) by using the Strauss-Vázquez scheme , the total length of the system being $`L=100`$, the steps being $`\mathrm{\Delta }x=0.05`$, $`\mathrm{\Delta }t=0.005`$ (or other values when indicated). Free boundary conditions were also used, and the final time in the simulations was $`\mathrm{25\hspace{0.17em}000}`$, enough to detect any possible resonance .

### A ac force

First, we have perturbed our system with an ac force of amplitude $`ϵ=0.01`$ and phase $`\delta _0=\pi /2`$. We have used a sG kink at rest as initial condition. We have computed the energy (as defined in ) during the evolution, and then the mean value of the energy in the time interval $`\mathrm{10\hspace{0.17em}000}t\mathrm{25\hspace{0.17em}000}`$, studying a range of values of $`\delta `$ around half of the possible internal mode frequencies. Figure 1 collects our results for the mean energy as a function of the frequency. We clearly see that the mean energy rapidly increases for some values of $`\delta `$ which do not coincide with any of the proposed internal modes or “quasimodes.” The first (leftmost) maximum lies below half the frequency of the lower phonon band edge; therefore, this one can not be associated with a phonon mode and we will discuss its origin in detail below. As for the other peaks, Table I shows that we have been able to identify them with half of the values of the first frequencies of the phonon band, $`\omega _n`$, which are given by

$$\omega _n=\sqrt{1+\left[(2\pi n)/L\right]^2},\text{ }n=1,2,3,\mathrm{},N=L/\mathrm{\Delta }x.$$
(5)
A comparison of the third and fourth columns of Table I, where we list the half values of the frequencies corresponding to the first $`8`$ radiational modes and the numerical peak frequencies $`\stackrel{~}{\omega }_n`$ in Fig. 1, counting from the second maximum from the left, allows to verify that these values are nearly identical in all cases (see also the relative differences $`|\omega _n/2\stackrel{~}{\omega }_n|`$ in Table I). We stress that identification of radiation modes as the origin of these resonances makes perfect physical sense: Indeed, an extended CC approach including radiation modes leads to the conclusion that the radiation modes should also be excited parametrically by ac drivings. In view of all this, the main conclusion we can draw from Fig. 1 is the absence of any resonance due to the internal quasi-mode in the phonon band reported by Boesch and Willis , which, in view of the frequency they reported, would lead to an extra peak located between the two rightmost ones in Fig. 1. In this regard, it is important to note that the accuracy of our simulations in detecting every radiation mode rules out the possibility that we have missed this quasi-mode (at either $`\omega _s`$ or $`\mathrm{\Omega }_R`$).

Interestingly, we do observe an internal mode: Let us recall that the first maximum of the energy occurs at $`\delta =\omega _d=0.4999`$, i.e., $`2\omega _d=0.9998`$ is below the lower phonon frequency $`\omega _{ph}=1`$ as we mentioned above. The reason for this peak is a novel internal mode of the sG kink, that appears as a bifurcation from the lowest phonon of the continuum equation, of frequency $`\omega _{ph}=1`$, due to the effect of discreteness . As has been recently shown , the frequency of this new internal mode is $`\omega _i\sqrt{1(4/2025)(\mathrm{\Delta }x)^4}`$ (the $`(\mathrm{\Delta }x)^4`$ dependence was already observed in ). In our simulations, $`\mathrm{\Delta }x=0.05`$, hence $`\omega _i=0.999998`$, approximately equal to $`2\omega _d`$. Therefore, the first peak in Fig. 1 corresponds to the parametric resonance at half the value of the discreteness induced internal mode, very much like the resonance with the $`\varphi ^4`$ kink internal mode reported earlier .

### B Deformed initial kink

In view of our negative result about the internal quasi-mode of the sG kink, we carried out more numerical simulations specifically designed to find it. Following , we chose initial conditions given by

$$\begin{array}{c}\varphi (x,0)=4\text{arctan}\left(\mathrm{exp}\left[\frac{xX(0)}{l(0)}\right]\right),\hfill \\ \\ \varphi _t(x,0)=2\frac{\left[{\displaystyle \frac{u(0)}{l(0)}}{\displaystyle \frac{xX(0)}{l^2(0)}}\dot{l}(0)\right]}{\text{cosh}\left[{\displaystyle \frac{xX(0)}{l(0)}}\right]},\hfill \end{array}$$
(6)
which correspond to a deformed kink if $`l(0)1`$ or $`\dot{l}(0)0`$. By starting their simulations with such initial conditions, Boesch and Willis claimed that they were able to excite the internal quasi-mode, and hence we hope that, if it is indeed present, we should find it.

At this point, we have to recall the work of Majerníková, Gaididei and Braun , who, as mentioned in the introduction, theoretically considered the problem of a deformed kink given by the expressions above driven by a external constant force. When applied to this problem, our CC approach yields the same two ordinary differential equations as in  for $`X(t)`$ and $`l(t)`$, \[Eqs. (2)-(3) without the sine term, i.e., with only a constant force $`ϵ`$\]. We have solved these equations exactly for arbitrary initial conditions  in terms of the Whittaker functions, concluding that the kink width oscillates only if $`l(0)l_s`$ (the natural width of the undistorted kink, $`1/\sqrt{1u(0)^2}`$) or $`\dot{l}(0)0`$ in (6), whereas external force and dissipation (if present) only damp out those oscillations. Hence, kink width oscillations should be easiest to detect without force and damping and, consequently, here we consider only free propagation of an initially distorted kink (a general analysis of the constant force problem will be presented elsewhere ). As a check of the validity of this procedure to probe internal modes, we have numerically verified that the $`\varphi ^4`$ kink indeed behaves in this fashion due to the excitation of the internal mode, exhibiting oscillations in the kink width very close to those predicted by the CC approach. We then turned to the sG case, and we found that the oscillations are not correctly described by the CC equations. In Fig. 2 we can see the evolution of the kink width computed from numerical simulation of the sG system, with $`l(0)=1`$, $`u(0)=0`$, and $`\dot{l}(0)=0.3`$. It is clear that $`l(t)`$ is not a simple oscillatory function as in the $`\varphi ^4`$ case, and Fourier analysis yields the following frequencies: $`\omega _d=0.9983`$, $`\stackrel{~}{\omega }_1=1.0034`$, $`\stackrel{~}{\omega }_2=1.0083`$, $`\stackrel{~}{\omega }_3=1.0184`$, $`\stackrel{~}{\omega }_4=1.0335`$, $`\stackrel{~}{\omega }_5=1.0511`$, $`\stackrel{~}{\omega }_6=1.0712`$, $`\stackrel{~}{\omega }_7=1.0963`$. As before, all these frequencies are very close to the discreteness induced internal mode and the first few radiational modes, with no evidence supporting the existence of the internal quasi-mode above the phonon band edge.

## IV Conclusions

The first conclusion of this work is that the existence of the proposed internal quasi-mode of the sG kink above the phonon band edge  is very unlikely. We have found neither resonances nor long-lived kink width oscillations arising from such a mode when the kink is subject to ac driving or to initial deformations, respectively. We are confident that we have carefully explored all the range of relevant frequencies, as we have accurately detected all possible resonances within that range. In addition, the comparison with the phenomenology observed for the $`\varphi ^4`$ kink, that possesses a true internal mode, shows that if the sG kink would have such an internal mode, we would have detected its influence through numerical simulations as for the $`\varphi ^4`$ equation . In fact, now that we have seen that radiation modes are parametrically excited by their corresponding half-frequency drivings, we believe we can suggest an explanation for the results of Boesch and Willis. Notice that they used Eq. (5) to calculate the lower phonon frequency, obtaining that $`\omega _1=\sqrt{1+(2\pi /1000)^2}=1.00002`$, whereas the internal “quasimode” frequency that they found in the numerical simulations was $`\omega _s=1.004`$. As those two values were largely different, they concluded that the quasi-mode was well above the lowest phonon frequency and different from it. However, as we have discussed in the previous sections, they used initial conditions in which the initial kink velocity was not zero; in this case, it can be easily shown that the corresponding frequencies of the phonons are given by

$$\overline{\omega }_k=\frac{\omega _kku(0)}{\sqrt{1u^2(0)}},\omega _k=\sqrt{1+k^2}.$$

If we now insert the parameters used in  to perform a numerical experiment similar to the one discussed here, we find $`\omega _1=1.004`$, i.e., the resonance observed by Boesch and Willis took place in fact with the lowest frequency phonon in the presence of a moving kink and not with any internal quasi-mode.

The second conclusion of the present paper relates to the recently proposed internal modes induced by perturbations, in our case by discreteness . By the two different types of numerical simulations used in our research, we have detected the internal mode bifurcating from the edge of the phonon band (of the continuum problem). The frequency of that mode agrees qualitatively very well with the one predicted in . Therefore, our research fully confirms the results of those papers and, in addition, shows that these perturbation induced internal modes behave very much like intrinsic ones. Our result, along with the finding of similar internal modes in other nonintegrable sG systems , strongly supports the generality of the phenomenon predicted in . Finally, the fact that we are able to detect this novel mode reinforces our previous conclusion of the non-existence of the internal quasi-mode, as indeed the sG kink responds to the ac driving or to deformations as it should if it had an internal mode, but for the discreteness induced one only. We thus believe that the present work definitely rules out the possibility of internal modes or “quasimodes” of sG kinks close to and above the lower phonon band edge.

## V Acknowledgments

We thank Yuri Kivshar and Yuri Gaididei for discussions. Work at GISC (Leganés) has been supported by DGESIC (Spain) grant PB96-0119. Travel between Bayreuth and Madrid has been supported by “Acciones Integradas Hispano-Alemanas”, a joint program of DAAD (Az. 314-AI) and DGESIC.


text:
# 1 Description of the model

The $`\rho `$ Meson in a Nuclear Medium

D. Cabrera, E. Oset, M. J. Vicente-Vacas

Departamento de Física Teórica e IFIC, Centro Mixto Universidad de Valencia-CSIC,
46100 Burjassot (Valencia), Spain

## Abstract

In this work, propagation properties of the $`\rho `$ meson in symmetric nuclear matter are studied. We make use of a coupled channel unitary approach to meson-meson scattering, calculated from the lowest order Chiral Perturbation Theory ($`\chi PT`$) lagrangian including explicit resonance fields. Low energy chiral constraints are considered by matching our expressions to those of one loop $`\chi PT`$. To account for the medium corrections, the $`\rho `$ couples to $`\pi \pi `$ and $`K\overline{K}`$ pairs which are properly renormalized in the nuclear medium.

PACS: 14.40.Aq; 14.40.Cs; 13.75.Lb

## 1 Description of the model

We study the $`\rho `$ propagation properties by obtaining the $`\pi \pi `$ and $`K\overline{K}`$ scattering amplitudes in the $`I=1`$ channel. Our states in the isospin basis are

$`|\pi \pi >`$ $`=`$ $`{\displaystyle \frac{1}{2}}|\pi ^+\pi ^{}\pi ^{}\pi ^+>`$
$`|K\overline{K}>`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}|K^+K^{}K^0\overline{K}^0>.`$ (1)
Tree level amplitudes are obtained from the lowest order $`\chi PT`$ and explicit resonance field lagrangians of refs. , . We collect this amplitudes in a $`2\times 2`$ $`K`$ matrix whose elements are

$`K_{11}(s)`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{p_1^2}{f^2}}\left[1+{\displaystyle \frac{2G_V^2}{f^2}}{\displaystyle \frac{s}{M_\rho ^2s}}\right]`$
$`K_{12}(s)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{3}}{\displaystyle \frac{p_1p_2}{f^2}}\left[1+{\displaystyle \frac{2G_V^2}{f^2}}{\displaystyle \frac{s}{M_\rho ^2s}}\right]`$
$`K_{21}(s)`$ $`=`$ $`K_{12}(s)`$
$`K_{22}(s)`$ $`=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{p_2^2}{f^2}}\left[1+{\displaystyle \frac{2G_V^2}{f^2}}{\displaystyle \frac{s}{M_\rho ^2s}}\right]`$ (2)
with the labels 1 for $`\pi \pi `$ and 2 for $`K\overline{K}`$ states. In equation (S1.Ex2) $`G_V`$ is the strength of the pseudoscalar-vector resonance vertex, $`f`$ the pion decay constant in the chiral limit, $`s`$ the squared invariant mass, $`M_\rho `$ the bare mass of the $`\rho `$ meson and $`p_i=\sqrt{s/4m_i^2}`$.

The final expression of the $`T`$ matrix is obtained following the N/D method, which was adapted to the context of chiral theory in ref. . It reads

$$T(s)=\left[I+K(s)g(s)\right]^1K(s),$$
(3)
where $`g(s)`$ is a diagonal matrix given by the loop of two mesons. In dimensional regularization it reads

$$g_i(s)=\frac{1}{16\pi ^2}\left[2+d_i+\sigma _i(s)\mathrm{log}\frac{\sigma _i(s)+1}{\sigma _i(s)1}\right],$$
(4)
where the subindex $`i`$ refers to the corresponding two meson state and $`\sigma _i(s)=\sqrt{14m_i^2/s}`$ with $`m_i`$ the mass of the particles in the state $`i`$. At this stage (vacuum case), the model has proved to be successful in describing $`\pi \pi `$ P-wave phase shifts and $`\pi `$, $`K`$ electromagnetic vector form factors  up to $`\sqrt{s}1.2`$ GeV. The $`d_i`$ constants contain the information of low energy chiral constraints. They are obtained by matching the expressions of the form factors calculated in this approach with those of one-loop $`\chi PT`$.

In our calculation, in which $`g(s)`$ is modified in the nuclear medium, we use cut-off regularization. The $`g_i(s)`$ function with a cut-off in the three-momentum of the particles in the loop can be found in Appendix A of ref. . By comparing the expressions in both schemes we can get the equivalent $`q_i^{max}`$ in order to keep the information of the $`d_i`$ constants.

Medium corrections are incorporated in the selfenergies of the mesons in the loops. All the graphs in fig. 1 are considered. In addition to the single $`ph`$ bubble we have to account for the medium modifications of the $`\rho MM^{^{}}`$ vertex via the $`\rho MN`$ contact term requested by the gauge invariance of the theory.

The pion selfenergy is written as usual in terms of the Lindhard functions. Both $`Nh`$ and $`\mathrm{\Delta }h`$ excitations are included. Short range correlations are also accounted for with the Landau-Migdal parameter $`g^{}`$, set to $`0.7`$. The final expression is

$$\mathrm{\Pi }_\pi (q,\rho )=\stackrel{}{q}^2\frac{(\frac{D+F}{2f})^2U(q,\rho )}{1(\frac{D+F}{2f})^2g^{}U(q,\rho )}$$
(5)
where $`U=U_N+U_\mathrm{\Delta }`$, the ordinary Lindhard function for $`ph`$, $`\mathrm{\Delta }h`$ excitations .

The $`\overline{K}`$ selfenergy has both S-wave and P-wave contributions. The S-wave piece is obtained from a self-consistent calculation with coupled channels ($`\overline{K}N`$, $`\pi \mathrm{\Sigma }`$, $`\pi \mathrm{\Lambda }`$, $`\eta \mathrm{\Sigma }`$, $`\eta \mathrm{\Lambda }`$, $`K\mathrm{\Xi }`$) in which both meson and baryon selfenergies in the medium have been considered. The P-wave piece includes $`\mathrm{\Lambda }h`$, $`\mathrm{\Sigma }h`$ and $`\mathrm{\Sigma }^{}(1385)h`$ excitations. The whole $`\overline{K}`$ selfenergy is borrowed from ref. .

At low energies the $`K`$ system interacts with nucleons only by S-wave elastic scattering. We use the expression for the selfenergy from ref. ,

$$\mathrm{\Pi }_K(\rho )0.13m_K^2\frac{\rho }{\rho _0}(MeV^2)$$
(6)
## 2 Results and discussion

We have plotted in fig. 2 the real and imaginary parts of $`T_{22}`$ for several densities. As can be seen, the resonance broadens significantly as density is increased, its width being around $`350`$ MeV at $`\rho =\rho _0`$. The zero of the real part experiences a downward shift which amounts to $`100150`$ MeV at normal density.

Another interesting result comes from the comparison between coupled and decoupled cases. This is done by setting $`K_{12}(s)`$ to zero, what automatically makes the $`T`$ matrix diagonal. We have found very small differences in $`T_{22}`$ when calculating in these two cases. This tells us that in our model the $`K\overline{K}`$ system has almost no influence on the $`\pi \pi \pi \pi `$ channel even at normal nuclear density.


text:
# On the integration accuracy in molecular density functional theory calculations using Gaussian basis sets

## I Introduction and theoretical background

Density Functional Theory (DFT) has in recent years become both a powerful and a very commonly used tool for studying molecular electronic structure. In addition to usually reliable energetics, the availability of analytical first and second derivatives for the DFT approaches allows for easy calculation of geometric structure and vibrational frequencies.

The rate-determining step in molecular DFT calculations using Gaussian basis sets is usually the evaluation of various integrals over the exchange-correlation functional and its geometric derivatives. In practice the three-dimensional numerical integration is carried out as a weighted sum of three-dimensional integrations in polar coordinates centered on atoms $`A`$:

$$F(𝐫)𝑑\tau =\underset{A}{\overset{\mathrm{atoms}}{}}\underset{g}{\overset{\mathrm{grid}}{}}w_{gA}p_A(𝐫_𝐠)F(𝐫_{gA})$$
(1)
in which the $`𝐫_g`$ and $`w_g`$ are numerical quadrature abscissae and weights, respectively, and the partitioning function $`p_A`$ satisfies $`_Ap_A(𝐫)=1`$ for all values of $`𝐫`$. In other words, $`p_A`$ acts as a ‘fader function’ between grids centered on different atoms.

Quite separate from issues involving the basis set and mathematical form of the exchange-correlation functional, which determine the accuracy of the calculation, the choice of atomic integration grid and partitioning can have a profound impact on its numerical precision (together with trivial parameters like convergence criteria). Especially since an increasing proportion of users of quantum chemical software in general — and DFT methods in particular — are experimental chemists without prior background in electronic structure theory, the tendency exists to simply rely on the built-in defaults of various codes.

In the present paper, we shall report a number of problems which may arise in this fashion. Because of its great popularity and widespread use, our discussion will focus on the Gaussian 98 program system. We wish to emphasize, however, that the issues raised in the present work could arise with other DFT codes. We note that some codes try to minimize possible numerical problems by the use of large default integration grids and/or by variation of the grid mesh with the atom type. While the use of very large grids certainly avoids numerical problems, it can lead to prohibitively expensive calculations for larger systems. Therefore, it is of interest to find the smallest grid that will yield precise results for molecular properties.

Before proceeding to discuss our results, we shall first briefly review some relevant mathematical aspects.

## II Mathematical background

The atomic integration grid is built up naturally as the direct product of radial and angular grids:

$$\underset{g}{}w_gF(\mathrm{r}_g)=\underset{a}{}\underset{b}{}w_aw_bF(r_a,\theta _b,\varphi _b)$$
(2)
The angular grid itself can be either a direct product of appropriate one-dimensional Gauss quadrature grids (e.g., Gauss-Legendre or Gauss-Chebyshev), or (more efficiently) can be a two-dimensional Lebedev grid which exactly integrates spherical harmonics up to a certain order. (In order to integrate spherical harmonics up to order $`L`$ exactly, the direct product of grids with $`(L+1)/2`$ abscissae in $`\theta `$ and $`L+1`$ abscissae in $`\varphi `$ requires $`(L+1)^2/2`$ angular quadrature points, compared to approximately $`(L+1)^2/3`$ for the corresponding Lebedev grid.)

The radial grids discussed in the present paper are all of the Murray-Handy-Laming (MHL) variety (following ideas first proposed by Boys and Handy), in which the infinite interval $`r=[0,\mathrm{}]`$ is remapped onto the finite interval $`q=[0,1]`$ by means of the change of variable $`r=\alpha (q/(1q))^m`$; the integration in $`q`$ is then approximated by Euler-Maclaurin summation. (The latter is facilitated by the fact that $`\left(\frac{d^nr}{dq^n}\right)`$ vanishes at $`q=1`$ for all $`n`$, and at $`q=0`$ up to $`n=3m1`$ inclusive.) Empirically, $`m=`$2 is found to yield the best results, and $`\alpha `$ should be dependent on the atomic radius for best results. The following working equations arise for abscissae and weights:

$`r_i`$ $`=`$ $`R(i/(N+1i))^2`$ (3)
$`w_i`$ $`=`$ $`2R^3(N+1)i^5/(N+1i)^7`$ (4)
in which $`N`$ is the number of radial points and $`R`$ is taken empirically as one-half the Bragg radius of the atom on which the grid is centered (except for hydrogen, where $`R`$ is set equal to the Bohr radius).

For molecular calculations, the partitioning function $`p(𝐫)`$ introduced in Eq.(1) needs to be defined. One common choice is due to Becke. Consider the confocal elliptical coordinate $`\mu _{AB}=r_Ar_B/R_{AB}`$ (in which $`r_A`$ and $`r_B`$ are distances towards nuclei $`A`$ and $`B`$, and $`R_{AB}`$ represents the internuclear distance), and define $`p(\mu )=\frac{1}{2}[1g(\mu )]`$, where $`g(\mu )`$ is the ‘cell function’. The extreme choice $`g(\mu _{AB})`$=$`1`$ if $`\mu _{AB}0`$; =+1 if $`\mu _{AB}>0`$ (i.e. a Heaviside-type step function) would actually represent a scheme in which space is simply partitioned into Voronoi polyhedra: this will obviously be a numerically unstable scheme, since molecules by definition have appreciable electron density $`\rho `$ (of which the integrand is a functional) near bond midpoints $`\mu _{AB}=0`$. Becke proposed to use a smoother approximation $`g(\mu )=f(f(f(\mu )))`$, where $`f(\mu )=\frac{3}{2}\mu \frac{1}{2}\mu ^3`$. Stratmann, Scuseria, and Frisch (SSF) proposed to use a piecewise function instead, namely $`g(\mu _{AB})`$=$`1`$ if $`\mu _{AB}a`$; =+1 if $`\mu _{AB}a`$; $`=\frac{1}{16}[35(\mu _{AB}/a)35(\mu _{AB}/a)^3+21(\mu _{AB}/a)^55(\mu _{AB}/a)^7]`$ otherwise, where the scale factor $`a=0.64`$ was determined empirically.***As pointed out by a referee, this latter expression is equivalent to eq. (24), with $`m`$=3, of Ref., although a piecewise $`g(\mu )`$ was not proposed there. The SSF partitioning is reported to be numerically somewhat more stable, and in addition lends itself well to linear-scaling implementations.

## III Computational methods

All calculations were carried out using Gaussian 98 Rev. A7 running on various workstation computers at NASA Ames Research Center and the Weizmann Institute of Science.

A variety of Gaussian basis sets was used, including the Pople group 6-31G\*, 6-31+G\*, and 6-31++G\** basis sets, the Dunning cc-pVDZ and cc-pVTZ basis sets (correlation consistent polarized valence double and triple zeta, respectively), and the Hay-Wadt LANL2DZ (Los Alamos National Laboratory 2 Double-Zeta) basis set-relativistic ECP combination.

In virtually all cases, the very popular B3LYP (Becke 3-parameter-Lee-Yang-Parr) hybrid exchange-correlation functional was used, except for some test calculations using the earlier BLYP (Becke-Lee-Yang-Parr) or BP86 (Becke-Perdew 1986) functionals.

We consider three types of grids:

* ($`nn`$,$`mmm`$), meaning the direct product of an $`nn`$-point Euler-Maclaurin radial grid with a two-dimensional Lebedev grid with $`mmm`$ angular points;
* ($`nn`$,$`mm`$,2$`mm`$), meaning the direct product of an Euler-Maclaurin grid with $`nn`$ points in $`r`$ and of Gauss-Legendre quadratures with $`mm`$ and $`2mm`$ points in $`\theta `$ and $`\varphi `$, respectively;
* ($`nn`$,$`mmm`$)p, i.e. pruned ($`nn`$,$`mmm`$) grids. In pruned grids, the number of angular points is varied with the radial coordinate, with the full number of $`mmm`$ only being used at $`r`$ values relevant for conventional chemical bonding. In Gaussian 98, the following pruned grids are implemented: (35,110)p (invoked with the keyword “Grid=Coarse”), (50,194)p (the SG-1 grid invoked with the keyword “Grid=SG1”), (75,302)p (invoked with the keyword “Grid=Fine”), and (99,590)p (invoked as “Grid=UltraFine”). Pruned grids are only implemented for elements up to and including krypton: for heavier elements, the corresponding unpruned grid is substituted. (We note that in the implementation in Gaussian 98, the cutoff radii between segments with different angular grid meshes are dependent on the row in the Periodic Table: this in effect does result in some grid quality adjustment with the atomic number.)

In frequency calculations using analytical second derivatives, Gaussian 98 permits the use of different grids for the evaluation of the energy and gradient on the one hand, and for the solution of the CPKS (coupled perturbed Kohn-Sham) equations on the other hand. The defaults of the program correspond to “Int(Grid=Fine) CPHF(Grid=Coarse)”; specifying “Int(Grid=UltraFine)” implies “CPHF(Grid=SG1)” by default.In Gaussian 98, the keyword CPHF controls options for both the Coupled Perturbed Hartree-Fock procedure in wavefunction based ab initio calculations and for CPKS in density functional calculations.

Unless noted otherwise, tightened convergence criteria were used in geometry optimizations.

## IV Results and discussion

### A CCH radical

The acetyl radical is linear with a $`{}_{}{}^{2}\mathrm{\Sigma }_{}^{+}`$ ground state. At the B3LYP/cc-pVTZ level with the default grid combination, however, one surprisingly finds (Table I) not only an imaginary bending frequency, but that the two “degenerate” components differ by 64 cm-1. Upon repeating the calculation with an extremely large (140,48,96) grid, a properly degenerate bending frequency of 312.5 cm-1 (intensity 3.84 km/mol) was found, irrespective of whether Becke or SSF partitioning is used. Taking the SSF (140,48,96) result as the reference, we find the default grid to be in error by only 0.3 and 0.1 cm-1 for the CC and CH stretches, respectively, while the two components of the bending frequency are off by a whopping $`421`$ and $`492`$ cm-1, respectively. (N.b.: we have treated the imaginary frequencies as negative numbers when we compute the error.) Likewise, while errors of about 5 and 1 % are seen for the CC and CH intensities, respectively, the bend intensity is too high by an order of magnitude. Using Becke instead of SSF weights does not materially affect these observations.

Specifying “Grid=Ultrafine” (i.e. (99,590)p for integration, and (50,194)p for CPKS) leads to a dramatic improvement in the bending frequency, which is now properly degenerate and only 7 cm-1 too low. (The intensities for the other two modes now agree with the reference values to two decimal places.) The error in the corresponding infrared intensity drops by two orders of magnitude. Specifying the (99,590)p grid for both integration and CPKS does not affect the stretching frequencies, but reduces the error on the bending frequency and intensity to 0.1 cm-1 and 0.01 km/mol, respectively, using either Becke or SSF partitioning. This clearly suggests that sensitivity is greatest to the CPKS grid.

Using the (75,302)p grid for integration and varying the CPKS grid from the associated (35,110)p default, we find that the spurious imaginary, nondegenerate bending frequencies disappear upon using either the unpruned (35,110) or the larger pruned (50,194)p grid. However, errors of about –10 and –7 cm-1 in the bending frequency remain, and are not reduced to the 1 cm-1 range until a (75,302)p or even (75,302) CPKS grid is substituted. Upgrading the energy grid to (99,590)p does not affect these conclusions.

Surprisingly, using a (99,974) grid for both integration and CPKS yields slightly different frequencies for the two components of the $`\pi `$ bending mode.

Upon finding an imaginary frequency, one is tempted to displace the molecular geometry in the direction of that mode and reoptimize the geometry. Following this procedure, using the Gaussian 98 defaults (an energy grid of (75,302)p, a CPKS grid of (35,110)p, and the SSF weights), one obtains a HCC angle of 171.1 degrees and three real frequencies. However, one of the six trivial vibrations (i.e., overall translation and rotation, which should have exactly vanishing frequencies in an infinitely precise calculation using a rotationally invariant method) has a frequency of $`508i`$ cm-1. Thus lowering symmetry has just made the numerical problem less obvious.

Unlike the frequencies which show a large sensitivity to the choice of grid, the total energy of CCH is in error by less than 1 microhartree for all grids studied. Considering the small effect on the energy, it is not surprising that the CCH results are most sensitive to the CPKS grid.

### B Harmonic frequencies of TaCl2

The ground state for this symmetric linear molecule is $`X{}_{}{}^{4}\mathrm{\Sigma }_{g}^{}`$. Results at the B3LYP/LANL2DZ level as a function of the integration and CPKS grids are summarized in Table II. Our largest grid (140,48,96) SSF results will again be considered as the reference.

None of the grid combinations considered in Table II lead to nonequivalent bending vibration components, but the default grid combination does produce a spurious imaginary bending frequency. Unlike CCH, the bending infrared intensities agree well with the largest grid value, but the symmetric and antisymmetric stretching frequencies are in error by $`6`$ and $`7`$ cm-1, respectively. Use of the Becke partitioning leads not only to increased errors in the stretching frequencies, but to an error of +0.002 Å in the optimum bond distance, compared to $`0.0003`$ Å with the SSF partitioning. Quite different from the behavior for CCH is the fact that the total energy is now in error by 310 microhartree with the SSF partitioning, and 369 microhartree with the Becke partitioning.

In this particular case, varying the size of the CPKS grid, and particularly of its angular component, does not appear to remedy the problem at hand. Using an unpruned energy grid, however Note that this change will only affect the Cl atomic grid in this case, as grids for elements heavier than Kr are always unpruned in Gaussian 98., or using a grid with more radial points, dramatically reduces the error in the computed total energy, bond distance, and frequencies, particularly using the SSF partitioning. In particular, the (99,590)p/(50,194)p integration/CPKS grid combination invoked by “Grid=UltraFine” incurs (with SSF partitioning) no error greater than 1 cm-1 on the frequencies, and reduces the errors in bond distance and total energy to 0.0001 Å and 1.5 microhartree, respectively. While it is possible that using the Becke atomic size adjustment procedure (Appendix of Ref.), in conjuction with the SSF weights, would reduce the sensitivity to the radial grid, we note that using the ‘UltraFine’ grid combination solves both the CCH angular grid problem and the TaCl2 radial grid problem, and therefore would seem to be a good choice to invesigate possible grid problems.

For TaCl2, performance of the SSF partitioning seems markedly superior to that of the Becke partitioning. As expected, differences between the two partitioning schemes are reduced as the grid is improved, and they yield essentially identical results for the (140,48,96) grid.

Lowering the symmetry to $`C_{2v}`$ and using the Gaussian 98 defaults results in a slightly bent structure with a ClTaCl angle of 171.1 degrees. As in the case of CCH, the only indication of the numerical problems is in the deviation from zero of the translational and rotational frequencies, but for TaCl2 the value is $`33i`$ cm-1, making the error for TaCl2 less obvious than in the case of CCH (see above).

If we now consider relative CPU times (default grids=1.00) as a function of grid size, we see that the (99,590)p/(50,194)p combination will approximately double CPU time. Using a (99,590)p grid for both integration and CPKS is found to be a factor of nine more expensive than the default, a ratio which goes up to a factor of 19 for the unpruned (99,974) grid, of 36 for the (96,32,64) grid often cited in benchmarks, and of 109 for the largest (140,48,96) grid considered here. The (99,590)p/(75,302)p combination found to be required for 1 cm-1 precision in the bending frequency of CCH would be four times more expensive than the default — rather more time-consuming than desirable, but still an order of magnitude less expensive than the (96,32,64) grid.

### C Harmonic frequencies of Ge2H5 and Ge2H6

Ge2H6 has $`D_{3d}`$ symmetry and a structure analogous to ethane. Ge2H5 is best viewed as removing one H from Ge2H6, which results in $`C_s`$ symmetry; there is only a small change in the structural parameters relative to the parent Ge2H6 molecule.

In Table III we summarize our results for Ge2H5. In addition to the effect of partitioning and grid size, we also consider the influence of the basis set and of the exchange-correlation functional on the three lowest harmonic frequencies of Ge2H5. We first consider using the Becke partitioning. Using the default grid in Gaussian, i.e. a pruned (75,302)p grid for energy and gradients and (35,110)p for CPKS, results in an imaginary frequency. Expanding the radial grid from 75 to 96 points and depruning the grid, i.e. the (96,302) grid, for the energy, gradient, and CPKS steps does not remove the imaginary frequency. Expanding the angular grid removes the numerical problem; with 590 angular points, the lowest harmonic frequency is essentially converged to a value of about 100 cm-1. Expanding the radial grid further has only an effect of a few cm-1.

The results using SSF partitioning are very different from those obtained using Becke partitioning; using SSF, even the default “Grid=Fine” is precise to better than 10 cm-1. Improving the grid makes much smaller changes for SSF and the “UltraFine” grid in Gaussian 98 (i.e., (99,590)p for energy and gradients, (50,194)p for CPKS) appears to be essentially converged. Note that for the Becke partitioning the “UltraFine” grid also yields essentially converged results.

Adding diffuse functions to the Ge atoms or diffuse and polarization functions to the H atoms does not eliminate the imaginary frequency for the default grid and the Becke partitioning. That is, the problem is not unique to the 6-31G\* basis set.

Nor is it unique to the use of a hybrid functional: substituting the BP86 functional for its B3LYP counterpart does not alter the trends with respect to partitioning and grid, the most notable point being the imaginary frequency for the default grid with Becke partitioning.

The Ge2H5 calculations are repeated using the LANL2DZ basis set/ECP combination, and the results are also given in Table III. Unlike the all-electron calculations, there is essentially no difference between Becke and SSF partitioning, and improving the grid has only a very small effect. This is also true for the four- and 14-valence electron Ge ECPs of Hurley et al.. Thus it appears that the difference observed in the all-electron calculations arises from the inner-shell orbitals.

For the default grids and the Becke partitioning, displacing the molecular geometry in the direction of the imaginary mode results in essentially the same geometry and frequencies as found in the $`C_s`$ treatment. That is, lowering the symmetry does not eliminate the imaginary frequency, so the numerical problems are much more obvious than for CCH and TaCl2.

Substituting Si for Ge leads to a similar structure, but as shown in Table IV, neither the choice of grid points nor the partitioning scheme make much of a difference for the Si2H5 system. The lack of sensitivity of the results to grid or partitioning is similar to the results reported by SSF, and is typical of many other tests that we have performed for systems containing first and second row atoms. The problem of the imaginary frequency does not exist for Ga2Cl5: the Becke and SSF partitioning agree to within 4 cm-1 for the 6-31+G\* basis set and the default grids. Given that Si2H5, Ge2Cl5, and Ge2H5 have similar shapes, we conclude that the shape of Ge2H5 does not lead to the numerical problems that result in the imaginary frequency.

In Table V we summarize the results for Ge2H6, which is a closed-shell, stable molecule. It shows the same variation in the frequencies with choice of partitioning and grid points as Ge2H5, and we conclude that the difference between results using Becke and SSF partitioning is not unique to Ge2H5. It is encouraging that the default SSF partitioning is less sensitive to the choice of grid, but it is also sobering to realize that it is possible to have an error of about 70 cm-1 with no obvious indication of an error.

### D Relative energies of cis- and trans-HIr\[PR3\]4Cl+

The peculiar cis-trans equilibrium of the HIr\[PR3\]4Cl+ (R=H, CH3) transition metal complexes was the subject of a very recent joint theoretical-experimental investigation. At the highest level of theory for R=CH3, the complexes are found to be essentially isoenergetic, with the equilibrium between the two being essentially wholly driven by the higher entropy of the cis form: experimentally, a 6:1 cis:trans equilibrium is found at room temperature.

Given the very small energy differences being considered here, as well as the rather flat potential surface of particularly the cis form, grid convergence (not discussed in Ref.) was obviously a cause for concern here. Upon first optimizing the cis-R=CH3 case using the default grids and computing vibrational frequencies, we found what appeared to be a local minimum. Upon recomputing the harmonic frequencies with the (99,590)p/(50,194)p integral/CPKS grid combination, we found this structure to be a saddle point of order 3, and only after thirteen geometry optimization cycles using analytical second derivatives (and two months of CPU time on an SGI Origin 2000) was the true minimum structure found. A superimposition of the two structures (created using the MOLDEN molecular viewing program) is displayed in Figure 1.

While the R=CH3 case is simply too large to do an exhaustive grid convergence study, we have summarized some data for the R=H case in Table VI. In particular, we have considered the absolute energies for the cis and trans isomers, and their relative energy, as a function of grid size. In this case, results seems to be quite dependent both on the radial and on the angular grid mesh. Using the (35,110)p grid employed by default in the CPKS step, the relative energy gets predicted with the wrong sign, and neither the use of an unpruned grid nor switching to a product grid of the same order remedy the problem — if anything, results get worse, suggesting that the (35,110)p results in fact benefit from a fortunate error compensation. Total energies are off by as much as 8 millihartree for the trans isomer. While a scatter in the energies does not necessarily imply a corresponding scatter in the CPKS solution, it does not seem entiely implausible that a grid that predicts isomer energy differences on a surface with the wrong sign might yield a Hessian with incorrect curvature if used in the CPKS step.

The smallest grid that at least predicts the correct sign for the cis-trans difference is (50,194)p (about 25% too high in absolute value): using an unpruned grid instead§§§ Again, since the Ir grid is always unpruned, the improvement in grid quality comes from additional grid points centered on the other atoms.  cuts the error in half, as does using the (75,302)p grid. Total energies are however still in error by amounts of about 60 microhartree: fortunately the errors largely cancel. The (99,590)p (“UltraFine”) grid agrees to two decimal places with the (140,48,96) reference results; even a grid as large as (99,974) cannot yet achieve 1 microhartree precision in the total energies, which is barely reached for a (99,27,54) grid, and more comfortably for a (99,32,64) grid.

## V Conclusions

Computed molecular properties using density functional theory in general, and vibrational frequencies in particular, tend to be quite sensitive to the integration grids and partitioning being used in the calculation. Our results corroborate the assertion of SSF that their proposed partitioning scheme is numerically more stable than the Becke partitioning.

While the use of coarse grids in the solution of the CPKS (coupled perturbed Kohn-Sham) equations certainly leads to significant speedups, it can cause spurious negative eigenvalues in the Hessian matrix, or even artifactual loss of vibrational mode degeneracy, in systems as small and uncomplicated as CCH radical. While this problem is alleviated by the use of unpruned grids, this offers no benefit (in terms of CPU time) over the use of a larger pruned grid. For systems involving very heavy elements (e.g. third-row transition metals), grids which are sufficiently fine-meshed for typical first-and second-row systems are no longer adequate. In situations where small pruned grids are inadequate, the use of the corresponding unpruned grid appears to be less efficient than the use of a larger pruned grid.

For the cases considered here, the most efficient grid combination that appears to be ‘immune’ to qualitatively incorrect results is (99,590)p for energy and gradients, (50,194)p for CPKS. (This is equivalent to the “Grid=UltraFine” option in Gaussian 98.) For greater quantitative precision, the use of a finer-meshed (75,302)p grid for the CPKS is recommended. While a (99,974), (96,32,64) or even finer grid in all steps of the calculation will certainly avoid this type of problems, this comes at an unacceptably high premium in terms of CPU time.

As quantum chemical methods mature and quantum chemical software becomes more user-friendly, an increasing proportion of their users are scientists from other fields than quantum chemistry, who are inevitably exposed to the temptation to treat such program systems as ‘black boxes’. Observations such as those made in the present paper illustrate that even at the present state of technology, ‘black box’ performance should not be taken for granted, and that acquiring some insight in both the quantum mechanical and numerical methodology is still an essential prerequisite for their effective and reliable scientific application.

###### Acknowledgements.

JM is the incumbent of the Helen and Milton A. Kimmelman Career Development Chair. Research at the Weizmann Institute was supported by the Minerva Foundation, Munich, Germany, and by the Tashtiyot program of the Ministry of Science (Israel). JM would like to thank Dr. S. Parthiban for bringing the problem with CCH to his attention, and Dr. Andreas Sundermann for critical reading of the draft manuscript. AR was supported by NASA Contract No. NAS2-14031.

## Note added

After acceptance of the present manuscript, and in response to a preprint thereof, we received a personal communication from Gary W. Trucks and Michael J. Frisch of Gaussian, Inc., informing us of a bug in Gaussian 98 which affects the integral (derivative) accuracy for the specific case of the (35,110)p grid. A fix will be incorporated in the next minor release after Rev.A9. The main effect on the data presented in this paper is that the errors in the first two entries in Table 1 become a good deal less dramatic: in particular, the degeneracy of the bend is restored and the errors in $`\omega _3`$ and $`I_3`$ are reduced to -19 cm-1 and 1.2 km/mol, respectively. In addition, we reoptimized the complex in Figure 1 using a patched version of the code with the default grid combination (and using analytical second derivatives as before). The resulting geometry is in much closer agreement with the Grid=UltraFine structure. (The cis-trans difference with the (35,110)p grid in Table VI is not affected materially.) The authors would like to thank Drs. Trucks and Frisch for helpful discussions.