text:
# Metropolis simulations of Met-Enkephalin with solvent-accessible area parameterizations

## I Introduction

In nature biomolecules exist in the environment of solvents, thus the molecule-solvent interactions must be taken into account. It is very computer time consuming to simulate models for which the molecules of the surrounding water are treated explicitly. Therefore, a number of approximations of solvent effects have been developed. In the solvent-accessible area approach lr71 ; c74 ; em86  it is assumed that the protein-solvent interaction is given by the sum of the surface area of each atomic group times the atomic solvation parameter (ASP). The choice of a set of ASPs (also called hydrophobicity parameters or simply hydrophobicities) defines a model of solvation. However, there is no agreement on how to determine the universally best set of ASPs, or a at least the best set for some limited purpose. For instance, eight sets were reviewed and studied by Juffer et al. jeh95  and it was found that they give rather distinct contributions to the free energy of proteins folding.

In this paper we investigate how different ASP sets modify the Metropolis simulations of the small brain peptide Met-Enkephalin (Tyr-Gly-Gly-Phe-Met) at 300 K. The reason for the choice of Met-Enkephalin is that its vacuum properties define a reference system for testing numerical methods, e.g. ls87 ; okk92 ; ho93 ; mmm94 ; hoo99 ; b03 . Therefore, Met-Enkephalin appears to be well suited to set references for the inclusion of solvent effects as well, but we are only aware of few articles ls88 ; koh98 , which comment on the modifications due to including a solvent model.

We set our simulation temperature to 300 K, because room temperature is the physical temperature at which biological activity takes place. Most of the previous simulations of Met-Enkephalin in vacuum were performed at much lower temperatures or employed elaborate minimization techniques with the aim to determine the global energy minimum (GEM). Only recently b03  it was shown that the GEM is well accessible by local minimization of properly selected configurations from an equilibrium time series at 300 K. Precisely this should be the case for a GEM which is of relevance at physical temperatures.

For our simulations we use the program package SMMP smmp  (Simple Molecular Mechanics for Proteins) together with parallel tempering g91 ; hn96 ; ha97  (PT) and the recently introduced b03  biased Metropolis technique RM1 (rugged Metropolis – approximation 1). SMMP implements a number of all-atom energy functions, describing the intramolecular interactions, and nine ASP sets em86 ; oons ; sch1 ; sch2 ; jrf ; we92 ; sch4 ; bm  to model the molecule solvent interactions. We use the ECEPP/2 sns84  (Empirical Conformational Energy Program for Peptides) energy function with fully variable $`\omega `$ angles and simulate all nine ASP sets. For comparison we simulate also Met-Enkephalin in vacuum and with the distance dependent electrostatic permittivity $`ϵ(r)`$ of Ref. hrf85 .

The paper is organized as follows: The energy functions and Metropolis methods used are explained in Sec. II. In Sec. III we present our results from simulations of the brain peptide Met-Enkephalin. Summary and conclusions are given in Sec. IV.

## II Models and Methods

### II.1 ASP sets

In all-atom models of biomolecules the total conformational energy of the intramolecular interactions $`E_I`$ is given as the sum of the electrostatic, the Lennard-Jones (Van der Waals), the hydrogen bond, and the torsional contributions,

$`E_I=332{\displaystyle \underset{i<j}{}}{\displaystyle \frac{q_iq_j}{ϵr_{ij}}}+{\displaystyle \underset{i<j}{}}\left({\displaystyle \frac{A_{ij}^{\mathrm{LJ}}}{r_{ij}^{12}}}{\displaystyle \frac{B_{ij}^{\mathrm{LJ}}}{r_{ij}^6}}\right)+`$
$`{\displaystyle \underset{i<j}{}}\left({\displaystyle \frac{A_{ij}^{\mathrm{HB}}}{r_{ij}^{12}}}{\displaystyle \frac{B_{ij}^{\mathrm{HB}}}{r_{ij}^{10}}}\right)+{\displaystyle \underset{k}{}}U_k[1\pm \mathrm{cos}(n_k\varphi _k)].`$ (1)
Here $`r_{ij}`$ is the distance between atoms $`i`$ and $`j`$, $`q_i`$ and $`q_j`$ are the partial charges on the atoms $`i`$ and $`j`$, $`ϵ`$ is the electric permittivity of the environment, $`A_{ij}`$, $`B_{ij}`$, $`C_{ij}`$ and $`D_{ij}`$ are parameters that define the well depth and width for a given Lennard-Jones or hydrogen bond interaction, and $`\varphi _k`$ is the $`k`$th torsion angle. The units are as follows: distances are in Å, charges are in units of the electronic charge and energies are in kcal/mol.

One of the simplest ways to include interactions with water is to assume a distance dependent electrostatic permittivity according to the formula hrf85 ; ok94

$$ϵ(r)=D\frac{D2}{2}\left[(sr)^2+2sr+2\right]e^{sr}.$$
(2)
Empirical values for the parameters $`D`$ and $`s`$ are chosen, so that for large distances the permittivity takes the value of bulk water, $`ϵ=80`$, and the value $`ϵ=2`$ for short distances, i.e. for the interior of the molecule. Approximating solvation effects in this way is implemented as an option in SMMP. It allows to include solvation effects without any significant slowing down over the vacuum simulation with $`ϵ=2`$. The approach is clearly an oversimplification, because atoms which are close to each other do not necessarily have to be simultaneously in the interior of the molecule. Reversely, two atoms which are separated by a large distance may still be in the interior of the molecule. More elaborated approaches are asked for.

If the molecule-solvent interaction is proportional to the surface area of the atomic groups, it is given by the sum of contributions of a product of the surface area of each atomic group and the atomic solvation parameter em86 ,

$$E_{\mathrm{sol}}=\underset{i}{}\sigma _iA_i.$$
(3)
Here $`E_{\mathrm{sol}}`$ is the solvation energy and the sum is over all atomic groups. $`A_i`$ is the solvent accessible surface area and $`\sigma _i`$ the atomic solvation parameter of group $`i`$. The choice of a set of ASPs $`\sigma _i`$ defines a model of solvation. There are nine sets of ASPs in the SMMP package, which we list in Table 1. Columns one and two of this table gives the notations used in SMMP to identify the different sets.

Eisenberg and McLachlan em86  were the first to determine a set of ASPs (itysol$`=4`$, EM86 in the SMMP notation). For this, they considered the process of transferring atoms or groups of atoms from the interior of a protein to aqueous solution and used transfer energies of amino acids from $`n`$-octanol to water as reported in fp83 . The ASPs are then determined by least-square fitting. Octanol is chosen, because it apparently resembles the interior of a protein. With the exception of sch4  and bm , the other authors used similar methods with the major variation that instead of transfer energies with respect to octanol-water (o/w) also transfer energies with respect to vacuum-water (v/w) were used (for early determination of v/w transfer energies see c81 ; w81 ). The last column of Table 1 indicates whether the transfer energy is o/w or v/w. In the chronological order em86 ; sch1 ; we92  Eisenberg and collaborators contributed the parameter sets EM86, SCH1, WE92 and SCH3. Scheraga and collaborators oons ; jrf  contributed the parameter sets OONS and JRF. Here it should be noted that some of the original ASP sets were modified in course of time. For itysol$`=1,\mathrm{},8`$ SMMP implements the parameters as reviewed and tabulated in Ref. jeh95 , where in turn the sets SCH1 to SCH4 are simply taken from Schiffer et al.sch4 . Table 1 of SMMP smmp  lists the implemented ASPs for $`\mathrm{𝚒𝚝𝚢𝚜𝚘𝚕}=1,\mathrm{},8`$.

Somewhat special cases are the ASP sets SCH4 sch4  and BM bm . SCH4 was determined by comparison of the crystal structure from molecular dynamics simulations of small peptides and proteins in explicit water with similar simulations using an ASP solvation term (v/ws). The BM set of SMMP relies on a specific classification (cla) of atomic groups, where for all non-hydrogen atoms the solvation coefficients are set to 1 kcal/mol per Å2.

### II.2 Metropolis methods

For the updating of our systems we use PT with two processors, one running at 300 K and the other at 400 K. This builds on the experience b03  with vacuum simulations of Met-Enkephalin for which the following observations are made:

1. The integrated autocorrelation time $`\tau _{\mathrm{int}}`$ (defined below in this section) increases from 400 K to 300 K by a factor of ten for the (internal) energy and by factors of more than twenty for certain dihedral angles.
2. The energy probability densities (pds) at 300 K and 400 K overlap sufficiently, so that the PT method works, leading to an improvement factor of about 2.5 in the real time needed for the simulation (see Table I of Ref. b03 ).

A brief description of the PT algorithm is given in the following. PT performs $`n`$ canonical MC simulations at different $`\beta `$-values with Boltzmann weight factors

$$w_{B,i}(E^{(k)})=e^{\beta _iE_i^{(k)}}=e^H,i=0,\mathrm{},n1$$
(4)
where $`\beta _0<\beta _1<\mathrm{}<\beta _{n2}<\beta _{n1}`$ and a configuration is denoted by $`k`$. PT allows the exchange of neighboring $`\beta `$-values

$$\beta _{i1}\beta _i\mathrm{for}i=1,\mathrm{},n1.$$
(5)
These transitions lead to the change

$`\mathrm{\Delta }H`$ $`=`$ $`\beta _{i1}\left(E_i^{(k)}E_{i1}^{(k^{})}\right)\beta _i\left(E_{i1}^{(k^{})}E_i^{(k)}\right)`$ (6)
$`=`$ $`(\beta _i\beta _{i1})\left(E_i^{(k)}E_{i1}^{(k^{})}\right)`$
which is accepted or rejected according to the Metropolis algorithm, i.e. with probability one for $`\mathrm{\Delta }H0`$ and with probability $`\mathrm{exp}(\mathrm{\Delta }H)`$ for $`\mathrm{\Delta }H>0`$.

