text:
# Variational study of the Holstein polaron

## I Introduction

Ever since LandauLandau  suggested that the electron can be trapped by the deformable lattice, strongly coupled electron-phonon systems have been the subject of intensive examination. Besides the investigations of those systems in which the lattice is coupled to the whole electronic band, there has been significant interest in the physics of a single polaron, in which the electron and the associated lattice deformation form a quasi-particle, spatially and spectrally decoupled from the rest of the system. Lattice degrees of freedom make even a single polaron problem a many-body one. The analytical and numerical examinations of most electron-phonon models are thus difficult. For this reason, even the simple Holstein modelHolstein  (suggested in 1959) is still being investigated in recent works. Various methods have been proposed in order to calculate the polaron ground state of the Holstein model. Almost exact results (except for the adiabatic limit) have been obtained with the quantum Monte Carlo calculations,Raedt ; Kornilovitch  the global-local method,Brown  the density matrix renormalization group method,Zhang  and some exact diagonalization methods.Wellein1 ; Bonca

The main goal of this paper is to determine the elements important for the qualitative description of the polarons in the whole range of electron-phonon coupling. Some of them, although already known, are found to be better understood when supplemented with additional details.

The ground state of the Holstein system changes from the delocalized polaron state, in which the electron is nearly free, to the small and self-trapped polaron state, as the electron-phonon coupling $`g`$ increases. These two opposite limits are usually identified as the weak and strong coupling regime, respectively. In the weak coupling regime the influence of the small lattice deformation on polaron dynamics is very small, which makes the energies of polaron and electron hopping to neighboring sites similar, $`t_{pol}t`$.

The exact ground state is an eigenstate of the system momentum, regardless of the coupling.Gosar  Therefore the polaron ground state is delocalized for all parameters. However, in the strong coupling regime $`t_{pol}`$ becomes negligible, leading to self-trapped polaron states. According to Ref. Alexandrov1, , the dynamics of the small self-trapped polaron can be separated into two time scales. On the short time scale, the lattice deformation is centered at some lattice site, and the electron can virtually hop among the neighboring lattice sites. Only after a certain number of such events (of the order $`t_{pol}/t`$), the polaron as a whole tunnels to a new central site. The localized polaron states thus may lead to a very accurate estimation of the polaron ground state energy in the strong coupling regime. Accordingly, for the self-trapped polarons, an exact diagonalization method of calculating the localized polaron states, rather than the translationally invariant ones, can be used. By comparison to the results of other methods, it is shown that this approach introduces only minor errors in the ground state energy of the self-trapped polaron. Moreover, the localized polaron functions permit, unlike the translationally invariant ones, a separate analysis of the electron and phonon properties of the polaron. The local electron density, the mean lattice deformation, and the on-site zero point motions of the self-trapped polarons can be calculated in this way.

In the crossover (intermediate) regime, which is between the weak and the strong coupling regime, no known perturbation calculation converges, which complicates the discussion of the polaron nature. Although it has been proved that the change of the polaron ground state with $`g`$ is smooth,Gerlach  the physics of the rapid crossover between two opposite limits of the electron-phonon coupling, in the small interval of $`g`$’s, is not completely clear. In Ref. Bonca,  it has been claimed that the phonon excitation associated with the first excited state of the system is uncorrelated to the electron in the weak coupling regime, while it is confined to the electron in the strong coupling regime. The transition occurs in the crossover regime in which, in addition, the energy difference between the ground and excited states becomes small. In the present paper, the excited polaron states are treated by the new method which uses variational approach in order to define and solve the generalized eigenvalue problem. Even if this method does not converge systematically, it does provide new results concerning the nature of the polaron ground and first excited state.

## II General

The Holstein Hamiltonian reads

$`\widehat{H}`$ $`=`$ $`t{\displaystyle \underset{n,s}{}}c_{n,s}^{}(c_{n+1,s}+c_{n1,s})+\mathrm{}\omega {\displaystyle \underset{n}{}}b_n^{}b_n`$ (1)
$``$ $`g{\displaystyle \underset{n,s}{}}c_{n,s}^{}c_{n,s}(b_n^{}+b_n).`$
It describes the tight-binding electrons in the nearest-neighbor approximation, coupled to one branch of dispersionless optical phonons. $`c_{n,s}^{}`$ is the creation operator for the electron of spin $`s`$ at lattice site $`n`$, and $`b_n^{}`$ is the creation operator for the phonon. $`t`$ is the transfer (hopping) integral of the electron. $`\mathrm{}\omega `$ and $`g`$ are the phonon and the electron-phonon coupling energies, respectively. The Holstein Hamiltonian depends only on two ratios of relevant energy parameters: $`g/\mathrm{}\omega `$ and $`t/\mathrm{}\omega `$, i.e., the results will use $`\mathrm{}\omega `$ as the energy unit. It is often convenient to express the lattice vibrations in terms of the nuclei space and momentum coordinates,

$$\widehat{x}_n=x_0(b_n^{}+b_n),\widehat{p}_n=ip_0(b_n^{}b_n),$$
(2)
where $`x_0=\sqrt{\mathrm{}/2M\omega }`$ and $`p_0=\sqrt{M\mathrm{}\omega /2}`$ are space and momentum uncertainties of the harmonic oscillator ground state. $`M`$ is the mass of a nucleus and $`\kappa `$ denotes the spring constant, $`\omega ^2=\kappa /M`$. By noting that the electron-lattice displacement coupling constant $`\alpha `$, in $`g=\alpha x_0`$, is independent of $`M`$, an alternative set of Holstein Hamiltonian parameters can be introduced,

$$t=\mathrm{}^2/2m_{el}a^2,M=\frac{\kappa }{\omega ^2},\epsilon _p=\frac{\alpha ^2}{2\kappa }=\frac{g^2}{\mathrm{}\omega },$$

which is convenient for the discussion of the adiabatic regime $`Mm_{el}`$. Here, $`m_{el}`$ is the electron effective mass and $`a`$ is the lattice constant. It is worth noting that $`t`$ and $`\epsilon _p`$ (the binding energy of the polaron for $`t=0`$) are independent of $`M`$ and thus they are the only parameters relevant in the adiabatic limit.

By using standard conventions for $`c_{k,s}^{}`$ and $`b_q^{}`$,

$$c_{k,s}^{}=1/\sqrt{N}\underset{n}{}e^{ikna}c_{n,s}^{},b_q^{}=1/\sqrt{N}\underset{n}{}e^{iqna}b_n^{},$$

Eq. (1) can be rewritten in momentum space and separated in two mutually commuting parts, $`\widehat{H}=\widehat{H}_{k,q0}+\widehat{H}_{q=0}`$,

$`\widehat{H}_{k,q0}`$ $`=`$ $`2t{\displaystyle \underset{k,s}{}}\mathrm{cos}(ka)\widehat{n}_{k,s}+\mathrm{}\omega {\displaystyle \underset{q0}{}}\widehat{n}_q`$
$``$ $`g/\sqrt{N}{\displaystyle \underset{k,s,q0}{}}c_{k+q,s}^{}c_{k,s}(b_q+b_q^{}),`$
$`\widehat{H}_{q=0}`$ $`=`$ $`\mathrm{}\omega \widehat{n}_{q=0}gN_{el}/\sqrt{N}(b_{q=0}+b_{q=0}^{}),`$
where $`\widehat{n}_{k,s}=c_{k,s}^{}c_{k,s}`$, $`\widehat{n}_q=b_q^{}b_q`$, and $`N_{el}_{k,s}\widehat{n}_{k,s}`$.

The part which involves only the $`q=0`$ phonon mode, $`\widehat{H}_{q=0}`$, can easily be transformed into the diagonal form, by using the unitary operator $`\widehat{S}_q(\xi )`$,

$`\widehat{S}_q(\xi )`$ $`=`$ $`\mathrm{exp}(\xi b_q^{}\xi ^{}b_q),`$
$`\widehat{S}_{q=0}^1(\xi )\widehat{H}_{q=0}\widehat{S}_{q=0}(\xi )`$ $`=`$ $`\mathrm{}\omega \widehat{n}_{q=0}(N_{el}g)^2/N\mathrm{}\omega ,`$
where $`\xi =N_{el}g/\sqrt{N}\mathrm{}\omega `$. Any eigenstate of $`\widehat{H}_{q=0}`$, $`\widehat{S}_{q=0}(\xi )|n_{q=0}`$, has the same mean total lattice deformation, $`\overline{x}_{tot}=_n\overline{x}_n`$,

$`\overline{x}_{tot}`$ $`=`$ $`x_0\sqrt{N}n_{q=0}|\widehat{S}_{q=0}^1(\xi )(b_{q=0}^{}+b_{q=0})\widehat{S}_{q=0}(\xi )|n_{q=0}`$ (3)
$`=`$ $`x_0\sqrt{N}(\xi ^{}+\xi )=2x_0N_{el}g/\mathrm{}\omega ,`$
which is independent of $`t`$. Since the $`q=0`$ part of the Hamiltonian commutes with $`\widehat{H}_{k,q0}`$, it can be concluded that Eq. (3) is also valid for all eigenstates of the total Holstein Hamiltonian. The phonon part of these eigenstates can be represented in the form of a direct product of two groups of states, the first one includes $`q0`$ phonon modes, while the second includes only the $`q=0`$ phonon mode. This is useful because one can always check approximate computations by calculating $`\overline{x}_{tot}`$, or include this property in the computation itself. This specific property of the $`q=0`$ mode is not restricted to the particular dimension of the system, nor to the number of electrons. Moreover, it can also be found in some other models in which the electron-phonon coupling consists of the lattice deformation linearly coupled to the local electron density. A hint in this direction was reported for the first time in Ref. Feinberg, .

