text:
# Scattering and Trapping of Nonlinear Schrödinger Solitons in External Potentials

abstract
Soliton motion in some external potentials is studied using the nonlinear Schrödinger equation. Solitons are scattered by a potential wall. Solitons propagate almost freely or are trapped in a periodic potential. The critical kinetic energy for reflection and trapping is evaluated approximately with a variational method.
keywords: soliton, nonlinear Schrödinger equation, Lagrangian

The nonlinear Schrödinger equation has been intensively studied as a soliton equation, which is derived as an envelope equation by several perturbation methods . Solitons in optical fibers are typical ones described by the nonlinear Schrödinger equation . Recent successful observation of solitons in Bose-Einstein condensates (BECs) raises the theoretical consideration of this equation again . The equation is called the Gross-Pitaevskii equation in the research field of BECs. In the BECs, matter waves obey the nonlinear Schrödinger equation with external potentials. For example, external magnetic fields are applied to confine the BECs. Interference patterns of laser beams generate external periodic potentials. Solitons in external potentials such as harmonic potential, periodic potential and stepwise potential have been studied by several authors . But solitons in a potential wall or potential well were not studied in detail. In this letter, we will study the scattering and trapping of solitons by a potential wall (or well) and a periodic potential. In particular, we focus on soliton behaviors like classical mechanical particles, and evaluate the critical energies for the transmission and the trapping.

The nonlinear Schrödinger equation with a potential term is written as

$$i\frac{\varphi }{t}=\frac{1}{2}\frac{^2\varphi }{x^2}|\varphi |^2\varphi +U(x)\varphi ,$$
(1)
where $`\varphi (x,t)`$ is the wave function, $`U(x)`$ is an external potential. Without the nonlinear term, this equation is equivalent to the Schrödinger equation with mass $`m=1`$ and $`\mathrm{}=1`$. When $`U(x)=0`$, the nonlinear Schrödinger equation has a soliton solution

$$\varphi (x,t)=A\mathrm{sech}\{A(xvt)\}e^{ip(xvt)i\omega t},$$
(2)
where $`v=p`$ and $`\omega =A^2/2+v^2/2`$.

Firstly, we investigate the nonlinear Schrödinger equation in a simple rectangular potential, that is,

$$U(x)=U_0,\mathrm{for}x_0<x<x_0+d,U(x)=0,\mathrm{for}x<x_0,x>x_0+d.$$
(3)
If $`U_0`$ is positive (negative), $`U(x)`$ represents a potential wall (well). We consider first the case of $`U_0>0`$. If the nonlinear term is absent in eq. (1), this is a typical potential problem of the quantum mechanics. The transmission coefficient $`T`$ of matter waves is expressed as

$$T=\left[1+\frac{U_0^2\mathrm{sin}^2(\alpha d)}{4E(EU_0)}\right]^1,\mathrm{for}E>U_0,T=\left[1+\frac{U_0^2\mathrm{sinh}^2(\alpha ^{}d)}{4E(EU_0)}\right]^1,\mathrm{for}E<U_0,$$
(4)
where $`E`$ denotes the kinetic energy $`E=v^2/2`$, and $`\alpha `$ and $`\alpha ^{}`$ are defined as $`\alpha =\sqrt{2(EU_0)}`$ (for $`E>U_0`$) and $`\alpha ^{}=\sqrt{2(U_0E)}`$ (for $`E<U_0`$). We study numerically the scattering of a soliton by the potential wall. Figure 1(a) displays the time evolution of a soliton passing through a potential wall with $`U_0=0.5`$ and $`d=1`$. The system size is $`L=100`$ and the potential wall locates at $`x=x_0=L/2=50`$. The initial velocity of the soliton is $`v=0.75`$ and the amplitude is $`A=1`$. The soliton is split into two pulses by the potential wall. The transmission coefficient is 0.527, which is calculated as $`T=_{L/2}^L|\varphi |^2𝑑x/_0^L|\varphi |^2𝑑x`$ after the scattering. Figure 1(b) displays the transmission coefficient as a function of the initial velocity $`v`$. The dashed curve is the transmission coefficient eq. (4). The transmission coefficient of a soliton is close to the coefficient eq. (4) for $`v>0.9`$, which might imply that the effect of the nonlinear term becomes relatively weaker as the kinetic energy is increased. The transmission coefficient decreases more sharply than the quantum mechanical particle as $`v`$ is decreased.

The transmission coefficients as shown in Fig. 1(b) change more and more sharply as $`U_0`$ is decreased. As $`U_0`$ becomes smaller, such a splitting or additional radiations are hardly observed and the soliton keeps its form and behaves like a particle. The coefficient jumps almost from 0 to 1 at a certain critical velocity for $`U_0<0.05`$. To understand the transition between the wavy behaviors including the splitting and the particle-like behavior, we have performed numerical simulations for another initial condition $`\varphi (x,0)=A\mathrm{sech}\{A(xL/2)\}`$, that is, the soliton is initially set at $`L/2`$ and the initial velocity is 0. The potential wall is located in $`[L/2d/2,L/2+d/2]`$. The mirror-symmetry around $`L/2`$ is maintained during the time evolution. For large $`U_0`$, the soliton is split into two pulses as shown in Fig. 2(a) for $`U_0=0.67,A=2,d=1`$ and $`L=80`$ owing to the repulsive potential. However, the soliton cannot be split for small $`U_0`$ owing to the attractive interaction by the nonlinear term. The soliton keeps its one pulse structure as shown in Fig. 2(b) for $`U_0=0.66,A=2,d=1`$ and $`L=80`$. The critical value of $`U_0`$ depends on the initial amplitude $`A`$ as shown in Fig. 2(c). Equation (1) has at least two invariants, that is, the total mass $`N=_0^L|\varphi |^2𝑑x`$ and the total energy $`E=_0^L\{1/2|\varphi _x|^21/2|\varphi |^4+U(x)|\varphi |^2\}𝑑x`$. If the soliton is split into two counterpropagating solitons: $`\varphi (x,t)=A^{}\mathrm{sech}\{A^{}(xvtL/2)\}e^{ip(xvtL/2)i\omega t}+A^{}\mathrm{sech}\{A^{}(x+vtL/2)\}e^{ip(x+vtL/2)i\omega t}`$, the amplitude $`A^{}`$ and the momentum $`p`$ satisfy $`A^{}=A/2`$ and $`2p^2A^{}=1/4A^3+2U_0A\mathrm{tanh}(Ad/2)`$ owing to the conservation laws. If $`U_0<A^2/\{8\mathrm{tanh}(Ad/2)\}`$, there is no real solution for $`p`$, that is, the splitting cannot occur and the one-pulse structure is maintained. Above the line $`U_0=A^2/\{8\mathrm{tanh}(Ad/2)\}`$, the potential energy is used to split a soliton against the attractive interaction. A soliton behaves like a particle even for more general initial conditions, if $`U_0`$ is sufficiently smaller than $`A^2/\{8\mathrm{tanh}(Ad/2)\}`$.

We focus on the particle-like behavior below. We return to the numerical simulations of the first type. Figures 3(a) and (b) display the time evolutions of a soliton passing through a potential wall of $`A=1,U_0=0.04`$ and $`d=1`$. The system size is 50 and the potential wall locates at $`x_0=L/2=25`$. The soliton is almost completely reflected when the initial velocity is $`v=0.18`$, and almost completely transmitted for $`v=0.2`$. The critical velocity for the soliton to pass through the

potential wall is 0.1925, and the critical kinetic energy is 0.0185, which is smaller than the critical value $`E_c=U_0=0.04`$ for a particle obeying the classical mechanics. Figure 3(c) displays the critical kinetic energy as a function of the amplitude $`A`$ of the soliton. As the amplitude $`A`$ is increased, the critical energy approaches $`U_0`$. A soliton with smaller amplitude can pass easily through the potential wall. This may be interpreted as coherent tunneling.

The nonlinear Schrödinger equation is rewritten with the Lagrangian form: $`/t(\delta L/\delta \varphi _t)=\delta L/\delta \varphi `$, where $`L=_{\mathrm{}}^{\mathrm{}}𝑑x\{i/2(\varphi _t\varphi ^{}\varphi _t^{}\varphi )1/2|\varphi _x|^2+1/2|\varphi |^4U(x)|\varphi |^2\}`$. If a solution to the nonlinear Schrödinger equation is approximated by

$$\varphi (x,t)=A(t)\mathrm{sech}\{(x\xi (t))/a(t)\}e^{ip(t)(x\xi (t))i\omega t},$$
(5)
the time evolution of $`a(t),\xi (t)`$ and $`p(t)`$ can be derived from the effective Lagrangian $`L_{eff}=_{\mathrm{}}^{\mathrm{}}L𝑑x`$  using

$$L_{eff}=2A^2ap\xi _t+2/3A^4a1/3A^2/aA^2ap^22A^2aU_{eff},$$

where $`U_{eff}=U_0/2[\mathrm{tanh}\{(x_0+d\xi )/a\}\mathrm{tanh}\{(x_0\xi )/a\}].`$ Equations of motion for the variational parameters $`q_i=\xi ,p`$ and $`a`$ are given by $`d/dt(L_{eff}/\dot{q_i})=L_{eff}/q_i`$ under the condition $`N=_{\mathrm{}}^{\mathrm{}}|\varphi |^2𝑑x=2A^2a=const`$. The results are

$$\frac{d\xi }{dt}=p,\frac{dp}{dt}=\frac{U_{eff}}{\xi },\frac{da}{dt}=0$$
(6)
and $`A=1/a=N/2`$. The soliton obeys the Newton equation with an effective potential $`U_{eff}`$. The peak amplitude of the effective potential is $`U_0\mathrm{tanh}(Ad/2)`$. The critical kinetic energy is therefore expected to be $`E_c=U_0\mathrm{tanh}(Ad/2)`$, which was drawn with the dashed curve in Fig. 3(c). Good agreement of numerical results with this approximation is seen in Fig. 3(c).