For the vacuum system the performance of the PT simulation was improved by an additional factor of two in Ref.b03  by using a first approximation, called RM1, to the rugged Metropolis scheme introduced there. A (short) simulation at 400 K was used to obtain estimates $`\overline{\rho }_j(v_j)`$, $`j=1,\mathrm{},24`$ of the pds of the 24 dihedral angles, which were then fed into the simulation. For a configuration change $`kk^{}`$ at temperature $`T_i`$ the new configuration accepted with the probability

$$p_{\mathrm{acpt}}=\mathrm{min}[1,\frac{\mathrm{exp}\left(\beta _iE_i^{(k^{})}\right)\underset{j=1}{\overset{24}{}}\overline{\rho }_j\left(v_j^{(k^{})}\right)}{\mathrm{exp}\left(\beta _iE_i^{(k)}\right)_{j=1}^{24}\overline{\rho }_j\left(v_j^{(k)}\right)}]$$
(7)
in the RM1 updating scheme. In the present paper we report the improvement due to this biased updating for some of the ASP sets.

In the vacuum simulation it is possible to determine the GEM by minimizing selected configurations of the 300 K time series. Here we apply the same procedure to our PT simulation of the ASP sets introduced in the previous subsection:

1. We determine the lower 10% quantile $`E_{0.1}`$ and the upper 10% quantile $`E_{0.9}`$ of the energy distribution of our time series. This is done by sorting all energies in increasing order and finding the values which cut out the lower and upper 10% of the data. For the statistical concepts see, e.g., Ref. brandt .
2. We partition the time series into bunches of configurations. A bunch contains the configurations from one crossing of the upper quantile $`E_{0.9}`$ to the next so that at least on crossing of the lower quantile $`E_{0.1}`$ is located between the two crossings of $`E_{0.9}`$. For each bunch we pick then its configuration of lowest energy. The idea behind this procedure is to pick minima of the time series, which are to a large degree statistically independent. In Fig. 1 the arrows indicate the energy values picked in that way from the first 600 configurations recorded in the RM1 simulation of Ref. b03 .
3. We run a conjugate gradient minimizer on all the selected configurations and thus obtain a set of configurations which are local energy minima. For the vacuum simulation b03  about 5% to 6% of the thus minimized configurations agree with the GEM.

To determine the speed at which the systems equilibrate, we measure the integrated autocorrelation time $`\tau _{\mathrm{int}}`$ for the energy and each dihedral angle. In particular the integrated autocorrelation times are directly proportional to the computer run times needed to achieve the same statistical accuracy for each system. They thus determine the relative performance of distinct algorithms. For an observable $`f`$ the autocorrelations are

$$C(t)=f_0f_tf^2$$
(8)
where $`t`$ labels the computer time. Defining $`c(t)=C(t)/C(0)`$, the time-dependent integrated autocorrelation time is given by

$$\tau _{\mathrm{int}}(t)=1+2\underset{t^{}=1}{\overset{t}{}}c(t^{}).$$
(9)
Formally the integrated autocorrelation time $`\tau _{\mathrm{int}}`$ is defined by $`\tau _{\mathrm{int}}=lim_t\mathrm{}\tau _{\mathrm{int}}(t)`$. Numerically, however, this limit cannot be reached as the noise of the estimator increases faster than the signal. Nevertheless, one can calculate reliable estimates by reaching a window of $`t`$ values for which $`\tau _{\mathrm{int}}(t)`$ becomes flat, while its error bars are still reasonably small. This is the method we employ in the next section, see Ref. sokal  for a more detailed discussion of the integrated autocorrelation time.

## III Results

### III.1 Autocorrelations

The PT simulations with temperatures $`T_0=400`$K and $`T_1=300`$K are performed on the system in vacuum ($`ϵ=2`$), with $`ϵ(r)`$ given by Eq. (2) and for the nine ASP sets of Table 1. The dihedral angles updated in our simulations are fully variable in the range from $`\pi `$ to $`\pi `$. We keep a time series of $`2^{16}=65536`$ configurations for each replica (i.e., each of the two processor), in which subsequent configurations are separated by 32 sweeps. A sweep is defined by updating each dihedral angle once sequentially. Before starting with measurements $`2^{18}=262144`$ sweeps are performed for reaching equilibrium. Thus, the entire simulation at one temperature relies on $`2^{21}+2^{18}=2,359,296`$ sweeps. On the Cray T3E, this takes about 14 hours for the vacuum system and $`5\times 14`$ hours for each ASP set.

Results of the average energy, acceptance rates and integrated autocorrelations times for the energy $`E`$ variable are shown in Table 2. For the vacuum simulations and the ASP sets OONS and EM86 the time-dependent integrated autocorrelations times (9) are shown in Fig. 2. In each case a window of $`t`$ values is reached for which $`\tau _{\mathrm{int}}(t)`$ does no longer increase within its statistical errors. In the case of the vacuum simulations it even does decrease, but this is not significant due to the statistical error. These windows are then used to estimate the asymptotic $`\tau _{\mathrm{int}}`$ values of Table 2. With the exception of the ASP sets JRF and BM, the integrated autocorrelations times of all other sets are determined in the same way.

From Table 2 we see that the acceptance rates of the solvent models JRF and BM are much lower than for the other models. In essence the simulations of these two models get stuck, which implies that their integrated autocorrelation times cannot be measured. This is illustrated in Fig. 3 for the time-dependent integrated autocorrelation time of the energy at 300 K. The function $`\tau _{\mathrm{int}}(t)`$ increases rapidly until it gets lost in the noise. The pds of the dihedral angles of these two models are also erratic and the conclusion is that they cannot be used to describe Met-Enkephalin in solvent.

The energy couples to all dihedral angles and its integrated autocorrelation time is characteristic for the entire system, while the integrated autocorrelation times of the single dihedral angles vary heavily from angle to angle. For all our systems, but JRF and BM, we show in Fig. 4 the integrated autocorrelation times at 300 K for the energy and all dihedral angles. The notation $`v_i`$, $`i=0,1,\mathrm{},24`$ is used, where $`v_0`$ is stands in for the energy $`E`$ and for $`i=1,\mathrm{},24`$ the $`v_i`$ are the dihedral angles as used in the SMMP computer program. The relationship of the $`v_i`$ angles to the conventional notation for dihedral angles and their residues is summarized in Table 3, where it should be noted that the SMMP notation smmp  differs from other literature ls87 ; mmm94 .

In Fig. 4 we see that for each dihedral angle $`v_i`$ the integrated autocorrelation times $`\tau _{\mathrm{int}}[v_i]`$ for the three solvent models OONS, WE92 and SCH2 are smaller than for the remaining systems, including the vacuum system. For the integrated autocorrelation time of the energy $`\tau _{\mathrm{int}}[E]`$ this observation is already obvious from Table 2. In particular, this means that the OONS, WE92 and SCH2 models require far less statistics than the vacuum run for achieving the same accuracy of results. Using the $`\tau _{\mathrm{int}}[E]`$ results of Table 2, we find a factor in the range 7 to 10, which more than offsets the factor of 5 by which the ASP model simulations are slower than the vacuum simulation. In the following the solvation models OONS, WE92 and SCH2 define the “fast class”, while the other models shown in Fig. 4 constitute the “slow class” (the models JRF and BM are omitted from this classification). “Good” behavior of the models OONS and WE92 has been previously observed ok98 .

The autocorrelation times in the fast class are so small that the resolution of 32 sweeps in our recorded time series becomes too crude. Namely, autocorrelations over less than 32 sweeps are then not measured and the integrated autocorrelation time approaches one as soon as autocorrelations stay within the range of 32 sweeps. To investigate this point further, we performed for the OONS, WE92 and SCH2 models simulations for which the configurations were recorded every four sweeps and the total statistics was reduced by the factor 1/8. In the new units of four sweeps the integrated autocorrelation time is larger by a factor which is bounded by $`8=32/4`$. The bound is assumed, if there is no improvement due to integrating additional small fluctuations out (i.e. due to the additional configuration in between the 32 sweeps, which are now kept in the time series).

For WE92 we report in the PT column of Table 3 the integrated autocorrelation times from the simulation with reduced statistics. For many dihedral angles the increase lies well below the factor of eight, showing that we gain in accuracy by averaging over small fluctuation within the range of 32 sweeps. On the other hand, nothing is gained by this extra averaging for several angles with large autocorrelations. In those case the simulations yield within their statistical errors the upper bound 8.

To supplement the vacuum results of Ref. b03 , we repeated the WE92 PT simulations by using estimates of the dihedral pds from 400 K as input for the biased updating of Eq. (7). These results are reported in the PT-RM1 column of Table 3. As in the case of the vacuum simulations, we find an improvement of the PT performance by a factor of approximately two, which is also obtained for the other models of the fast class. For the slow class we checked by direct improvement of the original simulations on the ASP models EM86 and SCH4 and find again an acceleration by a factor of about two when we are using RM1 updating.

### III.2 Structure

For all our simulations we applied the method outlined in subsection II.2 to determine local energy minima and some results are summarized in Table 4: $`E_{0.1}`$ and $`E_{0.9}`$ are the lower and upper 10% quantiles of the energy and $`N_{\mathrm{conf}}`$ denotes the number of minima of the time series prepared for further minimization. The lowest energy found in this minimization process is denoted by $`E_{\mathrm{min}}`$ and $`N_{\mathrm{hits}}`$ is the number of times the lowest energy configuration was hit. While the absolute values of $`E_{0.1}`$ and $`E_{0.9}`$ vary considerably from set to set, the differences $`E_{0.9}E_{0.1}`$ stay similar. The explanation is that the ASP sets differ by large additive constants to the energy.