The total momentum of the system, $`\widehat{K}`$, is the sum of the electron and phonon momenta,

$$\widehat{K}=\underset{k}{}k\widehat{n}_k+\underset{q}{}q\widehat{n}_q,$$

and it commutes with the Hamiltonian. In the present treatment, only the low energy polaron states (the low lying states of the system for which the electron and lattice part of the wave function are spatially bound), are explicitly calculated. In this case, the total momentum $`K`$ of the system is also the polaron momentum.

## III Methods

The eigenstate computations reported here are based on the variational approach. Still, from the physical and mathematical point of view, there are significant differences among them. Each method is therefore described separately in the present section. It is important to notice that two qualitatively different kinds of polaron functions are employed, the localized ones and the translationally invariant ones. The primary objective of the methods with localized functions are self-trapped polaron states, i.e., the strong coupling regime. On the other hand, translationally invariant states are devised in the first place with the weak coupling and crossover regime in mind. The methods with a small number of variational parameters are meant to help in understanding the basic properties of the polarons in different regimes of the Holstein model. Again, the methods with a very large number of variational parameters are necessary to obtain accurate results for the polaron states.

### III.1 $`L`$ method: A simple localized polaron function

Let us start with a simple localized polaron wave function formed as a product of the electron and the lattice part, centered at the lattice site $`j`$,

$$|\phi _j=(\underset{n}{}\eta _nc_{j+n}^{})(\underset{m}{}S_{j+m}(\xi _m))|0.$$
(4)
Here, $`\eta _n`$ is the normalized electron function at site $`j+n`$, $`_n\eta _n^{}\eta _n=1`$, while $`S_{j+m}(\xi _m)`$ denotes the coherent state operator acting on site $`j+m`$, with a complex amplitude, $`\xi _m=\mathrm{}(\xi _m)+i\mathrm{}(\xi _m)`$,

$$S_{j+m}(\xi _m)=\mathrm{exp}(\xi _mb_{j+m}^{}\xi _m^{}b_{j+m}).$$

It is easy to see that operator $`S_{j+m}(\xi _m)`$ shifts the space and momentum coordinates of lattice vibration at a site $`j+m`$ by $`2\mathrm{}(\xi _m)x_0`$ and $`2\mathrm{}(\xi _m)p_0`$, respectively,

$`S_{j+m}(\xi _m)`$ $`=`$ $`e^{[i2\mathrm{}(\xi _m)p_0\widehat{x}_{j+m}/\mathrm{}]}e^{[i2\mathrm{}(\xi _m)x_0\widehat{p}_{j+m}/\mathrm{}]}`$
$`\times `$ $`e^{[i\mathrm{}(\xi _m)\mathrm{}(\xi _m)]}.`$
The variational energy of the state (4), $`\overline{E}_\phi `$, is independent of $`j`$, and is given by

$`\overline{E}_\phi `$ $`=`$ $`t{\displaystyle \underset{n}{}}\eta _n^{}(\eta _{n1}+\eta _{n+1})+\mathrm{}\omega {\displaystyle \underset{n}{}}\xi _n^2`$ (5)
$``$ $`g{\displaystyle \underset{n}{}}\eta _n^2(\xi _n+\xi _n^{}).`$
The minimization of the energy with respect to $`\xi _n`$ establishes a simple relationship between the lattice mean deformation and electron density $`\varrho _n=|\eta _n|^2`$,

$$\xi _n=(g/\mathrm{}\omega )\varrho _n\overline{x}_n=(\alpha /\kappa )\varrho _n,$$
(6)
so that only the equation for $`\eta _n`$ has to be solved. The well known approximate solution to this problem is the large Holstein polaronHolstein  valid in the long wave limit,

$$\eta _n=\frac{1}{\sqrt{2d_{pol}}}\mathrm{cosh}^1(n/d_{pol}),d_{pol}=\frac{2t\mathrm{}\omega }{g^2}=\frac{2t}{\epsilon _p}.$$
(7)
The numerical scheme suggested in Ref. Kalosakas, , and denoted here by $`L`$ ($`L`$ for localized), is not restricted to long waves, and is used here in order to obtain the exact minimum of Eq. (5). The energy, henceforth referred to as $`E_L`$, depends only on two relevant Hamiltonian parameters, $`t`$ and $`g^2/\mathrm{}\omega =\epsilon _p`$, and therefore $`E_L`$ and $`\overline{x}_m^L`$ are both independent of $`M`$. For this reason Eq. (6) is sometimes referred to as the adiabatic locking of electron and lattice coordinates.

It is noteworthy that the lattice part of the function $`|\phi _j^L`$ is a simple product of coherent states with real amplitudes. Consequently, the state of the lattice at some site is defined by the ground state of the displaced harmonic oscillator. The mean lattice deformation $`\overline{x}_m^L`$ corresponds to the equilibrium position of that oscillator, while the lattice zero point motion is approximated in Eq. (4) by that of the free lattice.

### III.2 $`T`$ and $`CT`$ methods: Translational polaron functions and the first excited state

Next, we shall study a translationally invariant solution composed of a linear superposition of the localized states (4),

$$|\mathrm{\Psi }_K=\frac{1}{\sqrt{N}_\mathrm{\Psi }}\underset{j}{}e^{iKja}|\phi _j.$$
(8)
$`|\mathrm{\Psi }_K`$ describes the polaron state with the momentum $`K`$. A similar type of function was first proposed by Toyozawa.Toyozawa  In the present work, $`_m\xi _m=g/\mathrm{}\omega `$ is used so that the mean total deformation of function (8) satisfies Eq. (3),

$$\overline{x}_{tot}=x_0\underset{m}{}\mathrm{\Psi }_K|(b_m^{}+b_m)|\mathrm{\Psi }_K=2x_0g/\mathrm{}\omega .$$

The expectation value of the polaron energy, $`\overline{E}_\mathrm{\Psi }`$, may be written in terms of $`|\phi _j`$,

$$\overline{E}_\mathrm{\Psi }=\frac{_\mathrm{\Delta }e^{iK\mathrm{\Delta }a}\phi _j|\widehat{H}|\phi _{j+\mathrm{\Delta }}}{_\mathrm{\Delta }e^{iK\mathrm{\Delta }a}\phi _j|\phi _{j+\mathrm{\Delta }}}=\frac{_\mathrm{\Delta }e^{iK\mathrm{\Delta }a}E_\mathrm{\Delta }}{_\mathrm{\Delta }e^{iK\mathrm{\Delta }a}S_\mathrm{\Delta }}.$$
(9)
$`E_\mathrm{\Delta }`$ and $`S_\mathrm{\Delta }`$ are given in the Appendix. A simple method for calculating the minimum of the energy $`\overline{E}_\mathrm{\Psi }`$ has not yet been proposed, but accurate results have been obtained in Ref. Zhao,  by using the Toyozawa method, which includes a very large number of variational parameters. Some additional approximations may be found in Refs. Feinberg,  and Shore, ; Hock, ; Venzl, ; Klamt, ; Magna, . The approximation used here simplifies the general expression in Eq. (8) by introducing the exponential form for functions $`\eta _n`$ and $`\xi _m`$,Venzl

$$\eta _n=CG^n,\xi _m=AB^me^{iKma},\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}<G,B<1.$$
(10)
Equation (10) defines a polaron function $`|\mathrm{\Psi }_K(G,B)`$ which is completely determined by two parameters, $`G`$ and $`B`$.

In what follows, two different approaches are presented. The first, denoted by the index $`T`$ ($`T`$ for translational), treats $`G`$ and $`B`$ as the variational parameters for which the energy minimum $`E_T`$ has to be found, and its corresponding polaron function is $`|\mathrm{\Psi }_K^T`$. The $`T`$ method gives good results in the weak and the strong coupling regime. Namely, in both of these limits, the function $`|\mathrm{\Psi }_K^T`$ becomes similar to the polaron function obtained by the appropriate perturbative calculations.

In the second approach, presented here for the first time, the variational method is used to define a generalized eigenvalue problem as follows. The polaron wave function, denoted by the index $`CT`$ ($`CT`$ for combination of translational functions), is rewritten as a linear combination of functions $`|\mathrm{\Psi }_K(G_n,B_n)`$,

$$|\mathrm{\Phi }_K^{CT}=\underset{n=1}{\overset{p}{}}a_n|\mathrm{\Psi }_K(G_n,B_n).$$
(11)
It is understood here that the functions $`|\mathrm{\Psi }(G_n,B_n)`$ form a set of $`p`$ generally nonorthogonal functions, defined by $`p`$ different pairs of parameters $`(G_n,B_n)`$. Again, the coefficients $`a_n`$ should be determined from the requirement that the expectation value of the energy,

$$\overline{E}_{CT}=\frac{\mathrm{\Phi }_K^{CT}|\widehat{H}|\mathrm{\Phi }_K^{CT}}{\mathrm{\Phi }_K^{CT}|\mathrm{\Phi }_K^{CT}},$$

is minimal,

$$\overline{E}_{CT}/a_n^{}=0,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}np,$$

or,