Secondary, we consider the case of negative $`U_0`$. The potential acts as a potential well. We study the scattering of a soliton by the potential well. Complicated scattering including splitting is observed when $`|U_0|`$ is large, as is the case of positive $`U_0`$. We focus on the particle-like behavior, which is observed when $`|U_0|`$ is sufficiently small. Figure 4 displays the time evolutions of solitons for a potential well with $`U_0=0.04`$ and width $`d=1`$. The initial amplitude is $`A=1`$ and the initial velocities are $`v=0.035`$ and 0.015. For relatively large $`v`$, the soliton passes over the potential well like a classical mechanical particle obeying eq. (6). The velocity increases when the particle passes over the potential well, but its initial value is recovered after the passage. However, the soliton is trapped by the potential well for sufficiently small $`v`$ as shown in Fig. 4(b). The critical value of $`v`$ for the trapping is approximately $`v_c0.0225`$. In Fig. 4(a), a soliton passes through the potential well, but the velocity is decreased from 0.035 to 0.026 after the passage of the potential well. Some dissipative processes may occur at the potential well. Some small radiations may be generated and a part of the kinetic energy is changed into radiation energy. Another possibility is that some internal freedoms of a soliton such as breathing motion may be excited and the kinetic energy decreases. To investigate such motion, we assume a solution of the form

$$\varphi (x,t)=A(t)\mathrm{sech}\{(x\xi (t))/a(t)\}e^{ip(t)(x\xi (t))+i\sigma (t)\mathrm{log}\mathrm{cosh}\{(x\xi (t))/a(t)\}i\omega t},$$
(7)
where $`A,a,p,\xi `$ and $`\sigma `$ change in time, but $`N=2A^2a`$ is conserved in time. The term $`i\sigma (t)\mathrm{log}\mathrm{cosh}(x\xi (t))/a(t)`$ is a kind of chirp term and this form of chirp term appears in a solution to the complex Ginzburg-Landau equation where some dissipation and filtering effects are involved. Owing to the inclusion of this term, the width of the soliton changes in time and breathing motion becomes possible. The choice of this form of chirp term is not a unique one. We have studied another form of chirp term $`i\sigma (x\xi (t))^2`$, but better agreement with direct numerical simulation was seen for the form eq. (7). The substitution of eq. (7) to the Lagrangian and its spatial integration yield the effective Lagrangian

$$L_{eff}=N(p\xi _t\sigma _t(1\mathrm{log}2)+\sigma a_t/(2a)1/(6a^2)p^2/2\sigma ^2/(6a^2)+N/(6a)U_{eff}),$$

where $`U_{eff}=U_0/2[\mathrm{tanh}\{(x_0+d\xi )/a\}\mathrm{tanh}\{(x_0\xi )/a\}].`$ The equations of motions for $`p(t),\xi (t),a(t)`$ and $`\sigma (t)`$ are

$`{\displaystyle \frac{d\xi }{dt}}`$ $`=`$ $`p,{\displaystyle \frac{dp}{dt}}={\displaystyle \frac{U_{eff}}{\xi }},`$
$`{\displaystyle \frac{da}{dt}}`$ $`=`$ $`{\displaystyle \frac{2\sigma }{3a}},{\displaystyle \frac{d\sigma }{dt}}=2a\left({\displaystyle \frac{1+\sigma ^2}{3a^3}}{\displaystyle \frac{N}{6a^2}}{\displaystyle \frac{U_{eff}}{a}}\right).`$ (8)
The energy $`p^2/2+U_{eff}`$ is not conserved in this equation, and the width $`a`$ changes in time. If the energy decays, the soliton may be trapped in the potential well.

Figure 5 displays the time evolutions of $`\xi (t)`$ by eq. (8) for $`v=0.035`$ and $`v=0.015`$ with the dashed curves, and they are compared with the time evolutions of the peak position of the solitons shown in Figs. 4(a) and (b). The velocity of a particle changes from the initial value of 0.035 to 0.029 in the time evolution of eq. (8) in Fig. 5(a). The particle is trapped in the potential well and is oscillating in a region of $`23<x<28`$ for the initial velocity $`v=0.015`$ in Fig. 5(b). The trapping occurs for $`v<0.019`$ in eq. (8). The dissipative effect seems to be slightly weaker in eq. (8) than in the nonlinear Schrödinger equation eq. (1).

As the third case, we study soliton motion in a spatially periodic potential $`U(x)=U_0\mathrm{cos}(2\pi x)`$. Initially, a soliton locates at the position of the potential minimum $`x=L/2=25`$, and the amplitude is $`A=2`$. Figure 6 displays the time evolutions of solitons in the periodic potential with $`U_0=0.2`$ for two different initial velocities $`v`$. A soliton is almost steadily propagating for $`v=0.33`$ as shown in Fig. 6(a). For $`v=0.28`$, a soliton propagates in an initial stage, however, it is trapped in the periodic potential as shown in Fig. 6(b). The critical velocity for a soliton to be trapped for $`A=2`$ is approximately 0.312.

Figure 7(a) displays the critical kinetic energy $`E_c=v_c^2/2`$ as a function of $`A`$. The solid curve represents the critical kinetic energy, and the rhombi represent critical points below which a soliton cannot go over the nearest potential hump. (For the parameter of $`v`$ shown in Fig. 6(b), the soliton goes over several potential humps and then it is trapped.) The effective potential $`U_{eff}`$ for the periodic potential is expressed as

$$U_{eff}=_{\mathrm{}}^{\mathrm{}}\{U_0\mathrm{cos}(2\pi x)\}|\varphi |^2𝑑x/N=\frac{U_0\mathrm{cos}(2\pi \xi )\pi ^2a}{\mathrm{sinh}(\pi ^2a)}.$$

A particle, which obeys eq. (6) with the effective potential $`U_{eff}`$ and starts from a position of potential minimum, exhibits a transition from a trapped state to a propagating state. The critical initial kinetic energy is $`2U_0\pi ^2/(A\mathrm{sinh}(\pi ^2/A))`$, which is drawn by the dotted curve in Fig. 7(a). The critical curve is a good approximation for the critical values marked with rhombi. This is probably because eq. (6) is a good approximation to the short-time evolution of eq. (1), where dissipative effects are negligible, and which determines whether the soliton can pass over the first neighboring hump. In the region $`1.5<A<3`$, the critical curve $`2U_0\pi ^2/(A\mathrm{sinh}(\pi ^2/A))`$ is rather different from the true critical curve denoted by solid line. If the width of a soliton $`2/A`$ is comparable to the wavelength of the periodic potential, the interaction between the soliton and the potential will be enhanced, which may be an origin of the deviation. The trapping phenomena accompanying the dissipation of kinetic energy can be expressed with eq. (8) with $`U_{eff}=U_0\mathrm{cos}(2\pi \xi )\pi ^2a/\mathrm{sinh}(\pi ^2a)`$. The time evolution of $`\xi (t)`$ by eq. (8) for $`A=2`$ and $`v=0.2785`$ is displayed in Fig. 7(b) with the dashed curve. The solid curve in Fig. 7(b) is the time evolution of the peak position $`X_p(t)`$ of the soliton for $`A=2`$ and $`v=0.28`$. Similar time evolutions including damping of kinetic energy, trapping and chaotic motions are observed for $`\xi (t)`$ and $`X_p(t)`$, although the initial velocities are slightly different. Since the number of degrees of freedom in eq. (8) is four, it is not surprising that chaotic motion appears in eq. (8). The dashed curve in Fig. 7(a) is the critical line, above which steadily propagating solutions appear in eq. (8). There is fairly good agreement between the solid line and the dashed line, although the dashed line locates slightly below the sold line in the region $`1.5<A<3`$.

In summary, we have performed numerical simulations for the scattering and the trapping of nonlinear Schrödinger solitons in some external potentials. The transition from the reflection to the transmission by a potential wall depends strongly on the soliton amplitude. We have also found the trapping of solitons in a potential well and a periodic potential. The critical energy for the transitions can be approximately evaluated using eq. (6) or (8) derived from a variational method.


format: text

text:
# Dissipative tunneling and orthogonality catastrophe in molecular transistors

## I Introduction

The possibility of creating devices on a molecular level has opened up the field of single-molecule electronics in recent years. The transport properties of such mesoscopic systems have been investigated in a number of experiments. Reed et al. (1997); Park et al. (2000, 2002); Liang et al. (2002); Smit et al. (2002); Zhitenev et al. (2002); Pasupathy et al. ; Yu and Natelson (2004) Of particular interest has been the influence of strong electron-phonon coupling on electron transport, Park et al. (2000, 2002); Zhitenev et al. (2002) which manifests itself as emission and absorption of vibrational quanta observable in the excitation spectra. One example is the series of experiments by Park et al. Park et al. (2000) where it was shown that the current through a single $`C_{60}`$ molecule was strongly coupled to a single vibrational mode.

A large amount of theoretical works has dealt with the problem of tunneling through a single molecular electronic level that is coupled to phonon modes. Since in many experimental realizations the tunnel coupling of the molecule to the leads is rather weak compared to the other energy scales in the problem, transport is dominated by Coulomb-blockade effects and occurs sequentially. An approach based on rate equations can thus be justified and has been used in a number of recent papers. Boese and Schoeller (2001); McCarthy et al. ; Mitra et al. (a); Braig and Flensberg (2003) Physically, it is an essential question how the excited vibrational levels are allowed to relax, either through coupling to the environment, for example the phonons or plasmons of the metal substrate, or by virtue of the tunneling electrons. Aji et al. ; Flensberg (2003); Mitra et al. (b) In the case where the relaxation of the vibrational mode is faster than the tunneling rate one can assume an equilibrium phonon distribution.

The coupling between the vibrational mode of the molecule and the environment depends strongly on which vibrational mode is considered. For intra-molecular vibrations the lifetime can be very long. Gurevich and Schober (1998); Patton and Geller  However, we were able to show in a recent paper Braig and Flensberg (2003) that the vibrational mode associated with the center of mass motion of the molecule is coupled more strongly to the environment and is thus exposed to an effective damping mechanism. The electron-phonon effects seen in the experiments of Ref. Park et al., 2000 are suggested to be due to such a center of mass motion.

In the present paper, we investigate how screening of the charge on the molecule influences the damping mechanism. A sketch of the physical setup that serves as starting point is shown in Fig. 1.

If the molecule is not charged, it is held in place by an interaction of the van-der-Waals type. It can perform center-of-mass oscillations about its equilibrium position, and the characteristic energy $`\mathrm{}\omega `$ of such an oscillation was found to be $``$ 5meV. Park et al. (2000) These oscillations also influence the substrate to which the molecule is attached, and we will call this interaction the *van-der-Waals* interaction in the following. On the other hand, if the molecule is charged, there will be additional electric forces on the molecule, which are either due to static electric charges on the surface, e.g. because of impurities, or due to surface charges in the substrate that are induced by the charge on the molecule. Hence, when the electron tunnels onto the molecule, it results not only in a force on the molecule itself but also in electrostatic forces on the substrate to which it is attached, and the equilibrium positions of the vibrational modes in the substrate are shifted accordingly. In the following, we will focus on the case of surface charges only, but static impurity charges can be treated analogously and have the very same qualitative effect. We will refer to this type of interaction as the *electrostatic* interaction. foo (a)