Again, the results of the JRF and BM solvent models are erratic. The BM model is entirely frozen, $`N_{\mathrm{conf}}=2`$ at 400 K and $`N_{\mathrm{conf}}=1`$ at 300 K. Therefore, we do not give minimization results for BM. For JRF the $`N_{\mathrm{conf}}`$ numbers are more reasonable, but still by a factor of one third and less smaller than the $`N_{\mathrm{conf}}`$ numbers of each other system. JRF is also disregarded in the following discussion.

Only if $`N_{\mathrm{hits}}>1`$ holds we have an indication that we found the GEM. Interestingly, this happens for none of the ASP solvent models, while it is the case for the vacuum and the $`ϵ(r)`$ simulations (notably already at 400 K). Quite some time ago Li and Scheraga ls87 ; ls88  developed a Monte Carlo minimization method and applied it to Met-Enkephalin in vacuum and in solvent modeled by OONS. While for the vacuum system their method converged consistently to the GEM, all their five runs of the solvent model led to different conformations with comparable energies. They interpreted their results in the sense that Met-Enkephalin in water at $`20^{}`$C is likely in an unfolded state for which a large ensemble of distinct conformations co-exist in equilibrium. A consistent scenario was later observed in NMR experiments grcahi .

Although the minimization method of Li and Scheraga is entirely different from ours, they essentially tested for valleys of attraction to the GEM at room temperature, quite as we do in the present paper. So, we have not only confirmed their old result, but find that it is also common to a large set of ASP models implemented in SMMP. Neither the method by which an ASP set was derived, nor whether it belongs to the fast or slow class, appears to matter with this respect.

To give one example, the frequency of local energy minima of the WE92 solvation model as obtained by our minimization procedure from the 300 K time series is depicted in Fig. 5. $`N_{\mathrm{conf}}=2453`$ minimizations are performed. Our lowest energy state is only found once and the same holds for the close-by low energy states. Fig. 5 should be compared with Fig. 2 of Ref. b03 , where the frequency of the low energy minima of the vacuum simulation is shown. There the lowest energy state relies on 107 entries out of 1913 minimizations errata .

In a search for structural differences of Met-Enkephalin in vacuum, or in the $`ϵ(r)`$ system, versus the ASP models, we looked at the pds of the dihedral angles at 300 K. For all systems together there are $`9\times 24=216`$ figures to consider. At the first look the pds of the different systems are amazingly similar, independently of whether they are from systems of the fast or slow class, from an ASP model, from the vacuum or from the $`ϵ(r)`$ simulation. A more careful investigation reveals differences, which appear to relate to the distinct behavior under our minimization. For the dihedral angle $`v_7`$ this is illustrated in Fig. 6 and Fig. 7. Its probability densities are compared at 300 K and 400 K. For the vacuum simulation the pds are depicted in Fig. 6 and from 400 K to 300 K we observe an increase of the peak which is located close to the arrow, which indicates the vacuum GEM value of $`v_7`$. In contrast to this, the wrong peak increases in Fig. 7, where the pds are shown for the WE92 solvent model.

One may suspect that the difference between the models of our fast and slow class is simply due to an effectively higher temperature for the three models of the fast class. To gain insight into this question, we calculate the entropies of our pds. Each pd is discretized as a histogram of 200 entries, $`\rho _{ij}`$, where $`i=1,\mathrm{},24`$ labels the dihedral angles according to Table 3 and $`_{j=1}^{200}\rho _{ij}=1`$. The entropy of the pd of a dihedral is then defined by

$$S_i=\underset{j=1}{\overset{200}{}}\rho _{ij}\mathrm{ln}\rho _{ij}$$
(10)
and the total entropy of the pds of an APS model is $`S=_iS_i`$. In Fig. 8 the thus obtained entropies are depicted for all our models. The lines between the data points have no other meaning but to guide the eyes. The dips for the JRF and the BM model show, again, that their configurations are essentially frozen. For the other we see a decrease of entropy from 400 K to 300 K, but we find no larger entropy for the models of the fast class than for the models of the slow class. Therefore, the effective temperature scenario is rules out. Instead, it seems that for the models of the fast class the solvent has some kind of “lubrication” effect, which accelerates the simulation.

Strong similarities between the ASP models of the fast class on one side and the ASP models of the slow class on the other side are found for the solvation energies, the gyration radii and the end to end distances.

The solvation energies $`E_{\mathrm{sol}}`$ (3) measured during our simulations of ASP models are solvent-vacuum (s/v) transfer energies. There are structural difference between the typical configurations of an ASP model time series and the vacuum time series. Consequently, the average s/v transfer energies are not identical with the average vacuum-solvent (v/s) transfer energies, which are obtained by calculating $`E_{\mathrm{sol}}`$ of the solvent models on the configurations of the vacuum time series. The s/v as well as the $``$v/s average transfer energies are collected in Table 5. The averages are taken for the canonical time series at 300 K and at 400 K. At 300 K averages are also taken for the time series minima (indicated by arrows in Fig. 1) and for the local minima (which are obtained by running the conjugate gradient minimizer on the time series minima). For the gyration radii $`R_{\mathrm{gy}}`$ and the end to end distances $`R_{\mathrm{e}\mathrm{e}}`$ the same averages are collected in Table 6 (definitions and software are given in SMMP).

For the transfer energies the over-all effect is hydrophilic for the APS models of the fast class and hydrophobic for the APS models of the slow class. Within each class the values are quite similar, despite differences in the interaction coefficients (see Table 1 of Ref.smmp ). As expected the over-all transfer energies of the JRF and BM models are out of the reasonable range, JRF to the hydrophilic and BM to the hydrophobic side. Our Table 6 shows that we observe the extended structures found in previous simulations ls88 ; koh98  and in NMR experiments grcahi  only for the APS models of the fast class.

## IV Summary and Conclusions

We have performed Met-Enkephalin simulations at room temperature (300 K) for the solvation models of Table 1. Quantitative results obtained in that way cannot be trusted, apparently because the methods to derive the ASPs have been quite crude. Also our simulations do not give information that would allow us to pick a best ASP set for the intended purpose of simulating Met-Enkephalin at 300 K. Nevertheless, we obtain a qualitative overview of a number of interesting consequences, which one can expect by including solvation effects via an ASP model in Metropolis calculations.

Two of the ASP sets (JRF jrf  and BM bm  as implemented in SMMP smmp ) suffer from so large autocorrelations that for them Metropolis simulations at 300 K are in essence impossible. Their dihedral angles are essentially frozen. These two ASP sets are certainly erratic, as 300 K is a temperature at which thermodynamic fluctuations of the systems are expected (also these two sets perform badly at 400 K).

The remaining nine models, seven ASP sets, $`ϵ=2`$ vacuum, and an $`ϵ(r)`$ system hrf85 , fall into a fast and a slow class with respect to their integrated autocorrelation times, see Fig. 4. Vacuum simulations are in the slow class. This leads to the interesting feature that it takes less computer time to estimate physical observables at room temperature in the fast solvation models OONS oons , WE92 we92 , and SCH2 sch2 ; sch4 , than it takes for vacuum, despite the substantial increase of the computer time per sweep by a factor of about 5 for the solvation models over the vacuum system. We have no clear clue why some models have a fast and others a slow dynamics. To derive the parameters of OONS and WE92 vacuum-water (v/w) transfer energies were used, but for SCH2 it was octanol-water (o/w). Also the slow class features v/w as well as o/w ASP models.

We applied the minimization procedure of Ref. b03  in an attempt to locate the GEM for the nine systems, which are reasonably well-behaved under Metropolis simulations at 300 K. The GEM is unambiguously found for the vacuum system and for the simulation with a distance dependent electrostatic permittivity. No true GEM is found for any of the remaining seven ASP models. This confirms an old result of Li and Scheraga ls88 , who concluded that at room temperature Met-Enkephalin in water is likely in an unfolded state. To get a better understanding of this result, we studied at 300 K the dihedral pds in some details. At a first glance they look quite similar for all the models in the fast as well as in the slow class. Differences are found for a number of details, which may allow to explain why the 300 K configurations of the ASP models behave entirely different under our minimization procedure than the vacuum and the $`ϵ(r)`$ systems.

The central question, which remains to be settled, is whether ASP models will ultimately allow for accurate Metropolis simulations of biomolecules like Met-Enkephalin in solvent or not. In principle, this could be decided by determining whether ASPs exist which reproduce accurately mean energies of explicit solvent simulation around a large number of fixed Met-Enkephalin configurations.

###### Acknowledgements.

We would like to thank Prof. P. Grassberger for useful discussions and generous support of this work. The computer simulations were carried out on the Cray T3E of the John von Neumann Institute for Computing. BB acknowledges partial support by the U.S. Department of Energy under contract No. DE-FG02-97ER41022. In the final stage of this work we became aware of two articles where ASP models are studied for the helix-coil transition of polyalanine PeHa03 .


format: text

text:
# A correlation between the equilibrium and transport properties of intercalation systems.

## I Introduction

Intercalation processes find their application in many technologically important domains, such as the design of hydrogen-storage systems, rechargeable high-energy batteries, electrochromic devices, (see Ref. for a review), and superconductors .