$`{\displaystyle \underset{n^{}}{}}\mathrm{\Psi }_K(G_n,B_n)|\widehat{H}|\mathrm{\Psi }_K(G_{n^{},}B_n^{})a_n^{}`$
$`=E_{CT}{\displaystyle \underset{n^{}}{}}\mathrm{\Psi }_K(G_n,B_n)|\mathrm{\Psi }_K(G_n^{},B_n^{})a_n^{}.`$
The solution of this generalized eigenvalue problem is a set of $`p`$ orthogonal polaron functions, $`|\mathrm{\Phi }_{K,m}^{CT}`$, with corresponding energies $`E_{CT}^{(m)}`$. The ground state energy is $`E_{CT}^{(0)}`$. One may always include the function $`|\mathrm{\Psi }_K^T`$ in the sum (11) in order to ensure that the energy $`E_{CT}^{(0)}`$ is the same or better than $`E_T`$, the energy computed by the $`T`$ method. Moreover, by paying some further attention to the starting set of functions $`|\mathrm{\Psi }_K(G_n,B_n)`$ in Eq. (11) at the outset, one is able to investigate the first excited state $`|\mathrm{\Phi }_{K,m=1}^{CT}`$ of the system, when this state corresponds to an excited polaron. The best results for the $`CT`$ method are obtained when the number $`p`$ of $`|\mathrm{\Psi }_K(G_n,B_n)`$ functions in Eq. (8) changes with the Hamiltonian parameters. The special case where the $`CT`$ method is used with constant $`p=2`$ is denoted by the index $`CT^{\underset{¯}{2}}`$.

### III.3 $`eT`$ and $`eL`$ methods: Exact translational and exact localized polaron functions

Finally, this paper presents the results of two numerical exact diagonalization methods.Ranninger ; Marsiglio ; Alexandrov ; Wellein1 ; Bonca  In order to compute the low energy polaron states, one approximates the infinite dimensional Hamiltonian matrix with a finite one, and proceeds with the exact diagonalization of this matrix. The lowest eigenvalue and eigenvector in such a reduced Hilbert space correspond to the polaron energy and wave function, respectively. For a large sparse matrix, the energy and the wave function can be calculated very accurately, by using an appropriate numerical scheme, in the present case the Lanczos algorithm.

The two exact diagonalization methods used here differ in the choice of the basis of the Hilbert space. In the first method,Bonca  denoted by $`eT`$ ($`e`$ stands for exact diagonalization and $`T`$ for translational), the general orthonormal state is given by

$`|n_0,n_1,n_1,\mathrm{},n_m,\mathrm{}_K^{eT}=`$
$`{\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{j}{}}e^{iKja}c_j^{}|n_0,n_1,n_1,\mathrm{},n_m,\mathrm{}_j,`$ (12)
which describes an eigenstate of the system with momentum $`K`$. $`n_m`$ is the number of phonons at the $`m`$th lattice site away from the electron. For example, at the site $`j`$, and at the nearest neighbor lattice sites left and right from it, there are $`n_0`$, $`n_1`$, and $`n_1`$ phonons, respectively. The Hamiltonian (1) does not mix states (12) with different momenta. Therefore, the polaron function obtained by the $`eT`$ method has the same $`K`$ momentum as the basis states. The current implementation of the $`eT`$ method is highly accurate, and from a practical point of view may be treated as exact.Bonca  For this reason, the $`eT`$ results can be used to determine the numerical errors present in the other methods.

The minimal number of states of the reduced basis necessary to obtain accurate results depends on the Hamiltonian parameters. This number for the $`eT`$ method increases very rapidly for $`\mathrm{}\omega g`$, $`t`$, which prevents its use for both large $`g`$ and $`t`$.

In the second exact diagonalization method, denoted by $`eL`$ ($`e`$ stands for exact diagonalization and $`L`$ for localized), the general orthonormal state of the chosen basis is more complicated than for the $`eT`$ method,

$`|i,n_0,n_1,n_1,\mathrm{},n_m,\mathrm{};\xi _m_j^{eL}=`$
$`c_{j+i}^{}[{\displaystyle \underset{m}{}}S_{j+m}(\xi _m)]n_0,n_1,n_1,\mathrm{},n_m,\mathrm{}_j.`$ (13)
Here the $`i`$ and $`m`$ indices are given with respect to the center of polaron, which is placed at the site $`j`$. Thus $`c_{j+i}^{}`$ creates an electron at the $`i`$th site from the polaron center at $`j`$, $`n_m`$ is the number of extra phonons at the $`m`$th site away of the polaron center, when the lattice is already distorted by the coherent state operators $`S_{j+m}(\xi _m)`$, i.e., $`\overline{x}_{j+m}=2x_0\xi _m`$. The $`eL`$ method calculates localized polaron wave functions. Namely, Eq. (13) describes a localized state, with the polaron center at site $`j`$ kept constant. The $`eL`$ method is therefore accurate only in the strong coupling regime, in which the effects of polaron delocalization are negligible (self-trapped polarons).

For a given set of Hamiltonian parameters, $`\xi _m`$ in Eq. (13) are determined by the use of the $`L`$ method, i.e., by minimizing the energy (5), $`\xi _m=\xi _m^L`$. If only those states (13) with all $`n_m=0`$ are used in calculations, the $`eL`$ method gives the same polaron wave function as the $`L`$ method. The additional states (13), with $`n_m`$ phonon excitations, are necessary to obtain the actual equilibrium positions of the lattice in the exact localized polaron state and the zero point motion of the renormalized lattice vibrations. It is worth noting that the electron and phonon parts of the $`eL`$ polaron wave function cannot be completely separated, as in the case of the $`L`$ function.

In the case of $`eL`$ method, the mean lattice deformation of the localized polaron is approximately taken care of by the product of the coherent states operators, $`_mS_{j+m}(\xi _m)`$, which keeps the necessary number of states (13), in the $`eL`$ method, relatively small. In order to reduce the basis of the Hilbert space, the maximal allowed distance of the electron and phonons from the polaron center has been limited here by the choice $`i,mD^{max}`$. The distance $`D^{max}`$ has been determined from the condition that $`\xi _m/\xi _0<10^4`$ if $`m>D^{max}`$. The maximal total number of phonons has been kept limited, retaining the states with $`_mn_m4`$. As the sum $`_mn_m`$ does not include the phonons associated with the coherent state operators $`S_{j+m}(\xi _m)`$, the small value of $`_mn_m`$ is not a restriction on the overall amplitude of the lattice deformation. The accuracy of the results obtained by the $`eL`$ method, supplemented by the two abovementioned criteria, depends of course on the values of parameters.

For the purpose of clarity, it seems appropriate at this point to review briefly the notation $`L`$, $`T`$, $`CT`$, $`eT`$, and $`eL`$ of all five presented methods. All methods based on the localized polaron function are denoted by the letter $`L`$ ($`L`$ and $`eL`$ methods), whereas all methods based on the translational function are denoted by the letter $`T`$ ($`T`$, $`CT`$, and $`eT`$ methods). The letter $`e`$ denotes an exact diagonalization method ($`eT`$ and $`eL`$ methods), while a single letter notation ($`L`$ and $`T`$) suggests the simplest form of the method.

## IV Results

As the variational methods of the preceding section themselves, the results of the corresponding calculations may be best understood in terms of the weak and strong coupling limits and the crossover between them. It has been shown previously,Romero3  on the basis of the Global-Local method for $`0.1\mathrm{}\omega <t<10\mathrm{}\omega `$, that the empirical relation

$$g_{ST}=\mathrm{}\omega +\sqrt{t\mathrm{}\omega }$$
(14)
describes well the values of parameters for which the variation of the effective polaron mass with $`g`$ is the fastest. It will be also argued here that $`g_{ST}`$ of Eq. (14) describes accurately the crossover from the weak to the strong coupling limit with respect to the nature of the $`K=0`$ ground state. The latter changes continuously from the light delocalized state in the weak coupling limit to the heavy self-trapped state in the strong coupling limit, with the anticrossing of the two states at $`gg_{ST}`$. The physical content of Eq. (14) is best understood by considering the limits of small and large $`t`$ with respect to $`\mathrm{}\omega `$, when, respectively, $`g_{ST}\mathrm{}\omega `$ and $`\epsilon _p^{ST}=g_{ST}^2/\mathrm{}\omega t`$. Both these conditions were qualitatively explained in Ref. Capone, . In the present section we first discuss the strong and the weak coupling limit, and then devote most of our attention to the crossover between the two.

### IV.1 Strong coupling limit

The nature of the self-trapped polaron may be discussed by analyzing the properties of the two polaron wave functions given by Eqs. (4) and (8). When the localized polaron functions $`|\phi _j`$ at different lattice sites are not orthogonal,

$$S_\mathrm{\Delta }=\phi _j|\phi _{j+\mathrm{\Delta }}\delta _{0,\mathrm{\Delta }},$$

the local properties and delocalization effects of $`|\mathrm{\Psi }_K`$ in Eq. (8) are a complex mixture. However, in both the weak and strong coupling regimes, the translational polaron function can be approximately written in terms of orthogonal localized polaron functions. In the weak coupling regime, the orthogonality follows from the electron part of the wave function,

$$S_\mathrm{\Delta }\underset{n}{}\eta _n^{}\eta _{n+\mathrm{\Delta }}\delta _{0,\mathrm{\Delta }},$$

while in the strong coupling regime, it follows from the lattice part,

$$S_\mathrm{\Delta }Y_\mathrm{\Delta }=\mathrm{exp}[\frac{1}{2}\underset{m}{}(\xi _m^{}\xi _{m+\mathrm{\Delta }})^2]\delta _{0,\mathrm{\Delta }}.$$
(15)
In Eq. (15), $`Y_\mathrm{\Delta }`$ is the Debye-Waller factor. The condition (15) corresponds to the regime of self-trapped polarons. The negligible contribution of the lattice part to the overlap of any two localized polaron functions at different lattice sites results then in a negligible polaron hopping energy $`t_{pol}`$.

One may notice that in the limit $`Y_{\mathrm{\Delta }0}0`$ the translational form of the polaron function (8) has no consequences on polaron energy, and the minimal values of the variational energies (5) and (9) coincide, i.e., $`E_\mathrm{\Delta }`$, $`S_\mathrm{\Delta }\delta _{0,\mathrm{\Delta }}`$ in Eq. (9). The hopping of the self-trapped polaron occurs only at the time scale which is much larger than the scale relevant for the local interplay between the electron and the lattice deformation. Therefore an accurate description of the local polaron properties may be obtained even if the hopping of the polaron is completely omitted by using only localized polaron functions. The translational invariance of the polaron may be, however, always restored in the same way as the function (4) is used to obtain (8).