Under the influence of the Coulomb interaction, the surface is attracted to the molecule (and vice versa), and substrate atoms will be shifted away from their equilibrium positions until van-der-Waals forces and electrostatic forces are in equilibrium. This shift, however, occurs in an effective dipole field (or multipole field for charges that are not pointlike) since the field between molecule and substrate due to the additional charge on the molecule is equivalent to the case without surface charges but an image charge located inside the substrate, cmp. Fig. 2. Only in the limit of a very short ranged van-der-Waals interaction compared to the separation of charge and image charge could we approximate the electric field at the surface as constant, and thus monopolar. Such a mismatch of length scales is rather unlikely for a situation as depicted in Figs. 1 and 2(a). However, a different experimental geometry as in Fig. 2(b) can introduce monopolar effects. Braig and Flensberg (2003)

The shift of the equilibrium positions $`𝐚_𝐤`$ for each phonon mode $`𝐤`$ can be calculated by resorting to the dynamical matrix $`𝐃(𝐤,𝐤^{})`$ which appears in the harmonic expansion of the Hamiltonian (see also Ref. Louis and Sethna, 1995):

$$H_{\mathrm{harm}}=\frac{1}{2}\underset{𝐤𝐤^{}}{}𝐚_𝐤𝐃(𝐤,𝐤^{})𝐚_𝐤^{}$$
(1)
If we denote the electrostatic force on the substrate as $`𝐆(𝐫)`$ and choose a phonon eigenmode basis such that $`𝐃(𝐤,𝐤^{})`$ is diagonal, we can see that force equilibrium demands

$$𝐃(𝐤)𝐚_𝐤\underset{𝐫}{}𝐆(𝐫)\mathrm{e}^{i\mathrm{𝐤𝐫}}$$
(2)
where the right hand side is the Fourier component of the force $`𝐆`$ acting on the plane wave phonon state $`𝐤`$. The long-wavelength behavior can give rise to a version of orthogonality catastrophe, i.e. zero overlap between the ground states corresponding to the shifted and unshifted equilibrium positions. For small wavevectors $`k`$, we know that $`𝐃(k)k^2`$, see Ref. Ashcroft and Mermin, 1976. Expanding $`𝐆`$ in a multipole expansion, and expanding the exponential for small $`k`$, we can see that $`\mathrm{a}_𝐤k^1`$ for a dipole force since the volume integral of any multipole field is zero by definition, Jackson (1999) whereas $`\mathrm{a}_𝐤k^2`$ for a monopole field since this corresponds to the constant part in the multipole expansion. This difference in long-wavelength behavior gives rise to nonzero (dipole) or zero (monopole) overlap between the ground states in three dimensions, as one can see from the following simple argument. Let $`\eta `$ be the ground-state overlap between a set of harmonic quantum oscillators and its shifted counterparts, then

$`\eta \mathrm{exp}\left({\displaystyle \underset{𝐤}{}}𝐚_𝐤^2/2\mathrm{}_𝐤^2\right)`$ (3)
where $`\mathrm{}_𝐤^2=\mathrm{}/m_𝐤\omega _𝐤`$ is the oscillator quantum with $`\mathrm{}_𝐤^2k^1`$ for small $`k`$. Converting the sum to an integration, we see that the question of zero or nonzero overlap is decided by both dimensionality and the small-$`k`$ behavior of the equilibrium shifts $`𝐚_𝐤`$. In particular,

$`\eta \mathrm{exp}\left({\displaystyle dkk^d𝐚_𝐤^2/A}\right),`$ (4)
where $`d`$ is the dimensionality, and $`A`$ is a normalization factor. Thus, in one dimension, we encounter orthogonality catastrophe for both monopolar and dipolar shifts, whereas in two and three dimensions we see that dipolar equilibrium shifts feature nonzero overlap!

It is this effect that we investigate in the present paper in the context of how a single vibrational mode of the molecule dissipates energy into a phononic environment. Again referring to Fig. 2, if the range of the van-der-Waals interaction $`D`$ is comparable with the dipole length scale (set by the separation of charge and image charge), we expect nonzero overlap between the two ground states in accord with a dipolar equilibrium shift, and we should thus see a step in the $`I`$-$`V`$ characteristics, cmp. Fig. 3. If, on the other hand, the van-der-Waals range is very small compared to the dipole length scale, we expect an orthogonality catastrophe corresponding to the approximately constant (monopolar) force. The latter can also occur for particular geometries, cmp. Fig. 2(b). The current then behaves according to a power law at the onset of conduction. Braig and Flensberg (2003)

Another way of looking at the same problem is in terms of translational invariance. Suppose that the van-der-Waals range $`D`$ and the dipole length scale are comparable to each other. In this case, we could think of drawing a box around the whole system which contains both the van-der-Waals and the electrostatic interaction. Hence, all forces are internal, and we have a translationally invariant Hamiltonian. Small $`D`$, on the other hand, corresponds to the case that, from the point of view of the van-der-Waals interaction, the electrostatic interaction is extremely extended and spread out. Therefore, the electrostatic force acts approximately as an *external* force as far as the van-der-Waals interaction is concerned and translational invariance is broken. This case can be recovered in a physical setup where the molecule is located asymmetrically in between two substrates, cmp. Fig. 2(b). Then, image charges in the substrate farther away can approximately have the effect of an external force on the substrate closer to the molecule.

The paper is organized as follows. The model Hamiltonian is defined in Sec. II, and in Sec. III we recapitulate the expression for the current derived on the basis of rate equations in Ref. Braig and Flensberg, 2003. The function that describes the tunneling density of states is then solved in the case of equal spread of van-der-Waals and electrostatic interaction in the presence of the dissipative environment in Sec. IV. Finally, a discussion and summary of the results can be found in Sec. V.

## II Model Hamiltonian

We consider a model of one single spin-degenerate molecular level coupled to two leads. The single level is coupled to the vibrational mode of the molecule through the charge on the molecule, as are the vibrational modes of the substrate. The coupling between the oscillator and the environment is included as a linear coupling to a bath of harmonic oscillators. Similar to the model studied in our previous paper, Braig and Flensberg (2003) the model Hamiltonian then reads

$`H`$ $`=H_{LR}+H_D+H_B+H_{\mathrm{bath}}+H_{DB}`$
$`+H_{D\mathrm{bath}}+H_{B\mathrm{bath}}+H_T,`$ (5)
with

$`H_{LR}`$ $`={\displaystyle \underset{k\sigma ,\alpha =L,R}{}}\xi _{k\alpha }c_{k\sigma ,\alpha }^{}c_{k\sigma ,\alpha }^{}`$ (6a)
$`H_D`$ $`={\displaystyle \underset{\sigma }{}}\xi _0d_\sigma ^{}d_\sigma ^{}+Un_dn_d,`$ (6b)
$`H_\mathrm{B}`$ $`={\displaystyle \frac{p_0^2}{2m_0}}+{\displaystyle \frac{1}{2}}m_0\omega _0^2x_0^2,`$ (6c)
$`H_{DB}`$ $`=\lambda x_0{\displaystyle \underset{\sigma }{}}d_\sigma ^{}d_\sigma ,`$ (6d)
$`H_{\mathrm{bath}}`$ $`={\displaystyle \underset{j}{}}\left({\displaystyle \frac{p_j^2}{2m_j}}+{\displaystyle \frac{1}{2}}m_j\omega _j^2x_j^2\right),`$ (6e)
$`H_{B\mathrm{bath}}`$ $`={\displaystyle \underset{j}{}}\beta _jx_jx_0,`$ (6f)
$`H_T`$ $`={\displaystyle \underset{k\sigma ,\alpha =L,R}{}}t_{k\alpha }c_{k\sigma ,\alpha }^{}d_\sigma +h.c.`$ (6g)
where the $`c_{k\sigma ,\alpha }^{},c_{k\sigma ,\alpha }^{}`$ and the $`d_\sigma ^{},d_\sigma ^{}`$ are creation and annihilation operators for the leads and the molecule, respectively, $`x_0`$ is the oscillator degree of freedom, $`\{x_j\}`$ describes the set of environmental degrees of freedom and $`m_j`$, $`\omega _j`$ their respective masses and frequencies, $`\xi _0`$ is the onsite energy, and $`U`$ is comprised of the Coulomb interaction and the energy cost for having double occupancy on the molecule. The coupling constants for the electron-oscillator interaction and the oscillator-bath interaction are $`\lambda `$ and $`\{\beta _j\},`$ respectively. The lead electron energies are given by
$$\xi _{k\alpha }=\epsilon _{k\alpha }\mu _\alpha ,$$
(7)
where $`\mu _\alpha `$ is the chemical potential of lead $`\alpha `$. The tunneling amplitudes $`t_{k\alpha }`$ in the tunneling term $`H_T`$ could in principle also depend on the oscillator position, which should be a small effect for the experimental realizations in Refs. Park et al., 2000 and Park et al., 2002Braig and Flensberg (2003)

We can explicitly write down how the $`\beta _j`$ are connected to the displacement $`u(𝐫)`$ perpendicular to the substrate if we assume that the force acting between the molecular vibrational mode and the substrate phonons is spread over the substrate according to a normalized distribution $`f(𝐫)`$, i.e. $`d𝐫f(𝐫)=1`$. If we model the interaction as a harmonic spring potential, we can write

$$H_{\mathrm{int}}=\frac{1}{2}m_0\omega _0^2\left[x_0d𝐫f(𝐫)u(𝐫)\right]^2.$$
(8)
Expanding the square, we can combine the term $`\frac{1}{2}m_0\omega _0^2\left[d𝐫f(𝐫)u(𝐫)\right]^2`$ with the part of the Hamiltonian that describes the “free” environmental bath and diagonalize the combination by introducing new coordinates $`x_j`$. These are related to $`u(𝐫)`$ by a unitary transformation as

$$u(𝐫)=\underset{j}{}c_j(𝐫)x_j,$$
(9)
with appropriately chosen coefficient functions $`c_j(𝐫)`$. Defining

$$\beta _jm_0\omega _0^2d𝐫f(𝐫)c_j(𝐫),$$
(10)
we immediately see that we recover the structure of our original Hamiltonian $`H`$. The advantage of the form of $`H_{\mathrm{int}}`$ is that its translational invariance with respect to shifting both $`x_0`$ and $`u(𝐫)`$ by the same constant displacement is immediately apparent since the distribution function $`f(𝐫)`$ is normalized to 1.