The insertion process can be viewed (at least, in some aspects) as an adsorption of guest particles on the host lattice. For charged particles the ionic charge inside the matrix is compensated by the electrons. For that reason the electrochemical intercalation is also similar to a three-dimensional adsorption of neutral species. Based on this analogy the intercalation is traditionally described within the lattice gas (LG) model. In this approach all the properties (intercalation isotherm or capacity-concentration dependence) are connected with an ordering the guest on different adsorption sites of a rigid host lattice - the configurational transitions.

However, the host response to the accommodation of the guest species is not negligible. For instance, the insertion induces the stress into the host matrix. This may lead to segregation effects  or even to instabilities of the host-guest system. Also a loading path is shown to influence the guest uptake efficiency. Due to this the host may undergo an expansion or local distortion. Typical examples are the hydrogen sorption by metals , or the intercalation of $`Li`$ ions into layered materials . Quite often the host undergoes structural transformations  due to intercalation. These are the structural transitions. In such cases the standard LG approach also operates with several sublattices  corresponding to each configuration of the host. However, such a restructuring suggests that the elastic effects should be taken into account. A general thermo-mechanical theory of the stress-composition interaction is developed by Larche and Cahn . This implies the existence of a coupling between the elastic properties of the host material and the structure and/or the dynamics associated with the guest species. In other words the configurational and structural transitions should be considered to be coupled. In our previous papers we have investigated such a coupling for two-dimensional  as well as for three-dimensional systems. The main advantage of our approach is its capability in describing various (microscopically different) intercalation systems within a common framework. Despite of the fact that some microscopic details are omitted, the theory agrees well with experimental data on layered $`Li_xTiS_2`$ and crystalline $`Li_xWO_3`$ and $`Na_xWO_3`$ compounds at equilibrium conditions.

In this study we focus on a correlation between the equilibrium properties (isotherms and differential capacity) and the concentration dependence of the chemical diffusion coefficient. Our purpose is to determine the impact of the thermodynamics to the kinetic properties. Following our previous works we combine the LG model, describing the configurational transitions and the linear elasticity theory that accounts for the loading mechanism and the volume dilatation in the course of the structural transitions (expansion or restructuring). This allows us to derive an effective chemical potential, involving the concentration dependent stress and strain fields. Then, assuming that the diffusion flux is proportional to the gradient of the chemical potential, we investigate the behavior of the chemical diffusion coefficient.

## II Model

### A Microscopic formulation

The host material is described as a three-dimensional lattice of adsorbing sites with their positions given by a set vectors $`(𝐫_i)`$. Due to the elastic properties of the real host each site (for instance, an interstitial site) of this ”auxiliary” lattice may deviate from its equilibrium position $`𝐫_i^0`$, such that we deal with the displacements $`𝐮_i=𝐫_i𝐫_i^0`$. Therefore, the host properties are described by the Hamiltonian $`H_H[\{𝐮_i\}]`$. Note however that a connection between the elastic properties of the real matrix and those of the adsorbing lattice, is not straightforward.

A distribution of the intercalants on the host sites is given by a set of occupation numbers $`\{t_i\}`$, with $`t_i=0`$ or $`t_i=1`$. The guest subsystem is characterized by the chemical potential $`\mu `$ and the nearest neighbor interaction parameter $`W`$ (at equilibrium positions). For the electrochemical insertion, chemical potential $`\mu `$ gives a deviation of the electrode potential $`eV`$ from its standard value $`E_0`$. In this way an arrangement of the intercalants on a rigid host lattice is governed by the LG Hamiltonian

$$H_G=W\underset{ij}{}t_it_j\mu \underset{i}{}t_i$$
(1)
that describes the configurational transitions of the intercalated species. These could be the droplet formation for attractive interactions ($`W<0`$) or the order-disorder transition for repulsive interactions ($`W>0`$). In the latter case one introduces the sublattice concentrations as appropriate for the symmetry of a given system.

The coupling between the host and the guest is given by Hamiltonian $`H_C[\{𝐮_i\},\{t_i\}]`$ which takes into account a dependence of the binding energy on the site displacement and also the pairwise interaction between the guest particles through the host lattice. The overall Hamiltonian is now written as

$$H=H_H[\{𝐮_i\}]+H_G[\{t_i\}]+H_C[\{𝐮_i\},\{t_i\}]$$
(2)
The free energy $`F`$ corresponding to the above Hamiltonian is given by

$$F=F_H+F_G(x)+F_C(x)$$
(3)
where $`x`$ is the intercalant concentration. Here $`F_H`$ is the host free energy in the absence of intercalation, $`F_G`$ is the guest free energy in a case of the rigid host lattice. The latter term can be calculated using the mean field approximation to obtain the well-known relation for the chemical potential

$$\mu _0(x)=eVE_0=cWx+\frac{1}{\beta }\mathrm{ln}\left(\frac{x}{1x}\right)$$
(4)
where $`c`$ is the coordination number of the host lattice and $`\beta =1/kT`$. The coupling term $`\beta F_C=\mathrm{ln}(e^{\beta H_C}_{(𝐮_i)}_{(t_i)})`$ requires the averaging over the displacements and the occupation numbers, calculated with the reference terms $`H_H`$ and $`H_G`$. In fact, this is an infinite series including the correlations of all orders in the reference state. The main problem is to specify $`H_C[\{𝐮_i\},\{t_i\}]`$ coherently with the host symmetry and elastic properties. It is known that real host materials have rather complicated elastic properties (e.g. a strong anisotropy). Therefore, only some simplified model calculations are expected to give tractable results. On the other hand, such predictions (e.g. a rigid plane model) do not exhibit quantitative agreement with experimental data.

### B Approximation

The perturbative scheme developed previously  has shown that the coupling term is concerned with a concentration dependence of the host response to the intercalation. This involves at least two effects. First is a renormalization of the net pair interaction between the intercalants due to their indirect interaction through the matrix. The second is a change of the host volume upon insertion of the guest species. Depending on the host nature, a stress field may result if the lattice is not totally free to relax. In this situation it seems reasonable to estimate the coupling term based on the continuum theory of elasticity with the concentration dependent stress and strain fields. Then the host-guest free energy is a sum of the lattice gas and elastic part

$$F(x)=F_{LG}(x)+F_{el}(x)$$
(5)
where $`F_{LG}(x)`$ is the configurational (lattice gas) part, in which the pairwise interaction is renormalized due to the interaction through the matrix. This gives the chemical potential $`\mu _0(x)`$ (see eq. (4)) with a new interaction constant $`W`$. The elastic part is approximated by the free energy of a strained isotropic body under a loading stress $`\sigma (x)`$. Since the strain is assumed to be purely dilatational, we operate with traces $`\epsilon `$ and $`\sigma `$ of the corresponding tensors.

$$F_{el}(x)=\frac{\mathrm{\Lambda }}{2}\epsilon (x)^2\sigma (x)\epsilon (x)$$
(6)
with $`\mathrm{\Lambda }`$ being the effective elastic constant, independent of the concentration. The total stress $`S=S(x)`$ is given by

$$S(x)=\frac{dF_{el}}{d\epsilon (x)}=\mathrm{\Lambda }\epsilon (x)\sigma (x)$$
(7)
Therefore, we have two stress contributions. The internal, or self-stress $`\mathrm{\Lambda }\epsilon (x)`$ corresponds to the host reaction to the guest insertion. The second term $`\sigma (x)`$ describes a loading procedure, that may include the sample clamping or other effects which are not directly related to the strain. It is important that, in general, $`\sigma (x)`$ is a function of $`x`$ (not a function of $`\epsilon (x)`$). For instance, if the sample is clamped such that $`\epsilon (x)=0`$, then we have a stress accumulation proportional to the concentration $`\sigma (x)x`$.

## III Equilbrium properties

The guest chemical potential $`\mu (x)`$ is given by the concentration derivative of the total free energy.

$$\mu (x)=\mu _0(x)+S(x)\frac{d\epsilon (x)}{dx}\frac{d\sigma (x)}{dx}\epsilon (x)$$
(8)
Here the second term involves the so-called chemical expansion coefficient $`d\epsilon (x)/dx`$, while the last term is associated with the loading path. It is seen that the intercalation level depends on an interplay of the internal stress and the loading stress. The latter could be small, but its concentration derivative is not necessary small, so that the loading path may induce serious consequences. In particular, for a given material ($`\mathrm{\Lambda }`$) and a suitable loading path $`\sigma (x)`$, there may be a cancellation of the last two terms in eq. (8) in a given domain of $`x`$. This explains why in some cases the purely configurational description ($`\mu (x)=\mu _0(x)`$) works well.

For comparison to other theoretical approaches it is convenient to express the strain function in terms of other relevant quantities. The strain can be measured as a volume dilatation or as a change of the interlayer spacing in the course of insertion. If $`V(x)`$ is a composition dependent sample volume (or interlayer spacing for layered compounds), then by definition

$$\epsilon (x)=\frac{V(x)V(0)}{V(0)}=\delta p(x)$$
(9)
where $`\delta =[V(1)V(0)]/V(0)`$ is the relative volume variation and $`p(x)=[V(x)V(0)]/[V(1)V(0)]`$ is a modulating function. The latter varies in between $`0`$ and $`1`$. If $`p(x)`$ is a linear function, then the system is said to obey the Vegard’s law.

If the host-guest system forms a solid solution, then the guest partial molar volume $`V_m(x)`$ is related to the total volume

$$V(x)=xV_m(x)+V(0),$$
(10)
where $`V(0)`$ is the initial host volume (at $`x=0`$). Then we can find the following relation for the strain, the volume variation and partial molar volume

$$\epsilon (x)=\frac{V(x)V(0)}{V(0)}=x\frac{V_m(x)}{V(0)}$$
(11)
If the loading is composition independent $`\sigma (x)=\sigma `$ and the sample volume changes linearly ($`V_m=const,\epsilon (x)=xV_m/V(0)`$, like in the $`PdH_x`$ $`\alpha `$-phase), then we recover the well-known result