The $`eL`$ method provides very accurate results for the ground state energy of the self-trapped polaron. Except that it neglects the polaron hopping, the local polaron function is calculated exactly. The error of the $`eL`$ method may be estimated by using the $`eT`$ method, that is, by subtracting the exact energy of zero momentum polaron state $`E_{eT}`$ from $`E_{eL}`$. For the electron-phonon coupling which is greater than the critical electron-phonon coupling in Eq. (14) just by $`\mathrm{}\omega `$, i.e., $`g=g_{ST}+\mathrm{}\omega `$, the maximal error of the $`eL`$ method is $`E_{eL}E_{eT}<3\times 10^4\mathrm{}\omega `$. This is shown in Fig. 1, in which two qualitative estimates of $`t_{pol}`$ in the strong coupling regime are plotted as well. The first estimate gives $`t_{pol}`$ as one fourth of the polaron bandwidth computed exactly by the $`eT`$ method,

$$t_{pol}\frac{1}{4}W^{eT}=\frac{1}{4}[E_{eT}(K=\pi /a)E_{eT}(K=0)],$$
(16)
while the second estimate, based on the $`L`$ method, multiplies the electron hopping energy with the small Debye-Waller factor,

$$t_{pol}tY_1^L=t\mathrm{exp}[\frac{1}{2}\underset{m}{}(\xi _m^L\xi _{m+1}^L)^2].$$
(17)
The estimation of the Debye-Waller factor in Eq. (17) is based on the evaluation of the lattice part of the overlap of the two neighboring $`|\phi _j^L`$ functions (the $`L`$ method gives good results for the mean lattice deformation of the self-trapped polarons). All three curves in Fig. 1 are similar functions of $`g`$, which is not surprising since all plotted quantities are related to the polaron hopping energy $`t_{pol}`$ in the strong coupling regime. Moreover, it may be seen from Fig. 1 that the error of the $`eL`$ method is almost equal to one half of the polaron bandwidth,

$$E_{eL}E_{eT}W^{eT}/2.$$

One of the advantages of the $`eL`$ method is that it permits separate calculations of the electron and the lattice properties, in spite of the fact that the electron and the phonon part of the $`eL`$ wave function are not separable. For instance, for the polaron centered at the origin, the associated mean lattice deformation $`\overline{x}_n`$ is given simply by the expectation value of $`\widehat{x}_n`$ of Eq. (2). Besides $`\overline{x}_n`$, in the present paper the mean real space uncertainty of the on-site lattice vibration $`\overline{\mathrm{\Delta }x}_n`$ and the product $`\overline{\mathrm{\Delta }x}_n\overline{\mathrm{\Delta }p}_n`$, where $`\overline{\mathrm{\Delta }p}_n`$ is the mean momentum uncertainty of the on-site lattice vibration, are calculated.

Figure 2 shows the data for two sets of parameter, corresponding to a small self-trapped polaron and a self-trapped polaron extended over few lattice sites, respectively. For the second set of parameter (large $`g`$ and $`t`$), the $`eT`$ method is not accurate enough, and the question arises of whether the $`eL`$ polaron state is really a good approximation of the system ground state. It is difficult to prove that a polaron state, with non-negligible hopping energy $`t_{pol}`$ and an energy close or lower than $`eL`$ energy, does not exist. This question is currently under investigation.

Differences between the results of the $`L`$ and the $`eL`$ methods are also analyzed. It is found that they are more pronounced for small $`g`$. We may see from the results in Fig. 2, that the mean lattice deformation differs between these two methods. In the case of the $`eL`$ method, the mean lattice deformation is more extended and the width of the polaron is slightly larger. Additionally, the electron density of the $`eL`$ method remains approximately proportional to the mean lattice deformation (as in Eq. (6), valid for the $`L`$ method), since

$$[x_n^{eL}\frac{2g}{\mathrm{}\omega }\varrho _n^{eL}x_0]/x_n^{eL}<1\%.$$

As has already been mentioned, the lattice part of the $`L`$ function describes a set of displaced, but unrenormalized, harmonic oscillators, so $`\mathrm{\Delta }x_n^L=x_0`$, $`\mathrm{\Delta }p_n^L=p_0`$, where $`x_0`$ and $`p_0`$ are defined in Eq. (2). One may notice from Fig. 2 that real space uncertainties of the lattice vibrations on the sites occupied by the polaron are larger in the $`eL`$ case than in the $`L`$ case, $`\mathrm{\Delta }x_n^{eL}>x_0`$, but uncertainties of on-site momentum lattice vibrations are smaller, $`\mathrm{\Delta }p_n^{eL}<p_0`$.

It may be concluded, from the first plot in Fig. 2 that the electron affects mostly the central site of the small polaron. The product of uncertainties for this site stays close to the free lattice value, $`\mathrm{\Delta }x_{n=0}^{eL}\mathrm{\Delta }p_{n=0}^{eL}x_0p_0`$. Therefore the phonon mode at the central site of the small polaron may be treated, in a good approximation, as harmonic. In particular, the renormalized frequency of this mode, $`\stackrel{~}{\omega }_{n=0}^{eL}`$, can be roughly estimated from the relation

$$\mathrm{\Delta }x_{n=0}^{eL}\sqrt{\mathrm{}/2M\stackrel{~}{\omega }_{n=0}^{eL}}.$$

Since $`\mathrm{\Delta }x_{n=0}^{eL}>x_0`$, it follows that $`\stackrel{~}{\omega }_{n=0}^{eL}<\omega `$. Therefore in the strong coupling regime the first excited state should correspond to the excitation of the renormalized phonon mode, rather than to the excitation of the phonon of energy $`\mathrm{}\omega `$, which is uncorrelated to the polaron. One may also notice that the energy of the mean lattice deformation is larger for the $`eL`$ method than the $`L`$ method (for the $`L`$ method this energy is minimal). This is compensated for, however, with the lower energy associated with the zero point motion of the $`\stackrel{~}{\omega }_{n=0}^{eL}`$ phonon mode, which makes the total polaron energy of the $`eL`$ method lower.

For the more extended self-trapped polarons, the renormalized normal phonon modes are expected to be spread over a number of lattice sites. Consequently, a number of different phonon modes contribute to the lattice displacement at the lattice sites occupied by the polaron. Thus the analysis of the on-site vibrations cannot give direct information on the renormalized lattice modes. Nevertheless, from the second plot in Fig. 2 one may notice that the product of uncertainties $`\mathrm{\Delta }x_n^{eL}\mathrm{\Delta }p_n^{eL}`$ shows a minor deviation from that of the harmonic vibration. This suggests that the renormalized lattice modes of the extended self-trapped polaron are harmonic as well in a good approximation.

### IV.2 Weak coupling regime

In the weak coupling regime, the $`T`$ method gives results close to the $`eT`$ results. Since the form of the function $`T`$ is quite simple, it will be used in this section as a basis for further discussion. In the weak coupling regime, the minimum of the variational energy $`E_T`$ corresponds to small values of the variational parameters $`G`$ and $`A`$. The standard perturbative ground state of the system with momentum $`K`$ in terms of the $`T`$ method polaron function may be written as follows:

$`|\mathrm{\Psi }_K^T={\displaystyle \frac{1}{\sqrt{N}_\mathrm{\Psi }}}{\displaystyle \underset{j}{}}e^{iKja}c_j^{}`$
$`\times (1+A{\displaystyle \underset{m}{}}B^me^{iKma}b_{j+m}^{})|0,K<K_c,`$ (18)
$`b_K^{}|\mathrm{\Psi }_{K=0}^T=b_K^{}{\displaystyle \frac{1}{\sqrt{N}_\mathrm{\Psi }}}{\displaystyle \underset{j}{}}c_j^{}`$
$`\times (1+A{\displaystyle \underset{m}{}}B^mb_{j+m}^{})|0,K>K_c.`$ (19)
For $`K<K_c`$ the wave function (18) has two parts. The main part corresponds to the free electron of momentum $`K`$, and the smaller part, proportional to $`A`$, corresponds to the electron dressed by one spatially correlated virtual phonon. At the threshold $`K_c`$, the energy of such a polaron state intersects with the energy of the system consisting of the zero momentum polaron and one extra phonon with momentum $`K_c`$, Eq. (19). So, for $`K>K_c`$ the ground state is achieved with one real phonon in the system which carries the system momentum and which is spatially uncorrelated with the polaron.Zhao ; Robin  For $`K<K_c`$ this state becomes the first excited state of the system. The difference between the energies of the ground state and the first excited state is the largest for $`K=0`$, and is equal to $`\mathrm{}\omega `$.

The validity of the perturbative treatment requires that the weight of the second term in Eqs. (18) and (19) is small, i.e., the mean number of phonons associated with the lattice deformation has to satisfy

$`\overline{N}_{ph}^{pol}`$ $`=`$ $`A^2(1+B^2)/(1B^2)`$ (20)
$`=`$ $`{\displaystyle \frac{g^2}{(\mathrm{}\omega )^2}}{\displaystyle \frac{(1B)(1+B^2)}{(1+B)^3}}1.`$
Here, $`A`$ has been eliminated by using Eq. (3). There are two ways to satisfy condition (20). Either the electron-phonon coupling is small, $`g\mathrm{}\omega `$, or the lattice deformation is spread to a large number of lattice sites, $`1B1`$. In the latter case, the total mean polaron deformation does not have to be small, $`\overline{x}_{tot}/x_0=2g/\mathrm{}\omega `$, since $`g`$ can be larger than $`\mathrm{}\omega `$.