The difference to Ref. Braig and Flensberg, 2003 is the newly included coupling term $`H_{D\mathrm{bath}}`$ between the charge on the molecule and the bath of vibrational modes,

$`H_{D\mathrm{bath}}={\displaystyle \underset{j}{}}\lambda _jx_j{\displaystyle \underset{\sigma }{}}d_\sigma ^{}d_\sigma ^{}.`$ (11)
It is again instructive to rewrite this term as a function of $`u(𝐫)`$ and a normalized distribution function $`g(𝐫)`$ that characterizes the form and range of the electrostatic interaction between molecule and substrate:

$`H_{D\mathrm{bath}}=\lambda {\displaystyle d𝐫g(𝐫)u(𝐫)\underset{\sigma }{}d_\sigma ^{}d_\sigma ^{}},`$ (12)
so that the electrostatic force on the substrate $`𝐆(𝐫)`$ mentioned in Sec. I has the form $`𝐆(𝐫)=\lambda g(𝐫)_\sigma d_\sigma ^{}d_\sigma ^{}`$. Since $`u(𝐫)`$ can be written as a linear superposition of the eigenmodes $`x_j`$, see Eq. (9), we can express the coupling constants $`\lambda _j`$ in terms of the distribution function $`g(𝐫)`$ as

$`\lambda _j=\lambda {\displaystyle d𝐫g(𝐫)c_j(𝐫)}.`$ (13)
Given Eq. (12), translational invariance of the total Hamiltonian with respect to shifting both $`x_0`$ and $`u(𝐫)`$ by some constant displacement vector is ensured automatically if $`g(𝐫)`$ is normalized properly, i.e. $`d𝐫g(𝐫)=1`$.

We now want to relate the coupling constants in the boson-bath coupling to the finite damping of the vibrational mode, which can be accomplished by studying the classical equations of motion. After removing the bath degrees of freedom, we end up with the following equation of motion in the frequency domain, neglecting the coupling terms $`H_{DB}`$ and $`H_{D\mathrm{bath}}`$ which will be removed by a unitary transformation below:

$`[\omega ^2\omega _0^2𝒮(\omega )]x_0(\omega )=0,`$ (14)
where we have defined

$`𝒮(\omega )={\displaystyle \frac{1}{m_0}}{\displaystyle \underset{j}{}}{\displaystyle \frac{\beta _j^2}{m_j}}{\displaystyle \frac{1}{(\omega +i\eta )^2\omega _j^2}},`$ (15)
which is complex in general and gives rise to frictional damping and a frequency shift of the bare frequency $`\omega _0`$. In Ref. Braig and Flensberg, 2003, we explicitly calculated $`𝒮(\omega )`$ for the case of a molecule attached to a semi-infinite substrate. We will use these results throughout this paper, again assuming that the van-der-Waals force $``$ exerted by the molecule, cmp. Fig. 1(a), is directed perpendicular to the substrate surface and cylindrically symmetric. Also, we assume that both the van-der-Waals and the electrostatic forces only couple to the substrate surface, i.e., $`g(𝐫),f(𝐫)\delta (x)`$, where the $`x`$ axis is perpendicular to the substrate surface. The van-der-Waals force $``$ on the molecule is thus understood to be of the form

$`=m_0\omega _0^2\left[x_0{\displaystyle _0^{\mathrm{}}}2\pi rf(r)u_x(r)𝑑r\right],`$ (16)
where $`u_x(r)`$ is the perpendicular displacement on the surface of the substrate, and $`r`$ is the radial direction in the surface plane. For explicit numerical evaluations later in this paper, we will use $`f(r)=\mathrm{e}^{r/D}/2\pi D^2`$ just as in our previous paper, with $`D`$ on the order of the width $`D_0`$ of the molecule. Braig and Flensberg (2003)

We now eliminate the coupling terms $`H_{DB}`$ and $`H_{D\mathrm{bath}}`$ of the Hamiltonian (5) by a unitary transformation similar to the one used in the independent boson model, Mahan (1990) at the cost of introducing displacement operators in the tunneling term. However, since we have a somewhat more complicated system here due to the coupling to the bosonic bath, the unitary transformation in Ref. Mahan, 1990 has to be generalized. We define the transformation

$$\stackrel{~}{H}=SHS^{},S=e^{iAn_d},A=p_0\mathrm{}+\underset{j}{}p_j\mathrm{}_j,$$
(17)
where $`n_d=_\sigma d_\sigma ^{}d_\sigma `$. Using that

$$\stackrel{~}{x}_0=x_0\mathrm{}n_d,\stackrel{~}{x}_j=x_j\mathrm{}_jn_d,$$
(18)
both linear coupling terms, $`H_{DB}`$ and $`H_{D\mathrm{bath}}`$, cancel if we set

$`\mathrm{}={\displaystyle \frac{\lambda _j\frac{\beta _j\lambda _j}{m_j\omega _j^2}}{m_0[\omega _0^2+𝒮(0)]}},\mathrm{}_j={\displaystyle \frac{\lambda _j}{m_j\omega _j^2}}{\displaystyle \frac{\mathrm{}\beta _j}{m_j\omega _j^2}}`$ (19)
The Hamiltonian then transforms into

$$\stackrel{~}{H}=H_{LR}+\stackrel{~}{H}_D+H_B+H_{\mathrm{bath}}+H_{B\mathrm{bath}}+\stackrel{~}{H}_T,$$
(20)
where

$$\stackrel{~}{H}_T=\underset{k\sigma ,\alpha =L,R}{}t_{k\sigma ,\alpha }\left(c_{k\sigma ,\alpha }^{}e^{iA}d_\sigma +d_\sigma ^{}e^{iA}c_{k\sigma ,\alpha }\right)$$
(21)
and

$$\stackrel{~}{H}_D=\epsilon _0\underset{\sigma }{}d_\sigma ^{}d_\sigma +\stackrel{~}{U}n_dn_d,\epsilon _0=\xi _0\frac{1}{2}\lambda \mathrm{}.$$
(22)
Here we have defined

$`\stackrel{~}{U}`$ $`U{\displaystyle \frac{\kappa }{m_0[\omega _0^2+𝒮(0)]}},\epsilon _0\xi _0{\displaystyle \frac{1}{2}}{\displaystyle \frac{\kappa }{m_0[\omega _0^2+𝒮(0)]}},`$ (23)
where

$`\kappa `$ $`{\displaystyle \frac{1}{2}}\lambda ^2+{\displaystyle \frac{1}{2}}{\displaystyle \underset{j}{}}{\displaystyle \frac{\lambda _j^2}{m_j\omega _j^2}}\left(1{\displaystyle \underset{ij}{}}{\displaystyle \frac{\beta _i^2}{m_i\omega _i^2}}\right)`$
$`{\displaystyle \underset{j}{}}{\displaystyle \frac{\beta _j\lambda _j}{m_j\omega _j^2}}\left(\lambda {\displaystyle \underset{i>j}{}}{\displaystyle \frac{\beta _i\lambda _i}{m_i\omega _i^2}}\right)`$ (24)
Thus, the Coulomb repulsion $`U`$ is modified by both the phonon mediated interaction and the electrostatic interaction between molecule and substrate.

It is important to note that for the particular choice $`f(𝐫)=g(𝐫)`$, i.e. when the distribution functions of electrostatic and van-der-Waals interaction are equal, we obtain from Eqs. (10) and (13):

$`\lambda _j={\displaystyle \frac{\lambda \beta _j}{m_0\omega _0^2}},\mathrm{}={\displaystyle \frac{\lambda }{m_0\omega _0^2}},\mathrm{}_j=0.`$ (25)
Furthermore, we can quantify the deviation of $`f(𝐫)`$ from $`g(𝐫)`$ by defining

$$ϵ(𝐫)g(𝐫)f(𝐫),ϵ_jm_0\omega _0^2d𝐫ϵ(𝐫)c_j(𝐫),$$
(26)
where $`d𝐫ϵ(𝐫)=0`$. We can thus rewrite $`\lambda _j`$ as

$$\lambda _j=\frac{\lambda }{m_0\omega _0^2}\left(\beta _j+ϵ_j\right).$$
(27)
This immediately leads to

$`\mathrm{}={\displaystyle \frac{\lambda }{m_0\omega _0^2}}\left[1\mathrm{\Delta }(0)\right],\mathrm{}_j={\displaystyle \frac{\mathrm{}}{m_j\omega _j^2}}{\displaystyle \frac{ϵ_j+\mathrm{\Delta }(0)\beta _j}{1\mathrm{\Delta }(0)}},`$ (28)
where we defined

$`\mathrm{\Delta }(\omega ){\displaystyle \frac{1}{m_0\overline{\omega }_0^2}}{\displaystyle \underset{n}{}}{\displaystyle \frac{\beta _nϵ_n}{m_n}}{\displaystyle \frac{1}{(\omega +i\eta )^2\omega _n^2}}.`$ (29)
If the two distribution functions $`f(𝐫)`$ and $`g(𝐫)`$ deviate only slightly from each other, $`ϵ_j`$ is small compared to $`\beta _j`$, and we see that $`\mathrm{\Delta }1`$. Also note that the qualitative frequency dependence of $`\mathrm{\Delta }(\omega )`$ is exactly the same as of $`𝒮(\omega )`$.

## III Current formula

As mentioned in the Introduction, the most important assumption here is that the tunneling rate is much smaller than all other time scales, which means that we can assume the vibrational degrees of freedom and the Fermi seas in the two electrodes to be in equilibrium at all times. For simplicity, we consider only two charge states and therefore let $`U=\mathrm{}`$, which leaves us with only three states playing a role in the rate equations: empty, and occupied by either spin up or down. Treating the tunneling term $`H_T`$ as a perturbation and accordingly employing Fermi’s Golden Rule for the tunneling rates, the current $`I`$ through the molecule is given by

$$I=2e\frac{\mathrm{\Gamma }_{10}^R\mathrm{\Gamma }_{01}^L\mathrm{\Gamma }_{01}^R\mathrm{\Gamma }_{10}^L}{(\mathrm{\Gamma }_{01}^L+\mathrm{\Gamma }_{01}^R)+2(\mathrm{\Gamma }_{10}^L+\mathrm{\Gamma }_{10}^R)}.$$
(30)
Here, we defined $`\mathrm{\Gamma }_{ij}^\alpha `$ as the tunneling rate for tunneling from the state with occupation $`i`$ to the state occupation $`j`$ through lead $`\alpha `$ (see, for example, Ref. dev, ):