$$\mu (x)=\mu _0(x)\left(\sigma \mathrm{\Lambda }x\frac{V_m}{V(0)}\right)\frac{V_m}{V(0)}$$
(12)
The term liner in the concentration can be viewed as an additional interaction which can be combined with $`cWx`$ in $`\mu _0(x)`$. Although this additional term is proportional to the elastic constant $`\mathrm{\Lambda }`$, it is physically different from the interaction of elastic dipoles which is already absorbed into $`W`$ (see the discussion after equation (5)). This interaction takes place even if the sample is perfectly clamped ($`\epsilon (x)=0`$), while our term represents a cooperative effect due to $`\epsilon (x)0`$.

In order to analyze the role of elastic effects in the phase behavior of the guest species it is instructive to consider a simple example of vanishing loading stress $`\sigma (x)=0`$, $`S(x)=\mathrm{\Lambda }\epsilon (x)`$. In this case there is no difference between the volume expansion or contraction. Then we arrive at

$$\mu (x)=\mu _0(x)+\overline{\mathrm{\Lambda }}p(x)\frac{dp(x)}{dx}$$
(13)
where $`\overline{\mathrm{\Lambda }}=\mathrm{\Lambda }\delta ^2`$. In addition, $`p(x)`$ is taken to be that in the so-called layer rigidity model which describes a deviation from the Vegard’s law in layered compounds

$$p(x)=1(1x)^q$$
(14)
where $`q`$ is related to the rigidity of the host layers. The first order phase transition (of liquid-gas type) is manifested by a singularity in the differential capacity $`C=dx/d(\beta \mu )`$

$$\beta C=\left[\frac{d\mu _0(x)}{dx}+\overline{\mathrm{\Lambda }}\left(\frac{dp(x)}{dx}\right)^2+\overline{\mathrm{\Lambda }}\frac{d^2p(x)}{dx^2}\right]^1$$
(15)
It is clear that for $`q<1`$ both derivatives of $`p(x)`$ are positive and the only possibility for the critical behavior is the case of attractive interactions $`W<0`$, i.e. when the configurational part

$$\beta \frac{d\mu _0(x)}{dx}=\frac{1+\beta cWx(1x)}{x(1x)}$$
(16)
becomes negative. Then the last two terms in the denominator of (15) rescale the critical temperature and shift the critical composition from $`x=1/2`$. For $`q>1`$ the term $`d^2p/dx^2`$ becomes negative and the criticality may appear due to the elastic effects even if $`W=0`$. The coexistence curves for this case, obtained by a numerical analysis of the divergence in $`C`$, are plotted in Fig. 1. It is seen that the critical temperature (the maximum) increases with increasing $`q`$. Simultaneously the critical concentration decreases. The elastic effects break the hole-particle symmetry of the problem and consequently the curves are not symmetric and are not centered around $`x=1/2`$.

In our previous works  we have analyzed another extreme case of vanishing total stress $`S(x)0`$ and the loading stress proportional to the concentration $`\sigma (x)=\sigma _0x`$. Then

$$\mu (x)=\mu _0(x)\gamma p(x)$$
(17)
where $`\gamma =\delta \sigma _0`$. Note, however that, for $`\mathrm{\Lambda }=const`$, the above result is valid only for a linear behavior of $`p(x)=px`$, while is an approximation for an arbitrary $`p(x)`$. In reality the elastic constants could depend on the concentration, for instance the bulk modulus of $`Pd`$ is reduced up to $`20\%`$ due to the hydrogen sorption. Thus one can easily imagine a situation when $`\mathrm{\Lambda }(x)\epsilon (x)\sigma (x)`$ and then the approximation (17) is applicable for a nonlinear $`p(x)`$. In general, for a non-Vegard’s behavior and the linear loading $`\sigma (x)=\sigma _0x`$ we deal with

$$\mu (x)=\mu _0(x)+\left[\overline{\mathrm{\Lambda }}p(x)\gamma x\right]\frac{dp(x)}{dx}\gamma p(x)$$
(18)
The strain can be measured as a volume dilatation or as a change of the interlayer spacing during the transitions between different phases (staging in graphite, restructuring in $`Li_xWO_3`$, $`\alpha \beta `$ transition in $`PdH_x`$, etc). In any case the lattice parameters do not obey the linear Vegard’s law , exhibiting a sharp change near the transition compositions $`x_n^0`$. This can be described by the following approximation.

$$p(x)=\frac{1}{2}\left[1+\underset{n}{}p_n\mathrm{tanh}[\alpha _n(xx_n^0)]\right]$$
(19)
where $`n`$ counts the number of phase boundaries. In $`Li_xWO_3`$, $`n=1,2`$ corresponds to monoclinic-tetragonal and tetragonal-cubic transitions, respectively. For $`Li`$-graphite systems, $`n`$ corresponds to the stage transfer boundaries. The set of rigidity parameters $`\alpha _n`$ controls a local slope and $`p_n`$ are the weights corresponding to each phase, such that $`_np_n=1`$. This is consistent with experimental observations indicating that the structures are not completely pure, but contain some features that indicate a mixing of phases. The present form of $`p(x)`$ corresponds to continuous structural transitions, but can be easily modified to take into account the jump-like behavior (like for the staging ). Nevertheless, in that case the concentration derivatives of $`p(x)`$ would be singular, inducing the singularities in the thermodynamic quantities (as it should be in the neighborhood of a phase transition). Note however that the experimental dependencies are usually smoothed due to a finite concentration resolution. Then it is difficult to distinguish between the step-wise variation and a sharp (but continuous) transition. For that reason the criticality criteria, determined on the ground of rigorous statistical mechanical arguments ($`\mu (x)`$ loops or a divergent capacity $`C(x)`$), only approximately conform to experimental data. Therefore we stay with the continuous $`p(x)`$, performing the fitting to the experiment.

Following this methodology, in Fig. 3 we display the intercalation isotherm (inset) and the capacity curve, comparing them to the experimental data for $`Li_xWO_3`$. The fitting (eq. (18)) is performed assuming that the effective pair interaction inside the matrix is repulsive for any $`x`$. The curve modulation results from the strain behavior $`p(x)`$, associated with the volume dilatation. As is discussed previously, the peculiarities (inflection points and the peaks in $`C(x)=dx/d(\beta \mu )`$) mark the boundary between different host symmetries. From the experimental point of view, the peak height is usually associated with a transition sharpness. Nevertheless, as we have discussed above, the critical behavior of the model itself should be studied separately.

## IV Chemical diffusion coefficient

In many cases there are two mobile species inside host matrices: the host electrons and the guest particles. Therefore, the transport is characterized by their flux densities $`𝐉_e`$ and $`𝐉_G`$, related to the gradients of the corresponding chemical (or electro-chemical) potentials $`\mu _e`$ and $`\mu _G`$.

$$𝐉_e=L_{ee}\mu _e+L_{eG}\mu _G$$
(20)
$$𝐉_G=L_{GG}\mu _G+L_{Ge}\mu _e$$
(21)
where $`L_{ab}`$ are phenomenological transport coefficients, which are scalar quantities if the host is spatially isotropic. In matrices with a metallic conductivity the electrons are much more mobile than the guest species. Due to this the guest ionic charge is compensated by the electrons. For the same reason the electrons reach a uniform (equilibrium) distribution much faster than the intercalants ( $`\mu _e=0`$ on the time scale when $`\mu _G0`$). In addition, usually $`L_{GG}>>L_{eG}`$ . Therefore one may focus on the guest flux

$$𝐉=M(x)x\mu =D(x)x$$
(22)
where $`M(x)`$ is the mobility coefficient that must contain a blocking factor $`M(x)=D_0(1x)`$, determining the concentration-dependent chemical (or collective) diffusion coefficient $`D(x)`$. A comprehensive review of theoretical approaches to the determination of the chemical diffusion coefficient can be found in Ref. . In the simplest approximation we deal with

$$D(x)=D_0x(1x)\frac{d\beta \mu }{dx}=D_0x(1x)/C(x)$$
(23)
Under the assumptions above, the intercalant transport is described by the diffusion equation with an effective diffusion coefficient $`D(x)`$. The latter requires an information on the intercalation isotherm $`\mu (x)`$. Then, based on the equilibrium properties, at a given concentration gradient, one, at least in principle, can solve the kinetic problem.

Note, that such a simple scheme should be modified when the host conductivity changes remarkably upon the intercalation. Then the electronic impact should be taken into account. Also, we do not consider other driving forces, like external field or stress gradients, assuming that the guest concentration, $`x`$, is the only independent variable. And finally, near the phase coexistence (e.g, staging) a system becomes non-uniform because of the phase boundaries. Then the formulation of the kinetic problem must be coherent with the theory of critical phenomena (see  for a recent review).

Starting from (8) and (23) we obtain

$$D(x)=D_0(x)+D_{el}(x)$$
(24)
where the first term is the standard LG part $`D_0(x)=D_0x(1x)\frac{\mu _0}{x}`$ and the elastic contribution is given by

$`D_{el}=B(x)\left[S{\displaystyle \frac{d^2\epsilon }{dx^2}}+\mathrm{\Lambda }\left({\displaystyle \frac{d\epsilon }{dx}}\right)^22{\displaystyle \frac{d\epsilon }{dx}}{\displaystyle \frac{d\sigma }{dx}}{\displaystyle \frac{d^2\sigma }{dx^2}}\epsilon \right]`$ (25)
with $`B(x)=D_0x(1x)`$. The diffusion coefficient involves several competing factors (strains, stresses and their concentration derivatives). In general it is not trivial to see whether $`D`$ increases or decreases with $`x`$. The situation is even more complicated for systems in which the lattice spacing does not obey the linear Vegard’s law, but has inflection points separating different phases (e.g. staging in graphite or restructuring in $`Li_xWO_3`$). Then the derivatives above may change sign with the concentration. Also we see again that the loading path has a significant contribution.