The translationally invariant form of the wave function for the $`T`$ method, given by Eq. (8), provides an energy gain due to the polaron delocalization. At the same time, the spatial correlation between the electron and the lattice deformation has a finite length. For instance, in Eq. (8), this length is of the same order for $`|\mathrm{\Psi }_K`$ and $`|\phi _j`$. The perturbative calculation for $`B`$ in Eq. (18) gives

$$B=\mathrm{cos}(Ka)+\mathrm{}\omega /2t\sqrt{[\mathrm{cos}(Ka)+\mathrm{}\omega /2t]^21}.$$

$`B`$ measures the electron-lattice correlation length. $`B`$ is independent of $`g`$, which makes the correlation length finite, even in the limit $`g0`$ for which the lattice deformation vanishes, $`Ag/\mathrm{}\omega `$.

On the other hand, the electron-lattice deformation correlation length and the polaron delocalization range are of the same order for the localized functions. This may be easily seen from Eq. (4). In the limit $`g0`$ they both become infinite. Thus in the weak coupling regime one obtains a localized polaron state $`|\phi _j^L`$ of very large width, but associated with a tiny lattice deformation. This is specific for one dimensional systems in which an attractive symmetric potential has always a bound electron state. Therefore the corresponding polaron energy is less than the free electron energy $`2t`$. In higher dimensions, an arbitrary attractive symmetric potential does not have a bound electron state for sufficiently small g (or the electron binding energy is too small to balance the lattice deformation energy), and the total polaron energy is larger than $`2dt`$, where $`d`$ is the dimension of the system. This explains why, in the weak coupling regime, the adiabatic (localized) polaron functionsKalosakas ; Emin ; Kabanov  fail to have an energy lower than the free electron energy, for the dimension of the system greater than one. A detailed parallel perturbative analysis of polarons in one, two, and three dimensions is given in Ref. Romero6, .

### IV.3 Crossover regime

In the crossover regime the polaron hopping energy $`t_{pol}`$ is not negligible, and the translationally invariant polaron functions should be used in order to obtain the full physical picture of the polaron. In order to calculate the numerical errors of different methods accurately, the present discussion of the crossover regime is restricted to the values of $`t`$ smaller than $`25\mathrm{}\omega `$, for which results from the $`eT`$ method are available.

In order to examine the crossover regime it is instructive to calculate the energy difference between the ground and first excited state, $`E_{CT}^{(1)}E_{CT}^{(0)}`$. Let $`g_c`$ denote the value of $`g`$ for which this difference is minimal,

$$(E_{CT}^{(1)}E_{CT}^{(0)})/g|_{g=g_c}=0.$$

Our results indicate that $`g_c`$ is very close to the value of $`g_{ST}`$ given by Eq. (14). For example, for $`t=20\mathrm{}\omega `$, $`g_c=5.55\mathrm{}\omega `$, while $`g_{ST}=5.47\mathrm{}\omega `$. For smaller $`t`$, $`g_c`$ and $`g_{ST}`$ coincide even better. The analysis of the effective mass,Romero3  variational energy of the polaron ground state,Romero4  polaron size,Romero5  as well as the behavior of the first excitation energy, suggest that a dramatic change in the nature of the polaron ground state and the first excited state are intimately related in the crossover regime.

This can be well understood by considering the properties of the wave function in the $`T`$ method. Near $`g_c`$, this polaron wave function has two separate energy minima in the $`G`$-$`B`$ parameter space defined by Eq. (10), which become degenerated for $`g=g_c`$. Let the symbol $`<`$ denotes the lower minimum at $`g<g_c`$, and the symbol $`>`$ the lower minimum at $`g>g_c`$. $`|\mathrm{\Psi }^<`$ and $`|\mathrm{\Psi }^>`$ are the corresponding polaron wave functions. Even if they are not mutually orthogonal, they are still physically quite different. The numerical data show that the translational invariance of $`|\mathrm{\Psi }^<`$ contributes strongly to the polaron energy. On the other hand, the translational invariance of $`|\mathrm{\Psi }^>`$ has almost negligible energy contributions, i.e., $`|\mathrm{\Psi }^>`$ describes an almost self-trapped polaron. It is worth noting that the degenerate nature of the variational energy which has been reported for the Toyozawa methodZhao  is of the same kind as the one of the $`T`$ method discussed here. Namely, both $`T`$ and Toyozawa method are based on the same polaron function (8).

$`|\mathrm{\Psi }^<`$ and $`|\mathrm{\Psi }^>`$ can be combined to form new polaron functions,

$$|\mathrm{\Phi }^{CT^{\underset{¯}{2}}}=a_<|\mathrm{\Psi }^<+a_>|\mathrm{\Psi }^>.$$
(21)
It should be mentioned that a similar combination of two states has been already used in Ref. Cataudella,  in order to calculate the polaron ground state. $`|\mathrm{\Phi }^{CT^{\underset{¯}{2}}}`$ corresponds to Eq. (11) with $`p=2`$, which means that the $`CT^{\underset{¯}{2}}`$ method is implied. From this treatment the improved ground state $`|\mathrm{\Phi }_0^{CT^{\underset{¯}{2}}}`$ and the approximate first excited polaron state $`|\mathrm{\Phi }_1^{CT^{\underset{¯}{2}}}`$ are obtained. It may be seen from Fig. 3 that $`CT^{\underset{¯}{2}}`$ method describes an anticrossing of $`|\mathrm{\Psi }^>`$ and $`|\mathrm{\Psi }^<`$ states, which yields two orthogonal states, $`|\mathrm{\Phi }_0^{CT^{\underset{¯}{2}}}`$ and $`|\mathrm{\Phi }_1^{CT^{\underset{¯}{2}}}`$. In order to have a better illustration of that anticrossing in Fig. 3, the exact ground state energy $`E_{eT}`$ is subtracted from all plotted energy curves. For $`g<g_c`$ the light state $`|\mathrm{\Psi }^<`$ is lower in energy, and participates in the ground state more than the heavy state $`|\mathrm{\Psi }^>`$. However, as $`g`$ increases this balance changes continuously in favor of $`|\mathrm{\Psi }^>`$. The opposite trend is observed for the first excited state, which is heavier than the ground state for smaller $`g`$ and lighter for larger $`g`$.

The $`CT`$ method gives better results for the ground and first excited state when a large number of functions (large $`p`$) in Eq. (11) is used. From Fig. 3, one may estimate that $`E_{CT}^{(1)}E_{CT}^{(0)}\mathrm{}\omega /2`$ at $`g=g_c`$. Moreover, the energy $`E_{CT}^{(1)}`$, unlike $`E_{CT^{\underset{¯}{2}}}^{(1)}`$, satisfies $`E_{CT}^{(1)}<E^{(0)}+\mathrm{}\omega `$ for all $`g>g_c`$, which is an important improvement over the $`p=2`$ result.

## V Summary of results

Figure 4 shows the results for the ground state energy and the energy of the first excited state with total system momentum $`K=0`$, as functions of $`g`$, obtained by several different methods. Two plots correspond to $`t=5\mathrm{}\omega `$ and $`t=10\mathrm{}\omega `$, respectively. The ground state energy obtained by the $`eT`$ method is subtracted from all the other results. $`g_1`$ and $`g_2`$ are used to mark three different polaron regimes with respect to the strength of the electron-phonon coupling. For $`g<g_1`$ the mean number of phonons of the lattice deformation is less than one. Thus for $`gg_1`$ we recognize the weak coupling regime. For $`g>g_2=g_{ST}+\mathrm{}\omega `$ the $`eL`$ polaron energy has a negligible error (see Fig. 1), which means that for $`gg_2`$ polarons are self-trapped, and we recognize the strong coupling regime. The crossover regime is found in the interval $`g_1gg_2`$, with $`g_1<g_{ST}g_c<g_2`$.

It may be noted from Fig. 4 that in the weak coupling regime, the energy obtained by the $`L`$ method ($`E_L`$) is close to the free electron energy $`2t`$ (for $`t20\mathrm{}\omega `$, the absolute error of the $`L`$ method becomes greater than $`\mathrm{}\omega `$ in some intervals of $`g`$). In the strong coupling regime, for large $`g`$, $`E_L`$ approaches the exact polaron energy.

The error of the $`T`$ method, as may be seen from Fig. 4, is the largest in the crossover regime in which the results can be improved by a better choice of $`\eta _n`$ and $`\xi _m`$ in Eq. (10). In the strong coupling regime, the translationally invariant form of the function $`T`$ has no effect on the polaron energy, i.e., both the $`T`$ and $`L`$ methods give very similar results.

The energy of the first excited polaron state, $`E_{CT}^{(1)}`$, intersects the energy of the ground state plus one phonon, $`E_{eT}+\mathrm{}\omega `$, for $`g>g_1`$. After the minimum of $`E_{CT}^{(1)}E_{CT}^{(0)}`$ is reached in the crossover regime at $`g=g_cg_{ST}`$, $`E_{CT}^{(1)}`$ approaches $`E_{eT}+\mathrm{}\omega `$ asymptotically for $`g>g_2`$. As has already been pointed out in Sec. IV, the excitation of the renormalized phonon mode explains the nature of the first excited state in the strong coupling limit. In Ref. Kabanov,  perturbation theory was used to calculate the frequency of this local phonon mode, $`\stackrel{~}{\omega }^{}`$, to the lowest order in $`t\mathrm{}\omega /g^2`$,

$$\stackrel{~}{\omega }^{}=\omega \sqrt{1(t\mathrm{}\omega /g^2)^2}.$$
(22)
The obtained perturbative correction to the phonon frequency is adiabatic, i.e., the square root in Eq. (22) is independent of mass $`M`$. It is worth noting that for large $`g`$, $`E_{CT}^{(1)}`$ calculated here shows approximately the same behavior as the energy of the ground state with one additional phonon of frequency (22),