$`\mathrm{\Gamma }_{10}^\alpha `$ $`=\mathrm{\Gamma }_\alpha {\displaystyle \frac{d\omega }{2\pi }F(\omega )n_\alpha (\epsilon _0+\omega )},`$ (31a)
$`\mathrm{\Gamma }_{01}^\alpha `$ $`=\mathrm{\Gamma }_\alpha {\displaystyle \frac{d\omega }{2\pi }F(\omega )(1n_\alpha (\epsilon _0+\omega ))},`$ (31b)
with
$$F(\omega )=_{\mathrm{}}^{\mathrm{}}𝑑tF(t),F(t)=e^{iA(t)}e^{iA},$$
(32)
where the operator $`A`$ is defined in Eq. (17). Furthermore, we work in the wide-band limit with bare tunneling rates $`\mathrm{\Gamma }_\alpha `$=$`2\pi _k|t_{k\alpha }|^2\delta (\xi _k)`$, and we intoduced the Fermi distributions of the two leads $`n_\alpha (\epsilon )`$=$`(e^{\beta (\epsilon eV_\alpha )}+1)^1`$.

## IV Influence of the coupling to the environment

In presence of coupling to the environment, the evaluation of the function $`F(t)`$ in Eq. (32) is in principle straightforward since the Hamiltonian is quadratic in the oscillator and bath degrees of freedom. We obtain

$$F(t)=\mathrm{exp}\left(B(t)B(0)\right),B(t)=A(t)A(0)_0,$$
(33)
where the operator $`A`$ is defined in Eqs. (17) and (19). The expectation value $`\mathrm{}_0`$ is to be evaluated with respect to $`\stackrel{~}{H}`$ without the tunneling term. We follow the same procedure as in our previous paper Braig and Flensberg (2003) and use the fluctuation-dissipation theorem,

$$B(\omega )=2\mathrm{Im}[B^R(\omega )](1+n_B(\omega )),$$
(34)
to express $`B(t)`$ in terms of the corresponding retarded Green’s function

$$B^R(t)=i\theta (t)[A(t),A(0)]_0.$$
(35)
Here, $`n_B(\omega )=(e^{\beta \omega }1)^1`$ is the usual Bose function. We obtain

$`B^R(\omega )=`$ $`\mathrm{}^2m_0\left({\displaystyle \frac{\omega ^2}{\omega ^2\overline{\omega }_0^2\overline{𝒮}(\omega )}}1\right)`$
$`+2{\displaystyle \frac{\mathrm{}^2m_0}{1\mathrm{\Delta }(0)}}{\displaystyle \frac{\mathrm{\Delta }(0)\overline{𝒮}(\omega )\overline{\omega }_0^2\overline{\mathrm{\Delta }}(\omega )}{\omega ^2\overline{\omega }_0^2\overline{𝒮}(\omega )}}`$
$`+{\displaystyle \frac{\mathrm{}^2m_0}{[1\mathrm{\Delta }(0)]^2}}\left[\mathrm{\Delta }(0)^2\overline{𝒮}(\omega )2\mathrm{\Delta }(0)\overline{\mathrm{\Delta }}(\omega )\right]`$
$`+{\displaystyle \frac{\mathrm{}^2m_0}{[1\mathrm{\Delta }(0)]^2}}{\displaystyle \underset{j}{}}{\displaystyle \frac{ϵ_j^2}{m_j\omega _j^2}}{\displaystyle \frac{1}{(\omega +i\eta )^2\omega _j^2}}`$
$`+{\displaystyle \frac{\mathrm{}^2m_0}{\omega ^2[1\mathrm{\Delta }(0)]^2}}{\displaystyle \frac{[\mathrm{\Delta }(0)\overline{𝒮}(\omega )\overline{\omega }_0^2\overline{\mathrm{\Delta }}(\omega )]^2}{\omega ^2\overline{\omega }_0^2\overline{𝒮}(\omega )}},`$ (36)
where we have defined $`\overline{𝒮}(\omega )𝒮(\omega )𝒮(0)`$, $`\overline{\mathrm{\Delta }}(\omega )\mathrm{\Delta }(\omega )\mathrm{\Delta }(0)`$, and the experimentally observable renormalized frequency $`\overline{\omega }_0^2\omega _0^2+𝒮(0)`$.

If we choose a very wide-ranged distribution $`g(𝐫)`$, as compared to $`f(𝐫)`$, then the function $`ϵ(𝐫)`$ will be close to $`f(𝐫)`$, see Eq. (26). This leads to $`ϵ_j\beta _j`$, and upon insertion into Eq. (28) we recover the expressions for $`\mathrm{}`$ and $`\mathrm{}_j`$ that we found previously in Ref. Braig and Flensberg, 2003 since then $`\mathrm{\Delta }(\omega )𝒮(\omega )`$. In effect, the electrostatic interaction is thus so widespread compared to the van-der-Waals interaction that it can be considered as an external influence on the system. In this case, we are justified to neglect the term $`H_{D\mathrm{bath}}`$ in our Hamiltonian altogether since translational invariance is broken by this external force. It is a matter of straightforward algebra to show that $`B^R(\omega )`$ indeed reduces to the same expression as in Ref. Braig and Flensberg, 2003.

If, on the other hand, the spread of the two interactions over the surface is comparable to each other, i.e. $`g(𝐫)f(𝐫)`$, we have the situation where $`ϵ_j\beta _j`$, and $`\mathrm{\Delta }1`$ in Eq. (28). The first term in Eq. (36) is of zeroth order in $`ϵ_j`$, the second term is of first order, and all the remaining terms are of second order. Furthermore, it is important to see that the imaginary part of all terms except of the last one are convergent for small $`\omega `$ since the small-frequency dependence of Im$`\overline{\mathrm{\Delta }}(\omega )`$ is identical to that of Im$`\overline{𝒮}(\omega )\omega `$. However, the last term is of second order in $`ϵ_j`$, so that to leading order, Im$`B^R(\omega )`$ is convergent for small $`\omega `$ if the distribution functions $`f(𝐫)`$ and $`g(𝐫)`$ feature comparable ranges.

How small the $`ϵ_j`$ really have to be in comparison to the $`\beta _j`$ (i.e., how similar in range $`f(𝐫)`$ and $`g(𝐫)`$ have to be) in order for such a first-order expansion to make sense can be deduced from the following argument: from our previous analysis Braig and Flensberg (2003) we know that the divergent terms in $`\mathrm{Im}B^R(\omega )`$ for small $`\omega `$ behave like $`C/\omega `$, where $`C`$ is a numerical constant of second order in the $`ϵ_j`$. The divergence shows its predominant influence in the large-time behavior of $`F(t)`$:

$$F(t)=\mathrm{exp}\left(_{\mathrm{}}^{\mathrm{}}\frac{d\omega }{2\pi }(\mathrm{e}^{i\omega t}1)B(\omega )\right).$$
(37)
For large times $`t`$, the integral in the exponent diverges logarithmically, where the physical cutoff for the integration is given by the smallest frequency of the phonon spectrum $`\omega _{\mathrm{min}}\overline{\omega }_0`$foo (b) In order for a first-order expansion in the $`ϵ_j`$ to hold, we therefore have to demand that the following be true:

$$\mathrm{ln}\left[\frac{\overline{\omega }_0}{\omega _{\mathrm{min}}}\right]|C|1,\mathrm{i}.\mathrm{e}.\frac{\overline{\omega }_0}{\omega _{\mathrm{min}}}\mathrm{e}^{1/|C|}.$$
(38)
Even if $`ϵ_j\beta _j`$ is enforced only moderately, this condition is not very restrictive. We should therefore expect a step in the conductance even if $`f(𝐫)`$ and $`g(𝐫)`$ are more than just “slightly” different!

In the following, we will examine the particular case $`f(𝐫)=g(𝐫)`$, i.e. $`\lambda _j=\lambda \beta _j/m_0\omega _0^2`$, see Eq. (25), for which the small-frequency behavior changes quite dramatically. (This is the case most different from our previous analysis in Ref. Braig and Flensberg, 2003.) Then we have $`\mathrm{}_j=0`$, and thus

$`B^R(\omega )`$ $`=\mathrm{}^2m_0\left({\displaystyle \frac{\omega ^2}{\omega ^2\overline{\omega }_0^2\overline{𝒮}(\omega )}}1\right),`$ (39)
which leads to

$$B(\omega )=\frac{4g}{\overline{\omega }_0}\frac{\omega ^2[1+n_B(\omega )]\mathrm{Im}\overline{𝒮}(\omega )}{\left[\omega ^2\overline{\omega }_0^2\mathrm{Re}\overline{𝒮}(\omega )\right]^2[\mathrm{Im}\overline{𝒮}(\omega )]^2}.$$
(40)
Here, we have defined the experimental parameter

$$g\frac{1}{2}\mathrm{}^2m_0\overline{\omega }_0=\frac{1}{2}\frac{\lambda ^2\overline{\omega }_0}{m_0\omega _0^4},$$
(41)
which is a measure of how much the equilibrium position of the molecular vibrational mode is shifted compared to the oscillator quantum $`\sqrt{\mathrm{}/m_0\overline{\omega }_0}`$. Note the appearance of both the renormalized and bare frequency in Eq. (41) and thus Eq. (40).

One can then find $`F(t)`$ from Eq. (37). For large times $`t`$, the integral in the exponent is convergent and yields a nonzero finite number since Im$`\overline{𝒮}(\omega )\omega `$ for small frequencies in the case of coupling to a semi-infinite substrate. Braig and Flensberg (2003) In particular, we obtain

$`\underset{t\mathrm{}}{lim}F(t)=\mathrm{exp}\left({\displaystyle \frac{4g}{\overline{\omega }_0}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}{\displaystyle \frac{\omega ^2[1+n_B(\omega )]\mathrm{Im}\overline{𝒮}(\omega )}{𝒢(\omega )}}\right),`$ (42)
where we defined

$`𝒢(\omega )\left[\omega ^2\overline{\omega }_0^2\mathrm{Re}\overline{𝒮}(\omega )\right]^2[\mathrm{Im}\overline{𝒮}(\omega )]^2`$ (43)
for the convenience of shorter formulae. Thus, we expect a $`\delta `$-function at $`\omega =0`$, which in turn results in a finite steplike discontinuity in the $`I`$-$`V`$ curves at zero temperature according to Eqs. (30, 31). We define

$$f(t)=F(t)\underset{t\mathrm{}}{lim}F(t),$$
(44)
and thus obtain at zero temperature:

$`F(\omega )=`$ $`2\pi \left(\underset{t\mathrm{}}{lim}F(t)\right)\delta (\omega )+f(\omega )`$ (45)
with

$`f(\omega )=`$ $`{\displaystyle \frac{4g}{\overline{\omega }_0}}{\displaystyle \frac{\omega ^2\mathrm{Im}\overline{𝒮}(\omega )}{𝒢(\omega )}}\left(\underset{t\mathrm{}}{lim}F(t)\right)`$
$`{\displaystyle \frac{4g}{\omega \overline{\omega }_0}}{\displaystyle _0^\omega }{\displaystyle \frac{d\xi }{2\pi }}f(\xi ){\displaystyle \frac{(\omega \xi )^3\mathrm{Im}\overline{𝒮}(\omega \xi )}{𝒢(\omega \xi )}}.`$ (46)
This implies that $`f(\omega )\omega ^3`$ for small $`\omega `$. For small $`g`$, we can actually approximate $`f(\omega )`$ by the first term, since the other terms are of higher order in $`g`$.

The resulting discontinuity $`\mathrm{\Delta }I`$ in the $`I`$-$`V`$ characteristics calculated from Eq. (30) and measured in units of $`I_N=e\mathrm{\Gamma }_L\mathrm{\Gamma }_R/(\mathrm{\Gamma }_L+\mathrm{\Gamma }_R)`$ is given by

$$\frac{\mathrm{\Delta }I}{I_N}=2\mathrm{sign}(V)\left(\underset{t\mathrm{}}{lim}F(t)\right)\frac{1+\mathrm{\Gamma }_L/\mathrm{\Gamma }_R}{2+\mathrm{\Gamma }_L/\mathrm{\Gamma }_R}.$$
(47)
We explicitly calculate $`\mathrm{\Delta }I`$ for the case of a $`C_{60}`$ molecule attached to a gold substrate, which is shown in Fig. 4 as a function of the experimental parameter $`g`$. The curves are parametrized by the spread $`D/D_0`$ of the vibrational force over the surface of the substrate, see Eqs. (15, 42). However, the dependence of the jump $`\mathrm{\Delta }I`$ on $`D/D_0`$ is only weak, and the size is comparable to the step size encountered in the case without coupling to the environment, for which we had the same formula as above with $`lim_t\mathrm{}F(t)=\mathrm{e}^g`$Braig and Flensberg (2003) The discontinuity appears at the voltage $`eV=\pm 2\epsilon _0`$ for symmetric bias voltage $`V_{L,R}=\pm V/2`$, but is not present for $`\epsilon _0=0`$ for reasons of symmetry.

Unless suppressed for reasons of device geometry, such a discontinuity at onset of conduction should be visible experimentally for weak tunneling into and out of the molecule even when the effect of finite tunneling rates $`\mathrm{\Gamma }`$ is taken into account. Following the derivations by Meir and Wingreen, Meir and Wingreen (1992) the differential conductance line shape is given by $`A(\omega )\mathrm{\Gamma }`$, where $`A(\omega )`$ is the spectral function of the molecule. Taking the ground state phonon overlap $`\eta `$ from Eq. (3), it can be shown that Flensberg (2003)

$`A(\omega )\mathrm{\Gamma }{\displaystyle \frac{\mathrm{\Gamma }^2|\eta |^6}{\omega ^2+\mathrm{\Gamma }^2|\eta |^4/4}}.`$ (48)
Thus, the height of the differential conductance step scales as $`|\eta |^2`$, whereas the characteristic width of the step scales as $`\mathrm{\Gamma }|\eta |^2`$. For small enough $`\mathrm{\Gamma }`$ compared to the width of the first Franck-Condon step in the $`I`$-$`V`$ curve, we therefore see that we recover a sharp step-like discontinuity in the $`I`$-$`V`$ curve. We can estimate $`\mathrm{\Gamma }`$ by virtue of the maximum current through the molecule $`I_{max}e\mathrm{\Gamma }`$, which yields $`\mathrm{}\mathrm{\Gamma }0.1\mu `$eV since $`I_{max}`$ 100pA for the experimental setup by Park *et al.* Park et al. (2000) The Franck-Condon steps, on the other hand, feature a source-drain voltage width on the order of mV in the $`I`$-$`V`$ curves, see Ref. Park et al., 2000, so that the discontinuity in the $`I`$-$`V`$ curve will remain sharp even for finite tunneling.

## V Summary and discussion

We found that inclusion of the effects due to screening of the charge on the molecule can change the qualitative form of dissipation with which the molecular vibrational mode is faced. Such screening effects can be mediated by static surface charges on the substrate, or by surface charges that are induced in the substrate by the charge on the molecule. If the distribution functions that characterize the form of the electrostatic interaction and the van-der-Waals interaction of the molecule with the substrate are not too much different from each other, we encounter a step in the $`I`$-$`V`$ curves at the onset of conduction. Such a discontinuity should be observable experimentally in the weak tunneling limit even with smearing due to finite tunneling effects. This is a result qualitatively different from treating the electrostatic forces as an external influence on the molecule-substrate system. The step size is dependent on the experimental parameter $`g`$ that reflects the shift of the equilibrium position of the molecular mode compared to the oscillator quantum $`m_0\overline{\omega }_0`$ and features an exponential dependence on $`g`$ similar to what one would obtain for the case of no coupling to the environment. The exponent is weakly dependent on the spread $`D/D_0`$ of the van-der-Waals coupling between molecular mode and substrate modes, and we recover the decoupled limit for $`D/D_0\mathrm{}`$ as expected from our previous results in Ref. Braig and Flensberg, 2003.

However, we also find that the results of our previous analysis are not changed qualitatively (but can be changed quantitatively) if the electrostatic interaction is very long ranged compared to the van-der-Waals interaction, and we recover a power law for the current at onset of conduction. In particular, the limit $`\lambda _j0`$ taken in Ref. Braig and Flensberg, 2003 is explicitly seen to coincide with treating the electrostatic interaction as an external influence on the system that breaks translational invariance.

The data in Ref. Park et al., 2000 do not exhibit a step in the $`I`$-$`V`$ characteristics and thus feature agreement with the case of an external force as chosen in Ref. Braig and Flensberg, 2003. It is quite possible, however, that in other experimental realizations, van-der-Waals interaction and electrostatic interaction are of comparable range, without the geometry of the system suppressing the effects that we found to arise due to the absence of external forces. According to our calculations, one would then indeed observe a significant step in the $`I`$-$`V`$ curve at the onset of conduction.

###### Acknowledgements.

The authors would like to thank P. W. Brouwer, J. P. Sethna, and P. Sharma for discussions. The work was supported by the Cornell Center for Materials Research under NSF Grant No. DMR0079992, and by the Danish National Research Council.


text:
# Statistical Physics in Meteorology

## I  INTRODUCTION and FOREWORD

This contribution to the 18th Max Born Symposium Proceedings, cannot be seen as an extensive review of the connection between meteorology and various aspects of modern statistical physics. Space and time (and weather) limit its content. Much of what is found here can rather be considered to result from a biased view point or limited understanding of a frustrated new researcher unsatisfied by the present status of the field. Yet only to be found is a set of basic considerations and reflections expecting to give lines for various investigations, in the spirit of modern statistical physics ideas.

The author came into this subject starting from previous work in econophysics, when he observed that some ”weather derivatives” were in use, and some sort of game initiated by the Frankfurt Deutsche Börse in order to attract customers which could predict the temperature in various cities within a certain lapse of time, and win some prize thereafter. This subject was similar to predicting the S&P500 or other financial index values at a certain future time. Whence various techniques which were used in econophysics, like the detrended fluctuation analysis, the multifractals, the moving average crossing techniques, etc. could be attempted from scratch.

Beside the weather (temperature) derivatives other effects are of interest. Much is said and written about e.g. the ozone layer and the Kyoto ”agreement”. The El Ni$`\stackrel{~}{n}`$o system is a great challenge to scientists. Since there is some data available under the form of time series, like the Southern Oscillation Index, it is of interest to look for trends, coherent structures, periods, correlations in noise, etc. in order to bring some knowledge, if possible basic parameters, to this meteorological field and expect to import some modern statistical physics ideas into such climatological phenomena. It appeared that other data are also available, like those obtained under various experiments, put into force by various agencies, like the Atlantic Stratocumulus Transition Experiment (ASTEX) for ocean surfaces or those of the Atmospheric Radiation Measurement Program (ARM), among others.

However it appeared that the data is sometimes of rather limited value because of the lack of precision, or are biased because the raw data is already transformed through models, and arbitrarily averaged (”filtered”) whence even sometimes lacking the meaning it should contain. Therefore a great challenge comes through in order to sort out the wheat from the chaff in order to develop meaningful studies. I will mention most of the work to which I have contributed, being aware that I am failing to acknowledge many more important reports than those, - for what I truly apologize. There are very interesting lecture notes on the web for basic modules on meteorological training courses, e.g. one available through ECMWF website.

In Sect.2, I will briefly comment on the history of meteorology. The notion of clouds, in Sect. 3, allows for bringing up the geometrical notion of fractals for meteorology work, thus scaling laws, and modern data analysis techniques. Simple technical and useful approaches, based on standard statistical physics techniques and ideas, in particular based on the scaling hypothesis for phase transitions and percolation theory features will be found in Sect. 4.

## II  HISTORICAL INTRODUCTION

From the beginning of times, the earth, sky, weather have been of great concern. As soon as agriculture, commerce, travelling on land and sea prevailed, men have wished to predict the weather. Later on airborne machines need atmosphere knowledge and weather predictions for best flying. Nowadays there is much money spent on weather predictions for sport activities. It is known how the knowledge of weather (temperature, wind, humidity, ..) is relevant, (even $`fundamental`$ !), e.g. in sailing races or in Formula 1 and car rally races. Let it be recalled the importance of knowing and predicting the wind (strength and directions), pressure and temperature at high altitude for the (recent) no-stop balloon round the world trip. The first to draw sea wind maps was Halley, an admirer of Breslau administration. That followed the ”classical” isobaths and isoheights (these are geometrical measures !!!) for sailors needing to go through channels.

I am very pleased to point out that Heinrich Wilhelm Brandes(1777-1834), Professor of Mathematics and Physics at the University of Breslau was the first who had the idea of displaying weather data (temperature, air pressure, a.s.o.) on geographical mapsIt seems that H.W. Brandes left Breslau to get his Ph.D. thesis in Heidelberg in 1826. Alas it seems that the original drawings are not available at this time. Where are they?. Later von Humboldt (1769-1859) had the idea to connect points in order to draw isotherms. It is well known nowadays that various algorithms will give various isotherms, starting from the same temperature data and coordinate table. In fact the maximum or minimum temperature as defined in meteorology are far from the ones acceptable in physics laboratories. Note that displayed isotherms connect data points which values are obtained at different times! No need to say that it seems essential to concentrate on predicting the uncertainty in forecast models of weather and climate as emphasized elsewhere.