For $`PdH_x`$ ($`\alpha `$-phase) we recover the well known result . The concentration induced internal stress $`S_0(x)=\mathrm{\Lambda }xV_m/V_0`$ increases the diffusion

$$D=D_0(x)+B(x)\mathrm{\Lambda }\left(\frac{V_m}{V_0}\right)^2,$$
(26)
where $`D_0(x)`$ is the stress free contribution (it corresponds to $`\mu _0(x)`$, i.e. the lattice gas description, see above). Note that the non-local stress effects are not discussed here.

Since, according to (23), the diffusion coefficient is just an inverse of $`C(x)`$, $`D`$ as a function of $`x`$ has minima (Fig. 3), corresponding to the peaks in the capacity curve. For $`WO_3`$ (Fig. 3) the theory works well near the first minimum, but overestimates the diffusion at higher $`x`$. Assuming that eq.(23) is valid, we have also inverted the capacity data (Fig. 2) in order to demonstrate that there is no ”hidden” errors in the fit for $`C(x)`$. Since the equilibrium characteristics (Fig. 2) are predicted with a reasonable accuracy, then one of the source for the discrepancy could be our simplified assumptions on the mobility. For instance, one might suppose that eq.(23) does not work for higher $`x`$ because the hoping $`D_0`$ depends on the elastic properties. By the analogy with the two-dimensional diffusion, $`D_0`$ is related to the lattice spacing. The latter decreases with $`x`$ for $`WO_3`$, so that the resulting $`D(x)`$ would also decrease. Nevertheless the volume change does not exceed several percents. This cannot explain the two orders of magnitude difference in Fig. 3. The lattice anisotropy seems also to be irrelevant because the final (high $`x`$) structure is cubic, so the isotropic approximation for the diffusion and elastic properties is reasonable.

It is known that $`Li_xWO_3`$ changes its electronic conductivity near the structural instability composition. Then the ionic transport must correlate with the transport of the neutralizing electrons. Similar effects occur in superconductors. Therefore, the electronic mobility impact (that is absent from our approach) to the current should be taken into account by including the additional driving force $`\mu _e=d\mu _e/dxx`$. This contribution can be included into our approach, assuming some model for the electronic structure, e.g. - the rigid band model. This would correct our predictions for small $`x`$ (according to , for $`x0.1`$), but our high-$`x`$ estimation would remain unchanged. We believe that the electronic effects are not responsible for the decrease of $`D(x)`$ at high concentrations. Based on the fact that the electrical properties of $`Li_xWO_3`$ become more and more metallic with increasing $`Li`$ content, we may expect that our approximation $`\mu _e=0`$ becomes more reliable with increasing $`x`$. Moreover, if the conductivity effects were important, then the ionic diffusion would increase like this occurs in amorphous $`WO_3`$ films.

On our opinion we deal with an interplay of several facts. Although the model works well in predicting the equilibrium properties, nevertheless it is probably too simple for the description of the kinetics. In fact, the transport is assumed to take place in an infinite lattice through the hoping of neutral species ($`Li^++e`$), ignoring the kinetics other relevant phenomena, like the charge transfer (metal/host, guest/host), the exclusion (permselectivity) effect, the formation of the passive layer at the host/electrolyte boundary, etc. On the other hand, the experimental data on the diffusion are extracted from direct measurements employing a theoretical model, which could differ (in some details) from the one introduced here. Concerning the order of magnitude estimations, it should be noted that the experimental data for the same substance are usually differ significantly depending on the sample preparation, its size and the scanning time interval. Therefore, the quantitative description of the ionic transport is a delicate problem requiring mutual theoretical and experimental efforts. Nevertheless, our estimation of the elastic effects is qualitatively correct in predicting that the system characteristics changes coherently with the strain $`p(x)`$, and therefore, positions of the minima (Fig. 3) $`x_n^0`$ are predicted correctly, implying a correlation between the $`D(x)`$ behavior and the strain development.

Similar situation takes place for $`Li`$-graphite. In Fig. 4 the normalized interlayer spacing for graphite in the presence of intercalated $`Li`$ is plotted as a function of the guest composition. The fitting is preformed using eq. (19) under a suitable choice of the parameters $`p_n`$, $`\alpha _n`$, $`x_n^0`$. The experimental data are taken from Refs. . Our fit is rather reasonable, ”catching” the fact that the stage transfer compositions $`x_n^0`$ approximately correspond to the inflection points in $`p(x)`$, where the average interlayer spacing changes sharply. Note that we do not discuss the fine structure of the phase diagram, such as dilute and liquid-like phases for the same stage. These features are related to the the guest in-plane ordering which is indistinguishable in the behavior of the interlayer spacing.

Having a reliable approximation for $`p(x)`$ we calculate the diffusion coefficient in the framework outlined above. The fitting to the experimental data  is shown in Fig. 5. It is seen that the diffusion slows down significantly near the stage transfer compositions $`x_n^0`$. Although the magnitude of $`D(x)`$ is close to the experimental results, the width of the minima is underestimated. Note that the fitting is not optimized, that is, we did not try to find an optimal set of the parameters, which gives equally good agreement for the isotherm $`\mu (x)`$ and $`D(x)`$. Nevertheless, as in the case of $`Li_xWO_3`$ (Fig. 3), there is a clear correlation between the diffusivity and the strain behavior.

## V Conclusion

In summary, our approach implies that the main equilibrium and transport features of the intercalation systems (which differ in their microscopic details) can be well understood in terms of a concentration dependence of the hydrostatic parts of the stress and strain fields, associated with the internal and loading effects. It is shown that the elastic effects may induce the critical behavior even if the ”direct” interaction between the guest species is absent. Although the ionic transport is supposed to be equivalent to the diffusion of neutral (ion plus electron) species, the approach is flexible enough to incorporate other transport mechanism (migration, or electronic mobility effects).

The theory gives a quantitative description of different insertion processes, involving the volume expansion ($`Li_xTiS_2`$ (see) , $`\alpha PdH_x`$), staging ($`Li`$-graphite) or restructuring ($`Li_xWO_3`$). For all these processes the theory implies a correlation between the intercalation isotherm and a concentration dependence of the diffusion coefficient. The latter exhibit a set of characteristic minima, related to the boundaries between different phases (like different symmetry phases of $`Li_xWO_3`$, $`\alpha \beta `$ transition in $`PdH_x`$, etc). However the experimental dependencies are usually smoothed due to a finite concentration resolution. For that reason the criticality criteria, determined on the ground of rigorous statistical mechanical arguments (e.g. $`D(x)=0`$ at the transition concentrations) only approximately conform to the experimental data, exhibiting a sharp (but finite) decrease of $`D`$.

For a non-Vegard’s strain variation the diffusion coefficient is a nonlinear (and nonmonotonic) function of the concentration (25). Therefore the diffusion equation for the concentration profile would be strongly nonlinear. Then one can expect a rather complicated space-time variation, including , for instance, oscillations and other nonlinear effects.

Our results may have implication in various domains related to the insertion process, like hydrogen sorption, electrochemical intercalation, impurities in alloys, layered superconductors, volume transitions in hydrated gels, etc.


text:
# 1 Introduction

## 1 Introduction

In a series of papers  Sen has proposed ways to describe the classical time evolution of the open string tachyon field $`T(t)`$ on an unstable D-brane, at strictly zero string coupling $`g_s=0`$. Some of the results obtained using the BCFT method  were also reproduced from an effective field theory  and boundary string field theory . Since the closed string degrees of freedom, which are responsible for the dissipation of the energy liberated from the D-brane, are not taken into account there, the energy must be stored in the open string sector in spite of the fact that the D-brane is decaying. In order for the energy to be conserved in the settings of , it turns out that the ‘velocity’ of the tachyon field must approach a finite value $`\dot{T}(t)1`$ in the asymptotic future $`t\mathrm{}`$. This is possible only when the tachyon potential has its minimum at infinity. Then, what happens in theories where the minimum of the tachyon potential is at a finite point in the field configuration space, such as open string field theory and $`p`$-adic string theory? The authors who have tried to seek rolling solutions in such theories have found the oscillations with ever-growing amplitude, typically of the form $`T(t)=_na_ne^{nt}`$ . It is possible that these seemingly different behaviors of the tachyon field are related to each other through some field redefinition, but at present we do not know precisely how to reconcile them. Furthermore, it has been suggested that, once the non-zero string coupling $`g_s`$ is turned on, the energy which was initially stored in the open string sector is converted to the closed string radiation, no matter how small $`g_s`$ is .

In this paper, we will try to attack the problem of rolling tachyon from the viewpoint of open-closed string theory. For simplicity we will restrict ourselves to the spatially homogeneous decay of an unstable space-filling D-brane. And we assume that the tachyon vacuum is at a finite point in the field configuration space, in which case much more definite discussion is possible than in the case of the runaway potential. Then, we expect the following properties for the open string tachyon condensation in open-closed string theory:

1. There are two homogeneous and static (i.e. constant in spacetime) solutions corresponding to the unstable D-brane background and the open string tachyon vacuum where the D-brane has disappeared. These two configurations should have the same energy density because the energy initially stored in the D-brane (=the open string modes) is just transferred to the closed string modes.
2. The rolling of the tachyon is described by a time-dependent solution which interpolates between the above two vacuum configurations. In particular, since the open string sector can lose its energy, it would be natural to consider that the tachyon field will eventually relax in the (finite) tachyon vacuum.