$$E_{CT}^{(1)}E_{eT}+\mathrm{}\stackrel{~}{\omega }^{}.$$

The lattice part of the $`CT`$ function is spatially symmetric with respect to the electron. Therefore in the strong coupling regime the first excited state, which has been identified here as a local renormalized phonon mode of the self-trapped polaron, should be basically a symmetric oscillation of the lattice deformation around the central polaron site. In the vicinity of $`g_{ST}`$, on the other hand, the nature of the first excited state is currently explained by the anticrossing of the self-trapped and the delocalized polaron state. The question of how the excited self-trapped polaron state in the strong coupling regime and the first excited state near $`g_{ST}`$ may be linked together is a matter of further considerations. One may speculate that the anticrossing of the excited self-trapped polaron state and the delocalized state will give the answer.

## VI Conclusions

The present paper discusses the ground and first excited states of the polaron for three different regimes of the electron-phonon coupling parameter $`g`$. The results can be briefly summarized as follows. In the strong coupling regime the polaron hopping energy to neighboring sites is negligible, and the self-trapped polaron states are obtained. The results of the $`eL`$ method suggest that the adiabatic picture of the localized polaron state is valid, in which some of the local lattice vibrations are renormalized by the presence of the electron. The numerical data show no significant deviation from the adiabatic locking relation (6) of the electron site density and the mean lattice deformation. In the small polaron case, the predominant effect of the electron is the lowering of the frequency of the vibration at the central polaron site. The excitation of the renormalized phonon mode corresponds to the first excited state of the small self-trapped polaron.

The nature of the polaron ground state in the crossover regime has been discussed in a number of papers, and its rapid change with $`g`$ has been well established numerically. The difference between the energies of the first excited state and the ground state, as a function of $`g`$, has a minimum for $`g=g_c`$. It is shown, by using the $`CT^{\underset{¯}{2}}`$ method, that the anticrossing of the self-trapped and the delocalized polaron state can link the behavior of the ground and first excited polaron state. According to the $`CT`$ method, for $`g>g_c`$ the effective mass of the ground state is larger than the effective mass of the first excited state, while for $`g<g_c`$ the opposite is true. In addition, it is found that $`g_c`$, which characterizes the first excited state, and $`g_{ST}`$, obtained from ground state analysis (Eq. (14)) almost coincide for $`t<25\mathrm{}\omega `$.

Upon further reduction of $`g`$ the total mean number of phonons bound by the polaron becomes smaller than one, and the weak coupling regime is reached. The nearly free electron is dressed by a cloud of virtual phonons, and its mass is slightly renormalized. The first excited state of the system, with momentum $`K0`$, can be viewed as the ground state of the polaron plus one additional uncorrelated phonon, rather than as an excited polaron.

Finally, it is worth noting that there is a simple sum rule for the mean total lattice deformation, Eq. (3), which is valid for any number of electrons and is independent of the system dimension. This important sum rule may be extended to some other models in which the lattice deformation is linearly coupled to the local electron density.

###### Acknowledgements.

The author would like to thank I. Batistić for fruitful collaboration during the course of this work.

*

## Appendix A

In order to calculate $`S_\mathrm{\Delta }`$ and $`E_\mathrm{\Delta }`$ of Eq. (9), the following expressions may be used:

$`Y_\mathrm{\Delta }`$ $`=`$ $`\mathrm{exp}[{\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}(\xi _m^{}\xi _{m+\mathrm{\Delta }})^2],`$
$`S_\mathrm{\Delta }`$ $`=`$ $`Y_\mathrm{\Delta }{\displaystyle \underset{n}{}}\eta _n^{}\eta _{n+\mathrm{\Delta }},`$
$`E_\mathrm{\Delta }`$ $`=`$ $`tY_\mathrm{\Delta }{\displaystyle \underset{n}{}}\eta _n^{}(\eta _{n+\mathrm{\Delta }+1}+\eta _{n+\mathrm{\Delta }1})`$
$`+`$ $`\mathrm{}\omega S_\mathrm{\Delta }{\displaystyle \underset{m}{}}\xi _m^{}\xi _{m+\mathrm{\Delta }}`$
$``$ $`gY_\mathrm{\Delta }{\displaystyle \underset{n}{}}\eta _n^{}\eta _{n+\mathrm{\Delta }}(\xi _n^{}+\xi _{n+\mathrm{\Delta }}).`$


format: text

text:
# On the Energy Required to Eject Processed Matter from Galaxies

## 1 Introduction

The observational evidence of powerful starbursts in dwarf galaxies has led to the idea that, due to their rather shallow gravitational potential, supernova (SN) products and, even the whole of the interstellar medium (ISM), may be easily ejected from the host dwarf systems, causing the contamination of the intra-cluster medium. This issue has been addressed in several papers. De Young & Gallagher (1990) and De Young & Heckman (1994) concluded that typically 1/6 of the shock-wave’s swept-up mass may escape its host galaxy, and coherent SN explosions may even totally destroy a galaxy’s ISM. More recently, however, Silich & Tenorio-Tagle (1998) and D‘Ercole & Brighenti (1999) have shown that within a dark matter (DM) and a low-density extended gaseous halo, a galactic wind does not easily develop, making the removal of the galaxy’s ISM very inefficient. Recently, Mac Low and Ferrara (1999, hereafter MF) have used several observational relations to define the galactic dark matter component by means of the central DM density ($`\rho _c`$), its maximum extent ($`R_h`$) and its characteristic scale ($`R_c`$). They made quantitative conclusions regarding the energetics required for the ejection of metals and of the ISM (”blow-away”) from galaxies with an ISM mass in the range (MISM = 106 \- 1010 M).

In none of the above mentioned studies, however, have the roles of the extragalactic environment nor of the rotation of the interstellar gas been examined in full detail. Meanwhile, Babul & Rees (1992) have indicated that the extragalactic medium pressure ($`P_{IGM}/k`$) may vary between 10-1 and 104 cm-3 K, making it one of the key parameters in the galactic-extragalactic medium connection.

In the present paper we reexamine the efficiency of starbursts in ejecting their newly produced metals, and in ejecting their whole ISM, from dwarf and normal galaxies, taking into account new boundary and initial conditions. These follow from the pressure equilibrium condition at the galactic-extragalactic boundary, and the relative ability of rotation and interstellar gas pressure to balance different fractions of the galaxy’s gravity. In particular, we show that in a self-consistent model, the ISM’s gas velocity dispersion and its density at the galaxy center are no longer free parameters, but are defined by the total ISM mass ($`M_{ISM}`$), the intergalactic pressure ($`P_{IGM}`$), and the fraction of the radial component of gravity balanced by pressure gradients. It follows from our considerations that the models of MF which have a low velocity dispersion and a flattened ISM disk-like density distribution, represent a limiting case in the more general scenario of starburst dwarf galaxy evolution.

Section 2 presents the derivation of our model galaxy. Section 3 is devoted to a search for the threshold energy input rate required to expel either metals or the ISM from galaxies with a total ISM mass in the range of 106 to 109 M. Section 4 is a discussion of our results.

## 2 The Galaxy Model

Our initial model (Silich & Tenorio-Tagle, 1998; hereafter referred to as Paper I) included stars, dark matter, and several isothermal interstellar gas components related to the ionized, neutral and molecular ISM phases. It was in fact an extension to the model developed by Morita (1982), Tomisaka & Ikeuchi (1988), Tomisaka & Bregman (1993) and Suchkov et al. (1994). To be consistent with the approach of MF we have simplified our model by using only the gravitational field imposed by the dark matter (DM) component and using only one of the ISM phases with an effective velocity dispersion CISM.

The spheroidal DM mass distribution can be approximated by

$$\rho _{DM}(r)=\frac{\rho _c}{1+(r/R_c)^2}.$$
(1)
The DM halo radius ($`R_h`$), characteristic scale height ($`R_c`$) and central density ($`\rho _c`$) are given by MF:

$`R_h=0.016\left({\displaystyle \frac{M_{DM}}{M_{}}}\right)^{1/3}H^{2/3}kpc,`$ (2)
$`R_c=0.89\times 10^5\left({\displaystyle \frac{M_{DM}}{M_{}}}\right)^{1/2}H^{1/2}kpc,`$ (3)
$`\rho _c=6.3\times 10^{10}\left({\displaystyle \frac{M_{DM}}{M_{}}}\right)^{1/3}H^{1/3}M_{}kpc^3,`$ (4)
where $`H`$ is the Hubble constant in 100 km s-1units. We adopt $`H`$=0.65 throughout the paper.

On the other hand, the interstellar gas density and pressure distributions are then defined by (see paper I)

$`\rho _{ISM}=\rho _0\mathrm{exp}\left[{\displaystyle \frac{3}{2}}\left({\displaystyle \frac{V_h}{C_{ISM}}}\right)^2\chi \right],`$ (5)
$`P_{ISM}={\displaystyle \frac{1}{3}}\rho _{ISM}C_{ISM}^2.`$ (6)
where $`\rho _0`$ is the ISM density at the galactic center, $`V_h=\sqrt{\frac{2GM_{DM}}{R_h}}`$ is the escape velocity at the DM boundary $`R_h`$, $`G`$ is the gravitational constant and $`M_{DM}`$ is the DM mass, which is given by

$`M_{DM}=4\pi {\displaystyle _0^{R_h}}\rho (\omega )\omega ^2d\omega =`$
$`4\pi \rho _cR_c^3\left[y_h\mathrm{arctan}(y_h)\right].`$ (7)
The function $`\chi `$ is