## III  CLIMATE and WEATHER. The role of clouds

Earth’s climate is clearly determined by complex interactions between sun, oceans, atmosphere, land and biosphere . The composition of the atmosphere is particularly important because certain gases, including water vapor, carbon dioxide, etc., absorb heat radiated from Earth’s surface. As the atmosphere warms up, it in turn radiates heat back to the surface that increases the earth’s ”mean surface temperature”.

Much attention has been paid recently  to the importance of the main components of the atmosphere, in particular clouds, in the water three forms — vapor, liquid and solid, for buffering the global temperature against reduced or increased solar heating. This leads to efforts to improve not only models of the earth’s climate but also predictions of climate change, as understood over long time intervals, in contrast to shorter time scales for weather forecast. In fact, with respect to climatology the situation is very complicated because one does not even know what the evolution equations are. Since controlled experiments cannot be performed on the climate system, one relies on using ad hoc models to identify cause-and-effect relationships. Nowadays there are several climate models belonging to many different centers. Their web sites not only carry sometimes the model output used to make images but also provide the source code. It seems relevant to point out here that the stochastic resonance idea was proposed to describe climatology evolution.

It should be remembered that solutions of Navier-Stokes equations forcefully depend on the initial conditions, and steps of integrations. Therefore a great precision on the temperature, wind velocity, etc. cannot be expected and the solution(s) are only looking like a mess after a few numerical steps. The Monte Carlo technique suggests to introduce successively a set of initial conditions, perform the integration of the differential equations and make an average thereafter. It is hereby time to mention Lorenz’s work who simplified Navier-Stokes equations searching for some predictiability. However, predicting the outcome of such a set of equations with complex nonlinear interactions taking place in an open system is a difficult task.

The turbulent character in the atmospheric boundary layer (ABL) is one of its most important features. Turbulence can be caused by a variety of processes, like thermal convection, or mechanically generated by wind shear, or following interactions influenced by the rotation of the Earth. This complexity of physical processes and interactions between them create a variety of atmospheric formations. In particular, in a cloudy ABL the radiative fluxes produce local sources of heating or cooling within the mixed-layer and therefore can greatly influence its turbulent structure and dynamics, especially in the cloud base. Two practical cases, the marine ABL and the continental ABL have been investigated for their scaling properties

Yet, let it be emphasized that the first modern ideas of statistical physics implemented on cloud studies through fractal $`geometry`$ are due to Lovejoy who looked at the perimeter-area relationship of rain and cloud areas, fractal dimension of their shape or ground projection. He discovered the statistical self-similarity of cloud boundaries through area-perimeter analyses of the geometry of satellites,fractal scaling of the cloud perimeter in the horizontal plane. He found the fractal dimension $`D_p4/3`$ over a spectrum of 4 orders of magnitude in size, for small fair weather cumuli ($``$ 1021 km) up to huge stratus fields ($``$ 103 km). Cloud size distributions have also been studied from a scaling point of view. Rain has also received much attenion.

## IV  Modern Statistical Physics Approaches

Due to the nonlinear physics laws governing the phenomena in the atmosphere, the time series of the atmospheric quantities are usually non-stationary as revealed by Fourier spectral analysis, - whih is usually the first technique to use. Recently, new techniques have been developed that can systematically eliminate trends and cycles in the data and thus reveal intrinsic dynamical properties such as correlations that are very often masked by nonstationarities,. Whence many studies reveal long-range power-law correlations in geophysics time series in particular in meteorology. Multi-affine properties  can also be identified, using singular spectrum or/and wavelets.

There are different levels of essential interest for sorting out correlations from data, in order to increase the confidence in predictability. There are investigations based on long-, medium-, and short-range horizons. The $`i`$-diagram variability ($`iVD`$) method allows to sort out some short range correlations. The technique has been used on a liquid water cloud content data set taken from the Atlantic Stratocumulus Transition Experiment (ASTEX) 92 field program . It has also been shown that the random matrix approach can be applied to the empirical correlation matrices obtained from the analysis of the basic atmospheric parameters that characterize the state of atmosphere. The principal component analysis technique is a standard technique in meteorology and climate studies. The Fokker-Planck equation for describing the liquid water path  is also of interest. See also some tentative search for power law correlations in the Southern Oscillation Index fluctuations characterizing El Ni$`\stackrel{~}{n}`$o. But there are many other works of interest.

### A Ice in cirrus clouds

In clouds, ice appears in a variety of forms, shapes, depending on the formation mechanism and the atmospheric conditions. The cloud inner structure, content, temperature, life time, .. can be studied. In cirrus clouds, at temperatures colder than about $`40^{}`$ C ice crystals form. Because of the vertical extent, ca. from about 4 to 14 km and higher, and the layered structure of such clouds one way of obtaining some information about their properties is mainly by using ground-based remote sensing instruments. Attention can be focussed on correlations in the fluctuations of radar signals obtained at isodepths of $`winter`$ and $`fall`$ cirrus clouds giving (i) the backscattering cross-section, (ii) the Doppler velocity and (iii) the Doppler spectral width of the ice crystals. They correspond to the physical coefficients used in Navier Stokes equations to describe flows, i.e. bulk modulus, viscosity, and thermal conductivity. It was found that power-law time correlations exist with a crossover between regimes at about 3 to 5 min, but also $`1/f`$ behavior, characterizing the top and the bottom layers and the bulk of the clouds. The underlying mechanisms for such correlations likely originate in ice nucleation and crystal growth processes.

### B Stratus clouds

In stratus clouds, long-range power-law correlations  and multi-affine properties have reported for the liquid water fluctuations, beside the spectral density. Interestingly, stratus cloud data retrieved from the radiance, recorded as brightness temperature,§§§http://www.phys.unm.edu/ duric/phy423/l1/node3.html at the Southern Great Plains central facility and operated in the vertically pointing mode indicated a Fourier spectrum, $`S(f)f^\beta `$, $`\beta `$ exponent equal to $`1.56\pm 0.03`$ pointing to a nonstationary time series. The detrended fluctuation analysis (DFA) method applied on the stratus cloud brightness microwave recording indicates the existence of long-range power-law correlations over a two hour time.

Contrasts in behaviors, depending on seasons can be pointed out. The DFA analysis of liquid water path data measured in April 1998 gives a scaling exponent $`\alpha =0.34\pm 0.01`$ holding from 3 to 60 minutes. This scaling range is shorter than the 150 min scaling range for a stratus cloud in January 1998 at the same site. For longer correlation times a crossover to $`\alpha =0.50\pm 0.01`$ is seen up to about 2 h, after which the statistics of the DFA function is not reliable.

However a change in regime from Gaussian to non-Gaussian fluctuation regimes has been clearly defined for the cloud structure changes using a finite size (time) interval window. It has been shown that the DFA exponent turns from a low value (about 0.3) to 0.5 before the cloud breaks. This indicates that the stability of the cloud, represented by antipersistent fluctuations is (for some unknown reason at this level) turning into a system for which the fluctuations are similar to a pure random walk. The same type of finding was observed for the so called Liquid Water PathThe liquid water path (LWP) is the amount of liquid water in a vertical column of the atmosphere; it is measured in cm-3; … sometimes in cm !!!.

The value of $`\alpha 0.3`$ can be interpreted as the $`H_1`$ parameter of the multifractal analysis of liquid water content and of liquid water path . Whence, the appearance of broken clouds and clear sky following a period of thick stratus can be interpreted as a non equilibrium transition or a sort of fracture process in more conventional physics. The existence of a crossover suggests two types of correlated events as in classical fracture processes: nucleation and growth of diluted droplets. Such a marked change in persistence implies that specific fluctuation correlation dynamics should be usefully inserted as ingredients in ad hoc models.

### C Cloud base height

The variations in the local $`\alpha `$-exponent (”multi-affinity”) suggest that the nature of the correlations change with time, so called intermittency phenomena. The evolution of the time series can be decomposed into successive persistent and anti-persistent sequences. It should be noted that the intermittency of a signal is related to existence of extreme events, thus a distribution of events away from a Gaussian distribution, in the evolution of the process that has generated the data. If the tails of the distribution function follow a power law, then the scaling exponent defines the critical order value after which the statistical moments of the signal diverge. Therefore it is of interest to probe the distribution of the fluctuations of a time dependent signal $`y(t)`$ prior investigating its intermittency. Much work has been devoted to the cloud base height, under various ABL conditions, and the LWP. Neither the distribution of the fluctuations of liquid water path signals nor those of the cloud base height appear to be Gaussian. The tails of the distribution follow a power law pointing to ”large events” also occurring in the meteorological (space and time) framework. This may suggest routes for other models.

### D Sea Surface Temperature

Other time series analysis have been investigated searching for power law exponents, like in atmospheric or sea surface temperature (SST) fluctuations. These are of importance for weighing their impacts on regional climate, whence finally to greatly increase predictability of precipitation during all seasons. Currently, climate patterns derived from global SST are used to forecast precipitation.

Recently we have attempted to observe whether the fluctuations in the Southern Oscillation index ($`SOI`$) characterizing El Ni$`\stackrel{~}{n}`$o were also prone to a power law analysis. For the $`SOI`$ monthly averaged data time interval 1866-2000, the tails of the cumulative distribution of the fluctuations of $`SOI`$ signal it is found that large fluctuations are more likely to occur than the Gaussian distribution would predict. An antipersistent type of correlations exist for a time interval ranging from about 4 months to about 6 years. This leads to favor specific physical models for El Ni$`\stackrel{~}{n}`$o description.

## V Conclusions

Modern statistical physics techniques for analyzing atmospheric time series signals indicate scaling laws (exponents and ranges) for correlations. A few examples have been given briefly here above, mainly from contributed papers in which the author has been involved. Work by many other authors have not been included for lack of space. This brief set of comments is only intended for indicating how meteorology and climate problems can be tied to scaling laws and inherent time series data analysis techniques. Those ideas/theories have allowed me to reduce the list of quoted references, though even like this I might have been unfair. One example can be recalled in this conclusion to make the point: the stratus clouds break when the molecule density fluctuations become Gaussian, i.e. when the molecular motion becomes Brownian-like. This should lead to better predictability on the cloud evolution and enormously extend the predictability range in weather forecast along the lines of nonlinear dynamics.