As is clear from the appearance of a factor of $`\sqrt{\mathrm{}}`$ in the open-closed string action (3), we must go beyond the classical level to incorporate the interactions between open and closed strings.222For a recent discussion about the quantum aspects of Witten’s cubic open string field theory, see . Since it is quite difficult to deal with the full open-closed string field theory (see  and references therein), here we consider a simple model including tachyons only, as a first step toward the ultimate goal. In the next section we introduce the model and show that it admits a rolling solution satisfying the properties mentioned above. Section 3 is devoted to summary and discussion.

## 2 A Toy Model for Open and Closed String Tachyons

Let us consider a model described by the following action,333We are working in units where $`\alpha ^{}=1`$.

$`S`$ $`=`$ $`{\displaystyle d^Dx}={\displaystyle d^Dx(_o+_c+_{\mathrm{int}})};`$
$`_o`$ $`=`$ $`{\displaystyle \frac{1}{2}}\varphi \mathrm{}\varphi +{\displaystyle \frac{1}{2}}\varphi ^2,`$
$`_c`$ $`=`$ $`{\displaystyle \frac{1}{2}}\psi \mathrm{}\psi +2\psi ^2,`$
$`_{\mathrm{int}}`$ $`=`$ $`{\displaystyle \frac{1}{3}}\stackrel{~}{\varphi }^3+c_2\stackrel{~}{\varphi }\stackrel{~}{\psi }\stackrel{~}{\varphi }^2\stackrel{~}{\psi },`$ (1)
where $`\varphi `$ and $`\psi `$ are the open and closed string tachyon fields respectively, and we have defined

$$\stackrel{~}{A}(x)e^{(\mathrm{log}K)\mathrm{}}A(x)=K^{\mathrm{}}A(x),$$
(2)
with $`K=2`$. We take the flat spacetime metric to be $`\eta ^{\mu \nu }=\mathrm{diag}(1,1,\mathrm{},1)`$, so that the d’Alembertian is $`\mathrm{}=_t^2+^2`$. We do not need to specify the spacetime dimensionality $`D`$. $`c_2`$ is a constant which will be fixed below.

One may notice that the structure of the action (1), in particular the $`\stackrel{~}{}`$-operation (2), is reminiscent of the open-closed string field theory action

$`S_{\mathrm{OCSFT}}\left(\mathrm{\Phi }|Q_B|\mathrm{\Phi }+g_o\mathrm{\Phi }^3_{\mathrm{disk}}+g_o^2\mathrm{\Phi }^4_{\mathrm{disk}}\right)+\mathrm{\Psi }|(c_0\overline{c}_0)(Q_B+\overline{Q}_B)|\mathrm{\Psi }`$
$`+\sqrt{\mathrm{}}\left(g_c\mathrm{\Psi }^3_{\mathrm{sphere}}+\mathrm{\Psi }_{\mathrm{disk}}+g_o\mathrm{\Phi }\mathrm{\Psi }_{\mathrm{disk}}+g_o^2\mathrm{\Phi }^2\mathrm{\Psi }_{\mathrm{disk}}\right)+𝒪(g_o^3,\mathrm{}),`$ (3)
where $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ are open and closed string fields, respectively. By truncating the string fields as

$$|\mathrm{\Phi }=\frac{d^Dk}{(2\pi )^D}\varphi (k)c_1|k,|\mathrm{\Psi }=\frac{d^Dk}{(2\pi )^D}\psi (k)c_1\overline{c}_1|k,$$
(4)
and substituting them into the action (3), we would obtain the lagrangian whose structure is similar to (1).444We have approximately calculated the open-closed transition vertex with the stub and the strip omitted, and found that it is roughly estimated as $`\mathrm{\Phi }\mathrm{\Psi }_{\mathrm{disk}}\stackrel{ˇ}{\varphi }(x)\stackrel{ˇ}{\psi }(x)`$ with $`\stackrel{ˇ}{A}(x)\mathrm{exp}((\mathrm{log}1.352)\mathrm{})A(x)`$. It is quite difficult to compute the open-open-closed vertex $`\mathrm{\Phi }^2\mathrm{\Psi }_{\mathrm{disk}}`$ because it includes the integration over the moduli space. However, there should be some discrepancies between the resulting action and (1). In particular, there is no reason why $`\psi `$, $`\psi ^3`$ and $`\varphi ^4`$ terms can be neglected. Therefore we consider the action (1) as just a toy model, instead of having been derived from the open-closed string field theory action via level truncation.

The vacuum structure of this model can be studied by looking at the potential

$$𝒱=|_{\varphi ,\psi =\mathrm{const}.}=\frac{1}{2}\varphi ^22\psi ^2+\frac{1}{3}\varphi ^3c_2\varphi \psi +\varphi ^2\psi .$$
(5)
If we eliminate the closed string tachyon field by its equation of motion $`\psi =\frac{1}{4}(\varphi ^2c_2\varphi )`$, we get the effective potential for $`\varphi `$,

$$𝒱_{\mathrm{eff}}=\frac{1}{8}\varphi ^2\left\{\left(\varphi +\frac{4}{3}c_2\right)^2+\frac{8}{3}c_2\frac{52}{9}\right\}.$$
(6)
From this expression, we find that, if we choose $`c_2=\frac{13}{6}`$, the effective potential takes the form $`𝒱_{\mathrm{eff}}=\frac{\varphi ^2}{8}\left(\varphi \frac{5}{6}\right)^2`$, so that $`𝒱_{\mathrm{eff}}`$ has two degenerate vacua at $`\varphi =0`$ and $`\frac{5}{6}`$. With this value of $`c_2`$, the potential (5) has three stationary points,

$$\begin{array}{cccc}\text{solution (I)}& :& (\varphi ^\mathrm{I},\psi ^\mathrm{I})=(0,0),\hfill & 𝒱=0,\hfill \\ \text{solution (II)}& :& (\varphi ^{\mathrm{II}},\psi ^{\mathrm{II}})=(\frac{5}{6},\frac{5}{18}),\hfill & 𝒱=0,\hfill \\ \text{solution (III)}& :& (\varphi ^{\mathrm{III}},\psi ^{\mathrm{III}})=(\frac{5}{12},\frac{35}{192}),\hfill & 𝒱=\frac{625}{165888},\hfill \end{array}$$
(7)
where the value of $`𝒱`$ for each solution shows the height of the potential there.

Next we will explain that the solutions (I) and (II) can be regarded as representing the D-brane and the open string tachyon vacuum, respectively. First, note that there is a $`\varphi `$-$`\psi `$ mixing term in the action (1) at the quadratic level. In order to determine the perturbative spectrum around any one of the solutions, we must diagonalize it. Let us start with the solution (I). After the Fourier transformation, the quadratic part of the action (1) can be arranged as

$`S_{\mathrm{quad}}^\mathrm{I}={\displaystyle \frac{1}{2}}{\displaystyle \frac{d^Dk}{(2\pi )^D}(\varphi (k),\psi (k))^\mathrm{I}(k^2)\left(\begin{array}{c}\varphi (k)\\ \psi (k)\end{array}\right)},`$ (10)
$`^\mathrm{I}(k^2)=\left(\begin{array}{cc}k^21& \frac{13}{6}K^{2k^2}\\ \frac{13}{6}K^{2k^2}& k^24\end{array}\right).`$ (13)
The mass spectrum is found by looking for the values of $`k^2=m^2`$ at which the eigenvalues of the matrix $`^\mathrm{I}(k^2)`$ vanish. This can be equivalently accomplished by solving $`det^\mathrm{I}(k^2)=0`$. Since we cannot solve this equation analytically, we resort to the numerical study. We see from Fig. 1, where $`det^\mathrm{I}`$ is shown as a function of $`k^2`$, that there is a closed string tachyon state with $`m_c^24.000`$, and an open string tachyon state with $`m_o^20.863`$.

We therefore consider the solution (I) as the unstable D-brane background. We do not concern ourselves about an extra state around $`m^2=0.110`$ which is not important for our purpose.

Now we turn to the solution (II). When we expand the fields as

$$\varphi =\varphi ^{\mathrm{II}}+\varphi ^{},\psi =\psi ^{\mathrm{II}}+\psi ^{},$$

the action quadratic in the fluctuation fields $`\varphi ^{},\psi ^{}`$ becomes

$`S_{\mathrm{quad}}^{\mathrm{II}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{d^Dk}{(2\pi )^D}(\varphi ^{}(k),\psi ^{}(k))^{\mathrm{II}}(k^2)\left(\begin{array}{c}\varphi ^{}(k)\\ \psi ^{}(k)\end{array}\right)},`$ (16)
$`^{\mathrm{II}}(k^2)=\left(\begin{array}{cc}k^21+(2\varphi ^{\mathrm{II}}+2\psi ^{\mathrm{II}})K^{2k^2}& (2\varphi ^{\mathrm{II}}\frac{13}{6})K^{2k^2}\\ (2\varphi ^{\mathrm{II}}\frac{13}{6})K^{2k^2}& k^24\end{array}\right).`$ (19)
The determinant of the matrix $`^{\mathrm{II}}(k^2)`$ is plotted in Fig. 2.

We see that $`det^{\mathrm{II}}`$ vanishes only once near $`m_c^24.000`$. This means that the perturbative spectrum around the solution (II) contains the closed string tachyon state, but not the open string tachyon. Since there is no physical open string state, we consider that the solution (II) corresponds to the background with no D-brane, i.e. the open string tachyon vacuum.