$$\chi =F(\omega )ϵ^2F(r)(1ϵ^2)F(0),$$
(8)
where $`r=\sqrt{(x^2+y^2)}`$ and $`\omega =\sqrt{(x^2+y^2+z^2)}`$ are the cylindrical and spherical radii, respectively, and the dimensionless radii $`y_\omega =\omega /R_c`$ and $`y_h=R_h/R_c`$. The factor $`(1ϵ^2)`$ is the fraction of the radial component of gravity that is balanced by pressure gradients (Tomisaka & Ikeuchi, 1988). The case $`ϵ=1`$ implies a full balance between the radial component of gravity and the circular rotation of the galaxy, whereas $`ϵ=0`$ implies spherically symmetric systems without rotation. The function F was defined in paper I and, in the case of a potential exerted by DM only, it becomes:

$`F(\omega )=1+{\displaystyle \frac{R_h}{R_c}}{\displaystyle \frac{4\pi \rho _cR_c^3}{M_{DM}}}\times [{\displaystyle \frac{1}{2}}(\mathrm{ln}(1+y_h^2)`$ (10)
$`\mathrm{ln}(1+y_\omega ^2))+{\displaystyle \frac{\mathrm{arctan}(y_h)}{y_h}}{\displaystyle \frac{\mathrm{arctan}(y_\omega )}{y_\omega }}],`$
$`for\omega R_h,`$
$`F(\omega )={\displaystyle \frac{R_h}{\omega }},for\omega >R_h,`$
The DM mass is related to the visible mass by (MF)

$$M_{DM}=3.47\times 10^8\left(\frac{M_{ISM}}{10^7M_{}}\right)^{0.71}M_{}.$$
(11)
One can also include a transition from a rotating disk to a non-rotating spherical halo by expressing the $`ϵ`$ in equation (8) as

$$ϵ=ϵ_c/\mathrm{exp}[(z/H_z)^2+(r/H_r)^2],$$
(12)
where $`ϵ_c`$ is the value of $`ϵ`$ at galactic center. This prevents the funnel-like density distribution that develops around the symmetry axis (see Suchkov at al. 1994). However, to be consistent with the considerations of MF, we have set the scale-lengths Hz and Hr to be much greater than or comparable to the DM halo radius.

Following MF we assume that the interstellar gas extends out to a cutoff radius $`R_{ISM}`$ given by

$$R_{ISM}=3\left(\frac{M_{ISM}}{10^7M_{}}\right)^{0.338}kpc.$$
(13)
However, we further assume that the galaxy is stable and that the ISM is in pressure balance with the intercluster gas at the galaxy outer boundary ($`P_{ISM}`$ = $`P_{IGM}`$). This defines the ISM gas density at the galaxy outer edge to be

$$\rho _G=\frac{3P_{IGM}}{C_{ISM}^2}.$$
(14)
Thus, given the total mass of the ISM ($`M_{ISM}`$) and the pressure of the intergalactic medium ($`P_{IGM}`$), only $`\rho _0`$ and $`C_{ISM}`$ must be found for fully defining the ISM mass distribution. This can be achieved by solving a system of two nonlinear equations

$`4\pi {\displaystyle _0^{R_h}}{\displaystyle _0^{R_h}}\rho _{ISM}(r,z)rdrdz=M_{ISM}`$ (15)
$`P_{ISM}(R_{ISM},0)=P_{IGM},`$ (16)
where $`\rho _{ISM}(r,z)`$ and $`P_{ISM}(r,z)`$ are defined by the equations (5) and (6). Note that $`\rho _{ISM}(r,z)`$ should be set equal to zero whenever its calculated value falls below $`\rho _G`$. The equations can be solved by an iterative method of high accuracy. Thus, in a self-consistent model, the ISM gas velocity dispersion and gas density at the galaxy’s center are no longer free parameters, but are determined by the total ISM mass $`M_{ISM}`$, the extragalactic pressure $`P_{IGM}`$, and the fraction of the radial component of gravity balanced by pressure gradients.

Figure 1 presents as an example, the ISM density distribution for $`M_{ISM}=10^8`$ M models with different fractions ($`ϵ^2`$) of the radial component of gravity balanced by the centrifugal forces. The figure shows the smooth transition from a flattened disk-like system dominated by fast rotation (Panel a) to a spherical, non-rotating galaxy (Panel c).

## 3 Energy requirements

Once a nuclear coeval starburst has all its stars on the main sequence, it begins to build a superbubble. In a constant density medium, the shock enters a phase of continuous deceleration from the very beginning of the superbubble evolution. However, if the starburst is placed in a plane with a stratified decreasing gas density, it will expand furthest in the direction of least resistance, deforming the originally round superbubble into an elongated remnant along the symmetry axis. In the case of an exponentially-decreasing density distribution, the remnant will then break out if the leading shock reaches a couple of galactic scale heights ($`H_e`$) while maintaining a supersonic speed, $`V_S`$ (Kompaneets 1960, Koo & McKee 1992), while in a Gaussian atmosphere the shock acceleration starts once it reaches a galactic scale-height $`H_g`$ (Koo & McKee 1990). As $`V_S=(P_{bubble}/\rho _Z)^{0.5}`$, where $`P_{bubble}`$ is the interior pressure, and $`\rho _Z`$ is the ISM density that decreases sharply along the $`Z`$ axis, the shock is forced to accelerate. Once in the acceleration phase, neither the external gas pressure nor the galaxy’s gravitational potential can restrain the shock, which ultimately will acquire speeds larger than the galactic escape velocity ($`V_{esc}`$). This occurs even if $`V_s<V_{esc}`$ prior to the acceleration phase that begins at $`Z`$ the gaseous scale-height. It is also well known that the swept up matter immediately behind the accelerating shock, also accelerates, producing the onset of Rayleigh-Taylor instabilities that fragment the swept-up matter shell. The remnant is then able to ”blowout”, driving its high pressure gas between the fragments into the halo of the host galaxy, where it once again will push the shock (see Tenorio-Tagle & Bodenheimer 1988, and references therein).

The minimum starburst energy input rate ($`L_{cr}`$), that leads to breakout, continuous acceleration, and expulsion of matter from galaxies into the IGM, either from a Gaussian or an exponential atmosphere was studied by Koo & McKee (1992). The authors used the threshold luminosity criterion ($`L_{breakout}`$) which assures that the shock will reach one scale height $`H`$ while maintaining a supersonic speed:

$$L_{breakout}=17.9\rho _0H^2C_{ISM}^3ergss^1$$
(17)
For an exponential atmosphere, Koo & McKee (1992) demand an $`L_{cr}3L_{breakout}`$. In a Gaussian density distribution the acceleration phase starts earlier, and thus we shall assume $`L_{cr}L_{breakout}`$. Note, that in the flattened disk-like configurations, the characteristic scale height $`H`$ may exceed the thickness of the galaxy’s disk $`Z_{ISM}`$, requiring that $`H`$ in the equation (17) be replaced by $`Z_{ISM}`$.

### 3.1 The Density Distribution along the Axis of Symmetry

Applying the breakout criterion to our model galaxies required accurate fits to the resultant density distributions (see Section 2) along the symmetry axis.

#### 3.1.1 Disk-like density distributions

In the case of flattened disk-like distributions, one can show that near the galaxy’s plane the resultant gas density distribution along the symmetry axis is fitted well by a Gaussian function:

$$\rho _{ISM}=\rho _0\mathrm{exp}[(Z/H_g)^2],$$
(18)
with the characteristic scale height

$$H_g=\sqrt{\frac{C_{ISM}^2}{2\pi G\rho _c}}.$$
(19)
Further away from the galaxy’s plane, the density distribution can be well approximated by an exponential function

$$\rho _{ISM}=\rho _G\mathrm{exp}[(Z_{ISM}Z)/H_e],$$
(20)
where $`Z_{ISM}`$ is the disk thickness along the Z axis. As the match to the density distribution requires two different functions, we have restricted the parameter $`H_e`$ by the condition that the density gradient, as given by the functions (18) and (20), must be smooth and continuous at the crossing point ($`Z_f`$). This condition can be represented by the following two equations:

$`log\rho _0\left({\displaystyle \frac{Z_f}{H_g}}\right)^2log(e)=`$
$`log\rho _G+{\displaystyle \frac{Z_{ISM}Z_f}{H_e}}log(e),`$ (21)
$`{\displaystyle \frac{2Z_f}{H_g^2}}log(e)={\displaystyle \frac{1}{H_e}}log(e).`$ (22)
Solving for He and Zf gives:

$`H_e={\displaystyle \frac{1}{2}}{\displaystyle \frac{Z_{ISM}^2}{Z_f}},`$ (24)
$`Z_f=Z_{ISM}\times `$
$`\left[1\sqrt{1\left({\displaystyle \frac{H_g}{Z_{ISM}}}\right)^2{\displaystyle \frac{log(\rho _0/\rho _G)}{log(e)}}}\right],`$
which uniquely define the intersection point $`Z_f`$ and the fitting parameter $`H_e`$ required in estimating the minimum energy input rate for breakout.

The ISM distribution along the Z axis for different disk-like galaxies and different extragalactic gas pressures (solid lines) and their analytic fit (dashed lines) are shown in the Figures 2 a and b.

#### 3.1.2 Spherical density distributions

The ISM density distribution for different galaxies under the assumption of $`ϵ`$ = 0, and different extragalactic gas pressures are shown in Figures 2 c and d. In these cases the density distribution far from the galaxy’s plane cannot be fitted by a Gaussian or an an exponential function. Thus to estimate the minimum energy input rate required to eject the metals from these galaxies we have only used the numerical scheme of paper I. In all cases we have assumed a constant mechanical energy input rate during the starburst supernova phase (4 $`\times 10^7`$ yr).

#### 3.1.3 The energy requirements