Acknowledgments

Part of this studies have been supported through an Action Concertée Program of the University of Li$`\stackrel{`}{e}`$ge (Convention 02/07-293). Comments by A. Pȩkalski, N. Kitova, K. Ivanova and C. Collette are greatly appreciated.


text:
# The Fate of Luminous Compact Blue Galaxies: An Environmental Approach

## 1. Introduction

A key prediction of CDM hierarchical formation models (e.g., White & Frenk 1991) is the “down-sizing” of star formation sites, whereby the characteristic galaxy mass dominating the co-moving star-formation rate decreases with redshift (Cowie et al. 1996). At high redshift, Lyman-break galaxies dominate; we believe they are associated with today’s high-mass systems (Lowenthal et al. 1997; Steidel et al. 1996). Today, the most intensely star-forming systems are low-mass HII galaxies. Between $`0.1<z<1.0`$ the situation is less clear.

“Luminous compact blue galaxies” (LCBGs) appear to be the most rapidly evolving class of galaxy at intermediate redshift (Lilly et al. 1998; Mallen-Ornelas et al. 1999). These enigmatic galaxies, initially identified by Koo & Kron (1988), are luminous ($`M_B20`$), small ($`R_e23`$kpc), and massive engines of star formation (up to $`40M_{}yr^1`$) (Koo et al. 1995, Guzmán et al. 1998, Hammer et al. 2000). They appear to form a link in redshift, size, and luminosity between Lyman-break galaxies and HII galaxies today (Lowenthal et al. 1997; see our Figure 1). LCBGs are a major contributor to the observed enhancement of the star-formation density of the universe at $`z<1`$, and their mass and number density decline in concert with the rapid drop in the global SFR since $`z=1`$ (Guzmán et al. 1998).

However, a debate continues about today’s descendants of LCBGs. Due to their narrow emission line-widths and number density, Koo et al. (1995) proposed LCBGs are the progenitors of today’s lower-mass spheroidal (Sphs) galaxies: If their current burst of star formation terminates, LCBGs fade by 4 mag in a few Gyrs, yielding luminosities, surface-brightnesses, sizes, and colors similar to, e.g., NGC 205 (Guzmán et al. 1998). LCBGs have also been proposed as the period when bulges of spiral galaxies form due to their high gas phase metallicities (Hammer et al. 1999, Kobulnicky & Zaritsky 1999). In this picture, the bright, blue compact regions are embedded within a lower-surface brightness, disk of larger size, (Barton & Van Zee 2001).

## 2. LCBGs in Galaxy Clusters: New Clues on LCBGs Descendants

Because of the extreme density of galaxy clusters, they offer a distinctly different, yet heretofore unexplored environment in which to study the nature of LCBGs. For our purpose, there are two salient differences in galaxy populations in rich clusters and the field: First, the morphology-density relationship describes a decreasing fraction of spirals with increasing local surface density of galaxies (Dressler 1980). Second, the dwarf-galaxy population is much richer in clusters. Dwarf galaxies are the most numerous type of galaxy regardless of environment, but the ratio of dwarf to giant galaxies increases with the overall density (e.g., Trentham and Hodgkins 2002). Because the relative densities of the two proposed descendants of LCBGs (lower-mass spheroidals and spiral bulges) are different in galaxy clusters and the field, the number density of LCBGs in intermediate galaxy clusters make the connection between the past (LCBGs) and today (either low-mass spheroidals or spirals). Specifically, if LCBGs are relatively more numerous in intermediate-redshift clusters than in the field, this provides independent evidence that LCBGs are associated with today’s lower-mass spheroidal population.

We currently are analyzing imaging data on 10 intermediate redshift galaxy clusters from the WIYN Long Term Variability Survey (WLTV) to search for cluster LCBGs. Each cluster has been repeatedly imaged over 5 years in a 10 arcmin field with the WIYN 3.5m telescope. We have a total of $`25`$ hrs integration per target in each of the UBRIz bands taken under excellent seeing conditions (typically $`<`$ 0.75”). Deep, archival HST images are available of each cluster core. We are in the process of gathering rest-frame OII and continuum narrow-band images of the 6 deepest clusters, four of which are visible to SALT.

We have completed the first stage of an LCBG search in cluster MS0451, at $`z=0.53`$. We frame our analysis in the measurement of galaxy enhancement. The enhancement is the number of galaxies at a given redshift divided by the expected number in the field at that redshift (estimated from comparable deep images without rich clusters). Hence, if there is no galaxy over-density, the enhancement is one. Figure 2 shows the enhancement of LCBGs in the field of MS0451. This is a robust, differential measurement, which, for this cluster appears to support the hypothesis that LCBGs evolve into Sphs.

We also have calculated the photometric and structural characteristics of the LCBGs appearing in the inner 0.75 Mpc of MS0451 using F702 WFPC2 images. Cluster LCBGs also have similar colors, sizes, and profiles (based on image concentration) to field LCBGs, but somewhat lower luminosities or surface-brightnesses (see Figures 1 and 3). Are these differences indicative of a more rapid pace of “down-sizing” in clusters, akin to the environmental dependence of star formation rates (Balogh et al. 1998, Martin et al. 2000)? To differentiate between an accelerated downsizing scenario and environmental effects of truncated or quenched star-formation (due, e.g., to stripping), star formation rates and mass measurements of cluster LCBGs are required. These measurements, well suited to PFIS, will tell us what differs between cluster and field: the mass-function of star-bursts, or the star-formation rate at a given mass?

## 3. LCBG Properties: The Future with PFIS

Preliminary results from our survey raise several important questions about cluster LCBG requiring spectroscopic measurements on 10m-class telescopes to answer. SALT’s Prime Focus Imaging Spectrograph (PFIS) will enable the needed high throughput, medium-resolution spectroscopy. Here we outline how we plan to use PFIS for our LCBG studies.

* Mass Uniformity: Our initial evidence indicates photometric differences in cluster and field LCBGs. We can characterize this difference more meaningfully by comparing field to cluster LCBG masses. In making this measurement for several clusters, we will map how the mass of star-bursting galaxies evolves with redshift and environment. Because of their small sizes, LCBG virial masses are measured by combining spatially-unresolved emission line-widths and HST or WIYN size measurements. The highest spectral resolutions of PFIS are ample for LCBG line-width measurements. By using PFIS’s bank of narrow band ($`\frac{\lambda }{\mathrm{\Delta }\lambda }50`$) filters, we can limit the range of the dispersed spectra on the detector to $``$75Å (rest-frame) around the \[OII\]$`\lambda `$3727 emission line for each cluster. This allows us to increase the spatial multiplexing of our multi-object masks (see the MMS concept described by Bershady et al., these proceedings) – and hence also our survey efficiency – by about a factor of 4. The minimum source spatial separation in the dispersion direction is $``$1 arcmin at the highest spectral resolutions (grating angle $`\alpha =50^{}`$), and decreases linearly with resolution ($`sin\alpha /cos\alpha `$, i.e., the Littrow condition for the VPH gratings). At the resolution used by Guzman et al. (1998), the minimum separation is 10 arcseconds.
* Star Formation Rates: If field and cluster LCBGs have similar size and mass functions, luminosity or surface-brightness differences are caused by lower star-formation rates – plausibly due to gas stripping in the dense clusters environments. We can measure these rates via low-resolution PFIS MOS observations that capture \[OII\]$`\lambda `$3727 and other blue nebular lines that will allow us also to determine metallicity and estimate extinction. Combined mass, size, shape, color, star-formation characteristics will allow us to carefully quantify physical differences between field and cluster LCBGs, and determine if they share similar histories of stellar processing of baryonic matter.
* Origins: Cluster Accretion? There is ample evidence in low redshift clusters that blue, star-forming galaxies are on the periphery of these systems, and just falling in. One way to test for an accretion origin of intermediate-redshift cluster LCBGs is to compare their velocity dispersion to the red cluster population (for which some measurements already exist, and will be augmented). LCBG velocity dispersions can be calculated efficiently from the MOS measurements described above, or via Fabry-Perot (FP) scanning with the lowest-resolution PFIS etalon. The advantage of the latter is the gain in spatial multiplex over the full 8 arcmin PFIS field, thereby providing complete kinematic characterization of the emission-line population in the cluster out to several cluster radii.

## 4. Summary

We have presented the initial results from the WLTV survey that will provide an extensive catalog of $`200`$ cluster and field LCBGs with high quality photometric information in areas around 10 clusters between $`0.3<z<1`$. With cluster LCBGs and over-densities identified from these photometric data, we plan to use PFIS to answer deeper questions about LCBGs: Are cluster and field LCBGs the same? Does the evolution of cluster and field LCBGs differ? What is the origin of cluster LCBGs? The answers to these questions tie in with parallel studies of other cluster and field populations, and will yield a cohesive picture of galaxy evolution over the past 3-7 Gyr.

### Acknowledgments.

This program has been funded by NSF AST-9970780 and AST-0307417, and NASA LTSA NAG5-6032 and STScI/AR-9917. Graduate student travel to this conference was funded by the NSF.

## References

Balogh, M. L., et al. 1998, ApjL, 504, L75

Barton, E. J. & van Zee, L. 2001, ApJl, 550, L35

Bershady, M. A., Jangren, A., & Conselice, C. J. 2000, AJ , 119, 2645

Cowie, L. L., et al. 1996 AJ, 112, 839

Dressler, A. 1980, ApJ , 236, 351

Guzmán, R., et al. 1997, ApJ , 489, 559

Guzmán, R., et al. 1998, ApJl, 495, L13

Hammer, F., et al. 2001, ApJ , 550, 570

Kobulnicky, H. A. & Zaritsky, D. 1999, ApJ , 511, 118

Koo, D. C. & Kron, R. G. 1988, ApJ , 325, 92

Koo, D. C., et al. 1995, ApJl, 440, L49

Lilly, S. et al. 1998, ApJ, 500, 75

Lowenthal, J. D. et al. 1997, ApJ , 481, 673

Mallén-Ornelas, et al. 1999, ApJl, 518, L83

Martin, C. L., Lotz, J., & Ferguson, H. C. 2000, ApJ , 543, 97

Phillips, A. C., et al. 1997, ApJ , 489, 543

Pisano, D. J., et al. 2001, AJ, 122, 1194

Steidel, C. C., et al. 1996, AJ, 112, 352

Trentham N. & Hudgkin, S. 2002, astro-ph 0202437

White, S. D. M. & Frenk, C. S. 1991, ApJ, 379, 52