Now that we have obtained two vacuum configurations which can be interpreted as a D-brane and the tachyon vacuum, we look for a time-dependent solution interpolating between them. The equations of motion following from the action (1) are

$`(_t^2+1)K^{2_t^2}\stackrel{~}{\varphi }(t)\stackrel{~}{\varphi }(t)^2+{\displaystyle \frac{13}{6}}\stackrel{~}{\psi }(t)2\stackrel{~}{\varphi }(t)\stackrel{~}{\psi }(t)=0,`$
$`(_t^2+4)K^{2_t^2}\stackrel{~}{\psi }(t)+{\displaystyle \frac{13}{6}}\stackrel{~}{\varphi }(t)\stackrel{~}{\varphi }(t)^2=0,`$ (20)
where we have let $`\varphi `$ and $`\psi `$ be functions only of the time variable $`t`$. The differential operators appearing in (20) can be rewritten as the convolution form

$`(_t^2+\mu )K^{2_t^2}A(t)=𝒞^{(\mu )}[A](t)`$
$`{\displaystyle \frac{1}{\sqrt{8\pi \mathrm{log}K}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑s\left({\displaystyle \frac{(ts)^2}{(4\mathrm{log}K)^2}}+{\displaystyle \frac{1}{4\mathrm{log}K}}+\mu \right)e^{\frac{(ts)^2}{8\mathrm{log}K}}A(s).`$ (21)
Then we can numerically solve the ‘integral equations’ (20) using the iterative procedure: Given a pair of the $`n`$-th functions $`\varphi _n(t)`$ and $`\psi _n(t)`$, we can calculate the $`(n+1)`$-th functions as

$`\stackrel{~}{\varphi }_{n+1}(t)={\displaystyle \frac{6}{13}}\left(\stackrel{~}{\varphi }_n(t)^2𝒞^{(4)}[\stackrel{~}{\psi }_n](t)\right),`$
$`\stackrel{~}{\psi }_{n+1}(t)={\displaystyle \frac{6}{13}}\left(𝒞^{(1)}[\stackrel{~}{\varphi }_n](t)+\stackrel{~}{\varphi }_n(t)^2+2\stackrel{~}{\varphi }_n(t)\stackrel{~}{\psi }_n(t)\right).`$ (22)
If this iterative procedure converges, the set of functions $`\stackrel{~}{\varphi }_{\mathrm{}}(t),\stackrel{~}{\psi }_{\mathrm{}}(t)`$ gives a solution to the equations of motion (20). Starting with the initial configurations

$`\stackrel{~}{\varphi }_0(t)=\{\begin{array}{ccc}\varphi ^\mathrm{I}=0\hfill & & (t0)\\ \varphi ^{\mathrm{II}}=\frac{5}{6}\hfill & & (t>0)\end{array},`$ (25)
$`\stackrel{~}{\psi }_0(t)=\{\begin{array}{ccc}\psi ^\mathrm{I}=0\hfill & & (t0)\\ \psi ^{\mathrm{II}}=\frac{5}{18}\hfill & & (t>0)\end{array},`$ (28)
we have repeated the iterative procedure (22) 7000 times, and the resulting profiles for $`\stackrel{~}{\varphi }_{7000}(t)`$ and $`\stackrel{~}{\psi }_{7000}(t)`$ are shown in Fig. 3.

We have verified that $`\stackrel{~}{\varphi }_{7000}(t)`$ and $`\stackrel{~}{\varphi }_{6000}(t)`$ (and similarly for $`\stackrel{~}{\psi }_{7000}(t)`$ and $`\stackrel{~}{\psi }_{6000}(t)`$) completely agree with each other at least up to significant six digits for every value of $`t`$. We consider this as strong evidence that $`\stackrel{~}{\varphi }_n(t)`$ and $`\stackrel{~}{\psi }_n(t)`$ really converge in the limit $`n\mathrm{}`$. Note that, since $`\varphi =\stackrel{~}{\varphi }`$ holds within the range where $`\varphi `$ is a constant, the tachyon fields $`\varphi ,\psi `$ indeed approach their vacuum values with vanishing velocity in this solution.

We will now check whether the energy is conserved in time for the solution obtained above. The energy $`E(t)`$ is calculated by the following formula

$`E(t)=E_o(t)+E_c(t)+E_{\mathrm{int}}(t);`$
$`E_o(t)=_o+{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{2l1}{}}}(1)^m\left({\displaystyle \frac{}{\varphi _{2l}}}\right)_m\varphi _{2lm},`$
$`E_c(t)=_c+{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{2l1}{}}}(1)^m\left({\displaystyle \frac{}{\psi _{2l}}}\right)_m\psi _{2lm},`$
$`E_{\mathrm{int}}(t)=_{\mathrm{int}},`$ (29)
where the subscripts denote the number of time-derivatives, $`A_n(t)\frac{^n}{t^n}A(t)`$. Although this formula contains infinite sums, we can obtain a good approximation to it by ignoring higher derivative terms, because they have negligibly small orders of magnitude as compared to $`\varphi `$ and $`\psi `$ themselves,555This also shows that the profile of $`\stackrel{~}{\varphi }`$ ($`\stackrel{~}{\psi }`$) is almost the same as that of $`\varphi `$ ($`\psi `$). as shown in Fig. 4.

Here we keep only up to the fourth derivatives. Note that it is meaningful to compare $`\varphi `$ (or $`\psi `$) with its derivatives since the time variable $`t`$ has been made dimensionless by setting $`\alpha ^{}=1`$. The total energy $`E(t)`$ calculated this way is illustrated in Fig. 5A.

If we compare it with Fig. 5B where the constituents $`E_o,E_c,E_{\mathrm{int}}`$ defined in (29) are shown separately, it is clear that highly non-trivial cancellations do occur among $`E_o,E_c`$ and $`E_{\mathrm{int}}`$. We thus conclude that the total energy is constant ($`=0`$) in time. However, as pointed out in , we do not know how to separate the total energy into the open and the closed string sectors. It would be interesting if we can see the flow of the energy from open to closed strings.

We have also investigated whether we can construct a time-dependent solution which interpolates between the solutions (I) and (III) using the iterative procedure, starting with

$`\stackrel{~}{\varphi }_0(t)=\{\begin{array}{ccc}0& & (t0)\\ \varphi ^{\mathrm{III}}=\frac{5}{12}& & (t>0)\end{array},`$ (32)
$`\stackrel{~}{\psi }_0(t)=\{\begin{array}{ccc}0& & (t0)\\ \psi ^{\mathrm{III}}=\frac{35}{192}& & (t>0)\end{array}.`$ (35)
From the result shown in Fig. 6,

we see that the tachyon profiles tend to $`\varphi (t)\psi (t)0`$ as $`n`$ becomes large. Such a behavior is natural because there should be no energy-conserving solution connecting these two vacua in time.

## 3 Summary and Discussion

In this paper we have argued that in the full open-closed string theory, the spatially homogeneous decay process of a space-filling unstable D-brane should be described by a time-dependent solution which interpolates between two static and energetically degenerate field configurations representing the unstable D-brane and the open string tachyon vacuum. In particular, to realize the situation in which the tachyon field converges to some finite value, it was essential to incorporate the bulk closed string modes because the open string sector must dissipate the energy which amounts to the D-brane tension through the couplings to closed strings. We have constructed and investigated a simple model describing the open and closed string tachyons, and found that (i) this model has two degenerate vacua which can be interpreted as the unstable D-brane and the tachyon vacuum, (ii) there exists a time-dependent solution which interpolates in time between them, and (iii) the total energy is conserved during the decay process described by this time-dependent solution. Hence it seems that some features expected of the rolling tachyon in open-closed string theory are indeed captured by our model.

There are many things to be studied further. We list some of them below.

* Among others, our analysis based on the toy model (1) is far from complete. In particular, there is no plausible justification for neglecting various higher-order interactions and for truncating the string fields to the tachyons. It would be interesting if we find that the time-dependent solution obtained in section 2 persists to exist after reinforcing our model by including higher-order interactions and other massless or massive fields.
* We have focused on the spatially homogeneous decay of a space-filling D-brane. In the case of the inhomogeneous decay or the decay of a lower-dimensional D-brane, the closed strings emitted from the decaying D-brane will be propagating outward in the transverse space. In this sense, the asymptotic state of the closed string field is not static. Nevertheless, we do expect that the open string field reaches the tachyon vacuum configuration sometime and ceases to evolve in time thereafter.
* We have identified the solution (II) with the open string tachyon vacuum by looking at the fluctuation spectrum around it. It is desirable to see that the energy stored in the open string field is lower at the solution (II) than at the solution (I). As mentioned in section 2, however, it is not known how to split the total energy into open and closed string sectors in our tachyon model. It may become possible if we can incorporate the metric degrees of freedom in our model, as in .
* It seems that the convolution form of the equations of motion is difficult to relate to the ‘initial value problem’ . Our numerical solution corresponds to the rolling from the top of the potential with vanishing initial velocity, and has no parameter to be chosen freely. In fact, we have encountered a similar situation in seeking for static lump or kink solutions in (super)string field theory , where we did not treat it as a boundary value problem, and the width of the lump (kink) seemed to be determined automatically.

We hope that our investigation presented here sheds new light on the problem of the time-dependent decay of unstable D-brane systems.

## Acknowledgments

I am grateful to T. Eguchi, I. Kishimoto, N. Moeller, M. Schnabl and E. Watnabe for valuable discussions. I would also like to thank M. Fujii, Y. Matsuo, R. Nobuyama, K. Sakai, Y. Tachikawa and H. Takayanagi for useful conversations. The numerical computations in this work were performed using Mathematica. This work is supported by JSPS Research Fellowships for Young Scientists.