Our galaxy models are summarized in Table 1. The various models are labeled with several indexes that represent the logarithm of the ISM mass considered, $`ϵ`$ values that range between 0 and 0.9, and intergalactic gas pressure ($`P_{IGM}/k`$) values that are assumed to be either 1 or 100 cm-3 K. In the Table, column 1 gives the model identification. Columns 2 and 3 list the DM and ISM masses, column 4 gives the ISM cutoff radius, column 5 the intra-cluster gas pressure ($`P_{IGM}`$), and columns 6 and 7 show the interstellar gas velocity dispersion and escape velocity at the galaxy boundary $`R_{ISM}`$. The central gas number density is shown in column 8 and column 9 shows the ISM extent along Z axis ($`Z_{ISM}`$).

Figure 3 shows our energy estimates resultant from the numerical integration of the hydrodynamic equations for a coeval starburst with a constant energy input rate during the first 4 $`\times 10^7`$ yr of the evolution, and the initial density distributions shown in Figure 2 a - d. Figure 3 considers a range of galactic ISMs (from 106 to 109 M) and extreme values of $`ϵ`$ (= 0.9 for flattened disk-like density distributions and = 0 for spherical galaxies without rotation) and also for a range of values of $`P_{IGM}/k`$ (equal to 1 and 100 cm-3 K). Our results for disk-like gas distributions are in excellent agreement with the findings of MF, particularly for the case of $`P_{IGM}/k`$ = 1 cm-3 K. Note, however, that the threshold energy input rate in our models depends on the intergalactic gas pressure, $`P_{IGM}`$, and becomes rapidly larger for smaller $`M_{ISM}`$ and larger $`P_{IGM}`$.

In the case of the flattened disk-like configurations ($`ϵ`$ = 0.9; see Table 1), our numerical experiments indicate that the shock acceleration phase starts after 1.2 - 1.4 $`H_g`$, and thus it is the Gaussian rather than the exponential outer part of the disk-like density distributions the ones that define the onset of shock acceleration. Our results are in good agreement with the analytic criterion that ensures breakout (equation (17)) from the density distributions shown in Figure 2 a and b, particularly when the values of the Gaussian scale-height ($`H`$ = $`H_g`$) are used (filled symbols in Figure 3).

In Figure 3, values below each line imply total retention, while the region above each line indicates the expulsion of the hot superbubble interior gas (the new metals) out of disk-like distributions ($`ϵ`$ = 0.9), and of the new metals and the whole of the ISM in the spherical ($`ϵ`$ = 0) cases. Each line that separates the two regions marks the minimum energy input rate needed from a coeval starburst to reach the outer boundary of a given galaxy, regardless of the time that the remnant may require to do so. In the case of disk-like configurations, this energy input rate warrants breakout and the ensuing continuous acceleration. In most of these cases the ejection of the hot gas into the intergalactic medium occurs before the starburst supernova phase is over (4 $`\times 10^7`$ yr). In galaxies with $`ϵ`$ = 0, however, the evolution time considered largely exceeds the starburst supernova phase. This causes the decelerating shell of interstellar swept up matter to reach the galaxy edge by means of the momentum gathered during its early evolution. Figure 4 shows the results for a $`M_{ISM}`$ = 108 M galaxy under two extremes of intergalactic pressures ($`P_{IGM}/k`$ = 1 and 100 cm-3 K) as they become expelled by the lower limit mechanical energy input rate displayed in Figure 3. Throughout its evolution, the remnant in the high pressure galaxy undergoes a continuous deceleration, while the low pressure system is able to experience the acceleration promoted by breakout. In both cases however, the time needed to reach the galaxy outer boundary is much larger than the supernova phase and thus both remnants slow down by conservation of their own momentum, to reach the galaxy boundary and experience total disruption after 100 Myr and 200 Myr of evolution in the high and low intergalactic pressure cases, respectively.

The right-hand axis in our summary Figure 3 indicates the mass of the coeval starburst able to produce the energy input rate indicated in the left-hand axis. This has been scaled from the synthesis models of Leitherer & Heckman (1995) for a 106 M starburst under the assumption of a Salpeter IMF with a lower mass limit of 1 M, an upper stellar mass equal to 100 M, and a metalicity $`Z`$ = 0.1. In this case, the log of the mechanical energy input rate in erg s-1 is 40.5, and scales linearly with the starburst total mass.

In Figure 3, we have also indicated the results of another two calculations where we assumed a gravitational potential corresponding to a $`M_{ISM}`$ = 109 M galaxy. In these two models however, we have adopted a spherical halo with total masses of 108 M and 107 M and an intergalactic gas pressure equal to $`P_{IGM}/k`$ = 1 cm-3 K. The derived minimum energy input rate (or minimum mass of the coeval starburst required to expel the ISM) clearly indicates that it is the mass of the halo, instead of the disk component, what determines the minimum energy input rate necessary for ejection. Figure 3 shows that a halo with only 1/100 of the mass of the disk-like configuration ($`M_{ISM}`$ = 109 M) increases the value of the minimum energy input rate required for mass ejection by approximately an order of magnitude. An energy input rate nearly two orders of magnitude higher is needed if the spherical halo has a mass of 1/10 of the disk-like configuration.

## 4 Discussion

We have described a simple, self-consistent, approach to build the ISM gas distribution of galaxies, accounting for the gravitational potential exerted by a massive DM component, as well as for rotation of the ISM and the intergalactic gas pressure. The numerical experiments and analytical estimates have then showed that the final fate of the matter ejected from a starburst region, as well as that of the shocked ISM, is highly dependent on the boundary conditions. Following a continuous transition from a fast rotating, thin, disk-like to a spherical non-rotating interstellar gas distribution, we have found that the ISM to be more resistant to ejection than estimated in earlier papers.

We have shown (see Figure 3) that superbubbles evolving in galaxies that have a gaseous disk-like density distribution are likely to undergo the phenomenon of breakout. This allows them to accelerate and expel all of their newly produced metals, and perhaps even a small fraction of the interstellar medium, into intergalactic space. On the other hand, much larger energy input rates, or more massive coeval starbursts (up to 3 orders of magnitude larger), are required to provoke breakout or push a shell to the galaxy outer boundary for a spherically-symmetric ISM mass distribution. Even low mass ($``$ 1% - 10% of the total ISM mass), non-rotating subsystems increase the energy requirements by more than an order of magnitude. This makes the low density haloes, rather than DM itself, the key component in the evolution of dwarf galaxies.

The halo properties, on the other hand, are highly dependent on the gravitational potential. The presence of an extended DM component makes the halo mass distribution more smooth and extended. However, the halo parameters are also dependent on the properties of the intergalactic medium. To withstand a high intergalactic pressure, the halo has to support itself with a high random motion, which leads to a more homogeneous gas distribution compared to that found in low pressure surroundings.

Clearly the energy input rates derived here are lower limits to the amounts required for expelling matter from a galaxy. Particularly because only one component of the ISM was considered and because the central densities adopted are well below the values expected for the star forming cloud where the starburst originated. Our estimates thus neglect the effect of the starburst plowing into the parental cloud material. These are lower limits also because we adopted a constant energy input rate (see Strickland & Stevens 2000) and because our approach neglects an additional cooling by mass-loading process (Hartquist et al. 1986) and the presence of a magnetic field which also could inhibit expansion (Tomisaka 1998).

The indisputable presence of metals (in whatever abundance) in galaxies implies that the supernova products cannot be lost in all cases. Note in particular that many well known disk galaxies have a high metal abundance and a large number of centers of star formation. Most of these exciting star clusters are more massive than the $`10^4`$ M lower limit established by MF and the current paper as the minimum starburst mass required to cause mass ejection in the case of disk-like systems. This lower limit for disk-like galaxies with M$`{}_{ISM}{}^{}`$ 109 M (see Figure 3) implies that starbursts even smaller than the Orion cluster would break through the galaxy outer boundary and eject their supernova products into the intergalactic medium. Nevertheless, disk-like galaxies can avoid losing all their freshly produced metals by having a halo component, neglected in former studies, that acts as the barrier to the loss of the new metals.

The haloes, despite acting as the barrier to the loss of the new metal, have rather low densities ($`<n>10^3`$ cm-3) and thus have a long recombination time ($`t_{rec}=1/(\alpha n_{halo}`$); where $`\alpha `$ is the recombination coefficient) that can easily exceed the life time of the HII region ($`t_{HII}`$ = 107 yr) produced by the starburst. In such a case, the haloes may remain undetected at radio and optical frequencies (see Tenorio-Tagle et al. 1999), until large volumes are collected into the expanding supershells. Note that the continuous $`\mathrm{\Omega }`$ shape that supershells present in a number of galaxies, while remaining attached to the central starburst, and their small expansion velocity (comparable or smaller than the escape velocity of their host galaxy) imply that the mechanical energy of the star cluster is plowing into a continuous as yet undetected medium. Other authors (see MF 1999) have argued that blowout could leave fragments of the dense shell behind and that these observed with poor resolution may appear to form a continuous shell much smaller than the true extent of the shocked region. Clearly, X-ray observations of these supershells will help to decide if the hot processed material is enclosed or not by the expanding shell.

Thus the answer, the true limit for mass ejection from galaxies, must lie between the two extreme cases that we have investigated here. Note however, that in the presence of a halo, it is the mass of the halo that sets the limiting energy input rate required for mass ejection, and not the mass of the disk-like component. This argument applies to all galaxies whether spirals, amorphous irregulars or dwarfs.

We are grateful to our referee Dr. M-M. Mac Low for a speedy processing of our paper and the various comments that help to improve it. We also thank W. Wall for his careful reading of our manuscript. GTT also thanks CONACYT for the grant 211290-5-28501E which allowed for the completing of this study.